Nonlinear Processes in Geophysics The effect of volatile bubble growth rate on the periodic dynamics of shallow volcanic systems

Many volcanic eruptions exhibit periodic behavior. For instance, periodic ground inflations and deflations in proximity to a volcano are the consequences of periodic overpressure variations in the magma conduit and periodic magma flow rate. The period varies from a few hours to many years, depending on the volcano parameters. On the other hand, volatile components exsolve from an ascending magma by forming bubbles. The strong dependence of the melt viscosity with the volatile concentration generates a positive feedback on the magma flow. We consider here the effect of the growth of volatile bubbles on the dynamics of a magmatic flow in a shallow volcanic system. Various expressions for the bubble growth rate are treated, thus generalizing previous work. In particular, a growth rate law derived from a recent many-bubble theory is considered. It is seen that, for a range of flow rate values at the base of the magma conduit, the system undergoes a Hopf bifurcation. Periodic solutions compatible with the observations are generated. This work shows that measurements of volcanic activity have the potential to test various bubble growth models in magmatic systems.


Introduction
Volcanic eruptions are the consequences of a succession of complex processes acting over a wide range of timescales and understanding their dynamics is an important, but challenging task.One example of intriguing volcanic eruptions dynamics consists in the occurrence of an oscillatory behavior.For instance, tilt meters installed on the crater of the Soufrière volcano (Montserrat Island) in 1996 and 1997 detected periodic ground inflations that correlated with periodic Correspondence to: I. L'Heureux (ilheureu@uottawa.ca)seismic activity and periodic magma flow activity (Voight et al., 1998(Voight et al., , 1999;;Wylie et al., 1999).The period was of the order of 10-20 h.Other volcanic systems also exhibit cyclic episodes.Examples are found in Mount St. Helens (Washington State, USA) between 1980 and 1986 (with an average period of 74 days or 230 days, depending on the cycle) (Barmin et al., 2002), Santiaguito (Guatemala) from 1922 to the present (with an average period of 10.7 years) (Harris et al., 2003), Mount Unzen (Japan), Karymsky (Russia) and Merapi (Java) (Nakada et al., 1999;Ozerov et al., 2003;Voight et al., 2000).
On the other hand, various volatile species are commonly found in magmatic systems (mainly H 2 O, CO 2 and sulfide species) and volatile bubble dynamics is an important contributing factor to the understanding of volcanic eruptions (Sparks, 1978;Bottinga and Javoy, 1990;De Vivo et al., 2005;Behrens and Gaillard, 2006;Proussevitch andSahagian, 1996, 1998).The object of this contribution is to examine in details one simple nonlinear dynamical model in which an oscillatory behavior due to the growth of volatile bubbles occurs in a shallow volcanic system; the effects of various bubble growth rate laws are also investigated.The treatment is limited here to the typically dominant volatile species, water (De Vivo et al., 2005;Behrens and Gaillard, 2006).
Similar models are available in the literature.In Wylie et al. (1999) (thereafter referred to as WVW), the authors presented a constant temperature model that is mostly applicable to shallow magmatic systems (such as Montserrat).In that model, the magma flow dynamics is coupled to changes in melt viscosity induced by variations in the dissolved volatile content (specifically water).In turn, through the formation and growth of bubbles, volatile exsolution controls the volatile concentration dynamics.
Melt viscosity depends not only on the volatile concentration but also on the magma temperature and its chemical composition (Hess and Dingwell, 1996).In Costa and I. L'Heureux: Effect of volatile bubble growth rate on periodic dynamics of shallow volcanic systems Macedonio (2002), a model is presented in which temperature variations (rather than volatile content) induce changes in the melt viscosity, that feedback, in turn, on flow dynamics.No dynamics are presented by these authors, but the possibility of a bursting effect (whereby the flow rate changes abruptly back and forth between two values) is described. Mourtada-Bonnefoi et al. (1999) presents a compartment model in which crystal magma content and temperature are coupled through two ordinary differential equations.Again, bubble bursting effects are described.However, the model ignores transport processes.In Barmin et al. (2002) and Melnik andSparks (1999, 2005), a more detailed transport model is described that takes into account the viscosity changes induced by variations in temperature and volatile content in a volcanic conduit of circular cross section.A crystal growth factor is also included in that model as the magma crystal content affects the permeability of the system and the flow dynamics.This model was further generalized (Costa et al., 2007) to the situation where the lower part of the conduit has an elliptical cross section.
Other mechanisms that do not involve directly the presence of volatile have been proposed to explain oscillatory eruptive behaviors.Ozerov et al. (2003) considered the magma as a visco-elastic fluid that undergoes stick-slip transitions in a boundary layer adjacent to the conduit walls.They showed that an oscillatory behavior with short periods (seconds to minutes) may ensue.
In this paper, we generalize the WVW model in four ways: 1.In that model, it was assumed that the volatile concentration field adjusts quickly to its steady state value.Strictly speaking, this is valid only in the limit where the magma pressure relaxes very slowly in response to changes in magma flow rate.This is not true in general and we will not make this assumption here.
2. WVW introduced an approximation to simplify the mathematical treatment of their model: it was assumed that the integral of the viscosity has the same functional dependence as the integral of the volatile concentration field.As the relation between viscosity and volatile concentration is not linear, this is not a correct assumption in general.We do not use this simplification here.
3. WVW uses a simple linear bubble growth law without taking into account the degree of supersaturation.Here, we consider also more realistic bubble growth laws: the single bubble growth rate law in the diffusion-limit case (Navon and Lyakhovsky, 1998) and a recently published (L'Heureux, 2007) many-bubble growth law that takes into account (in the mean field sense) the collective effect of the other bubbles on the growth of a given bubble.
4. Although a linear stability is presented by Nakanishi and Koyaguchi (2008) for the original Barmin et al.'s model (2002), no detailed linear stability analysis is available for the WVW model, generalized to an arbitrary bubble growth rate law.Such an analysis and a dynamical phase diagram in parameter space are presented here.
The theory presented here is applicable to shallow magmatic systems (such as Montserrat) as it considers only the effect of volatile content on the melt viscosity and neglects temperature changes during the magma short ascent.Indeed, the time scale associated with conducting cooling of a shallow ascending magma via heat exchange through the walls can be estimated from the magma thermal diffusivity (about 10 −6 m 2 /s, Costa and Macedonio, 2002) and the conduit lateral dimension (about 10 m); its value is very large compared to the relevant period of the oscillations (Voight et al., 1999).
Over the time scales of interest here, we can thus neglect the cooling effect.
For simplicity, we also neglect the viscosity change due to increase in the magma crystal content.Previous models (Melnik and Sparks, 2005;Barmin et al., 2002) have parameterized the viscosity dependence on crystal volume fraction in the following way: for a crystal volume fraction smaller than a critical value of about 0.7, the viscosity is basically constant; otherwise it increases rapidly with the crystal volume fraction.For simplicity, we assume here that the crystal volume fraction is limited to a value smaller than 0.7, so that the viscosity does not depend much on crystallinity.The same assumption was used by Wylie et al. (1999).A more complete theory would also include both crystallization kinetics and temperature variations (through conduction, heat exchange at the conduit walls and latent heat production).However, we feel that investigating the effects of various volatile bubble growth laws in more details is a necessary first step in developing an understanding of periodic dynamics in shallow magmatic systems.Indeed, we plan to investigate the effects of parametric random noise driving the model system through its boundary conditions and a detailed study of the corresponding deterministic system -as presented in this contribution -is needed first.
The paper is structured as follows.In Sect.2, we present the details of the model and how it generalizes the WVW model.A linear stability analysis is presented in Sect. 3 for the three bubble growth rate models considered here.In Sect.4, we comment on the numerical method used to solve the model, and we present and discuss some numerical results.A comparison with the tilt angle measurement data of Fig. 1 of Wylie et al. (1999) is made using various bubble growth rate expressions.We then make some concluding statements in Sect. 5. Finally, Appendix A presents a list of the model variables and their corresponding notation and Appendix B estimates a limiting value of the cycle period.

Model
As in many previous models of periodic volcanic activity (Wylie et al., 1999;Barmin et al., 2002;Melnik andSparks, 1999, 2005), we describe the volcanic system as a thin vertical conduit in which magma flows from a deep magma reservoir to the surface.It is in the upper part of the conduit that the relevant degassing dynamics occurs.It is thus convenient to divide the conduit into two sections.The lower section of the conduit connects the deep magma reservoir to the upper conduit section.No significant degassing occurs in the lower conduit.The upper conduit is characterized by a length L and a constant radius r.Observations of the Soufrière shallow volcanic conduit are indeed compatible with a simple cylindrical geometry, as resulting from previous explosive eruptions (Sparks and Young, 2002).
We call Q(z,t) the magma flow rate (m 3 /s) in the upper conduit, p(z,t) the total magma pressure and c(z,t) the dissolved volatile concentration (in units of mass of dissolved volatile reported per melt mass).Here z is the spatial coordinate measured upwards from the base of the upper conduit (z=0) and t is time.The melt viscosity in the upper conduit is denoted η(c) and depends on the magma's dissolved volatile content.
A few words qualitatively describing the mechanism that generates oscillatory behavior are in order.The lowest section of the conduit connects the deep magma reservoir to the upper conduit section.No significant degassing occurs in the lower conduit and the volatile concentration is kept constant at a value c o .We will assume that the magma flow rate from the reservoir is constant and has a value Q o .We also define the magma overpressure in the lower conduit P as the actual fluid pressure, from which are subtracted the hydrostatic pressure and the atmospheric pressure p a .It is also assumed that P (t) is spatially homogeneous in the lower conduit, although it depends on time in general.Now suppose that the flow rate in the upper conduit is smaller than the value Q o at its base.Then, the lower conduit overpressure will increase, which will increase the pressure in the upper conduit and will increase the upper conduit flow rate.The dissolved concentration will then increase in the upper conduit through this enhancement in advective input.At some time however, bubble growth will occur as a result of the increasing volatile supersaturation and, as a result of this exsolution, the concentration of dissolved volatile will decrease.However, melt viscosity is a strongly dependent function of the concentration of dissolved volatile: as the concentration decreases, viscosity increases rapidly.This, in turn, slows down the flow rate back to its initial value, and the cycle has the potential to start over again, thus generating repeated cycles.

Basic equations
The total magma density is ρ g φ + ρ(1 − φ) where ρ g is the volatile gas density, ρ is the melt density (assumed constant) and φ is the vesicularity (the fraction of the magma volume occupied by the gas).We also neglect the motion of the gas bubbles relatively to the ascending melt.The magma continuity equation in the upper conduit reads: Rapid degassing often leads to magma fractionation (whereby the magma becomes a foam containing a dispersion of crystals and melt drops) and to explosive eruptions, in which the gaseous mixture can reach supersonic velocities (Mader, 1998).Here, we assume that degassing is sufficiently slow, so that the magma mixture does not reach the fragmentation stage and the flow remains subsonic.Thus, as was done in previous models of periodic volcanic activity (Wylie et al., 1999;Barmin et al., 2002;Melnik andSparks, 1999, 2005), the magma can be considered incompressible.Equation ( 1) then states that the flow rate Q(t) depends on time only.The Navier-Stokes equation for magma in the upper conduit can be written for a flow in a cylindrical pipe.The ratio of the inertial term to the viscous term scales as the Reynolds number Re∼ ρQ/Lη, which is typically small (∼10 −4 ).The inertial term can therefore be neglected and one can use Poiseuille law to describe the magma laminar flow.The pressure gradient can be written in terms of the flow rate as: where g is the acceleration of gravity.The boundary conditions are such that the pressure must match the total pressure at the bottom of the upper conduit and it must be equal to the atmospheric pressure at the upper conduit outlet: The overpressure P in the lower conduit varies in response to the difference between the input flow rate Q o and the upper conduit flow rate (Landau and Lifshitz, 1970): Here, Y is the average Young's modulus of the surrounding rock, σ its Poisson ratio and V the volume of the deep magma reservoir and the lower conduit.
The volatile mass conservation equation is now described.The volatile is present in the upper conduit in two possible phases: a dissolved phase (of concentration c) and a bubble phase.Let R denotes a typical bubble radius and N the bubble number density (per melt volume).The total volatile concentration m (volatile mass reported per melt mass) is Here, the ideal gas law has been used to express the volatile gas density, T o being temperature, R the ideal gas constant and M the volatile molar mass.Surface tension effects between the bubble and the surrounding fluid have been neglected, so that the gas pressure has been set equal to the local fluid pressure p.The mass conservation equation for total volatile takes the form: We have neglected here volatile losses through the sides of the conduit.Combining Eq. ( 7) with Eq. ( 1) (with ρ g φ neglected compared to ρ(1 − φ)) gives: Similarly, the continuity equation for the bubble number is: where J is the bubble nucleation rate reported per melt volume.For the shallow systems of interest here, we assume that bubble dynamics is in a post-nucleation regime: growth of bubbles that have previously been nucleated in the lower conduit dominates over the formation of new bubbles.We will thus set J =0.This can be justified since the surface tension between water vapor and silicic melt increases substantially as the pressure decreases (Proussevitch et al., 1993), so that bubble nucleation occurs preferably at lower depths.
Finally, the kinetics of the volatile dissolved phase is described by where G is a term proportional to the bubble growth rate (reported per melt mass), as detailed below.The boundary condition results from matching the concentration of dissolved volatile with its value in the lower conduit: Various bubble growth models lead to different expressions for G.We will consider here three models: 1. WVW model: In Wylie et al's model, G is taken as a linear function of the concentration, for simplicity: where k is a rate constant (in s −1 ).
2. Single bubble growth (SBG) model: From the well known diffusion-limited single bubble growth theory (Navon and Lyakhovsky, 1998), we learn that the radial bubble growth rate v (m/s) is proportional to the volatile diffusion coefficient D and to the supersaturation (c−c eq ), where c eq is the value of the concentration in thermodynamic equilibrium with the fluid: Multiplying this by 4π R 2 ρ g gives the bubble mass increase per unit time.Multiplying again by N/ρ gives the relative mass change in dissolved volatile per unit time over the N bubbles per unit volume.Thus An expression for the equilibrium concentration is found from a generalized Henry law: where K H is the Henry constant and n is an empirical constant that depends on the identity of the volatile.For water, n= 1 / 2 to a good approximation (Sparks, 1978;Burnham, 1975).The fact that c eq is proportional to the square root of the pressure in the appropriate pressure range is indicative of the presence of two water species in a silicate melt: molecular H 2 O and hydroxyl OH associated with the silicate framework (Behrens and Gaillard, 2006).
3. Multiple-bubble growth (MBG) model: In the single bubble growth theory, the competitive effects of the other bubbles are neglected.In L'Heureux (2007), a model was proposed that approximately relaxes this constraint.There, the growth of randomly located bubbles is treated in the framework of a mean field approximation that generalizes the approach of Marqusee and Ross (1984) to leading order in the bubble volume fraction.One of the approximate consequences of this treatment in is that the steam bubble growth rate decays exponentially with time under isobaric conditions: where the rate constant K is given by Here, c eq,0 is the initial equilibrium concentration at a given position and R ∞ is the large-time limit of the bubble radius, which is related to m and p via Eq.( 6):  16), ( 17) have been derived under constant pressure conditions.In general, the pressure changes in time at a given position.Nevertheless, for simplicity, we will assume a pseudo-steady state approximation for the bubble growth, so that we can use the instantaneous pressure in the growth rate expression.Finally, the melt viscosity is needed.Empirical expressions for the melt viscosity as a function of temperature, water content and magma composition have been used regularly by igneous petrologists (see for instance Hess and Dingwell, 1996).These expressions are based on a generalization of the Arrhenius law for viscosity in the form log where b 1 ,b 2 and b 3 are functions of the water content c.We use here a simple linear expression for log η (Clemens and Petford, 1999;Shaw, 1972), which is a good approximation for c between 2% and 5% and for temperatures around 900 • C: Here, β defines a temperature-dependent viscosity response coefficient and η o is the melt viscosity at the base of the lower conduit, where c = c o =5%.In fact, the widely used empirical expression proposed by Hess and Dingwell (1996) reduces to Eq. ( 19) when c is not too far from c o .

Dimensionless formulation
It is convenient to express the model in a dimensionless form.
We scale the position coordinate by the length L of the upper conduit and time by a scale t, to be determined later.We now introduce the scaled concentrations c and m , the scaled bubble number density N , the scaled viscosity η , the scaled flow rate Q and the scaled pressures p and P as: Here, N o is the value of N at the base of the upper conduit.
The relevant equations then take the form: Here, are dimensionless parameters representing the melt density, the atmospheric pressure, the magma reservoir elastic response and the reservoir flow rate, respectively.The time scale t depends on the adopted bubble growth model and will be selected so as to simplify the scaled bubble growth rate G =Gt/c o as much as possible.Explicitly, one obtains: For the WVW model: For the SBG model (using Eq. 6 to eliminate R): where For the MBG model, t has the same value as in the SBG model.Adopting n = 1 / 2 in Henry's law, the scaled growth rate becomes: In summary, the problem is defined by Eq. ( 21) for the variables p , c , m , N , P and Q , with the parameters α, β, γ , δ, ε,Q o and A. Unless stated otherwise, we will drop the symbols.

Simplified formulation
It is possible to rewrite the model as an integro-differential system involving the variables P , Q and c only and which will be more amenable to a numerical solution.We now introduce a time-like variable τ as (Barmin et al., 2002): I. L'Heureux: Effect of volatile bubble growth rate on periodic dynamics of shallow volcanic systems The equations for the dissolved volatile concentration and for the overpressure become, respectively: The pressure equation can be integrated to in which the boundary condition at z=1 imposes the following integral constraint on the overpressure: Finally, a change of variable ξ = z − τ , τ * = τ transforms Eq. (21c) to This indicates that the total volatile concentration m is constant and equal to its initial value.Similarly, the scaled bubble number density can be taken equal to unity.
In their model, WVW introduce two simplifying assumptions, besides taking G = c.(i) By integrating Eq. ( 28) over z from 0 to 1, the quantity c(1,τ ) appears.In order to estimate it, WVW assume that the concentration profile adjusts itself quickly to its stationary value, so that c(1,τ ) = exp(−1/Q(τ )).(ii) Instead of Eq. ( 19), they make the approximation which is correct only for small β.Eliminating P , they then obtain a simple system of two coupled ordinary differential equations for the variables Q and 1 0 c(z ,τ )dz .We will not make these simplifications here, but consider the full system (21).

General formulation
It is straightforward to solve the system for its steady state (Q s , c s (z), P s , p s (z)).Setting the derivatives with respect to τ equal to zero gives: c s can be solved, at least numerically.For example, for the WVW growth rate, G = c and c s (z) = exp(−z/Q o ).
We now perform a linear stability analysis.For this purpose, we simplify the pressure dependence in the growth rate by using its hydrostatic value in the argument of G: The numerical results performed on the exact formulation (Sect.4) confirm that this approximation leaves unchanged the nature of the dynamical instability and does not change the stability phase diagram substantially.We then set where δc(z), δP and δQ are small perturbations and a complex frequency, to be determined.Also, the boundary condition c s (0)=1 implies that δc(0)=0.Substituting in Eqs. ( 29) and ( 31), it is easy to eliminate δP : where the integrated viscosity is Equation ( 36) requires the derivative of the viscosity with respect to the volatile concentration.With the relation (21e), this is: The approximation (34) indicates that G does not depend on overpressure but only on the concentration profile and on z, through p s = α(1 − z) + δ.We let G s be the bubble growth rate evaluated with the stationary concentration profile.After linearization, Eq. ( 28) becomes in turn: where The solution of Eq. ( 39) with the boundary condition δc(0)=0 reads: where It is convenient to simplify the function G as follows.The total differential of the stationary growth rate G s is: where use has been made of dc s /dz=−G s /Q o .Thus, in Eq. ( 42): In consequence, Finally, using Eqs.(41) in Eq. ( 36), and dividing by δQ, we obtain the general frequency equation, which reduces the linear stability problem to solving the following transcendental equation for : Since ε > 0, it is not possible for the stationary solution to lose its stability ( = 0) with a real frequency.A Hopf bifurcation is characterized by a loss of stability with Re( )=0, Im( ) = 0 and is the signature of the emergence of a limit cycle (periodic) solution, which itself could be stable or unstable.In order to find the locus of Hopf bifurcations in parameter space {ε, Q o , β}, we set ≡ i in Eq. ( 46) with real.Separating the real and imaginary parts and using Eq. ( 45), one finds For a given value of Q o and β, Eq. ( 47) is used to solve for (which can be done numerically using a Newton-Raphson solver coupled to an integration code).Then, Eq. ( 48) gives directly the Hopf bifurcation curve ε H = ε(Q o , β).

Examples
We now present Hopf bifurcation curves for the WVW, SBG and MBG growth rate laws.For the simplest case (WVW), we have The stationary overpressure P s reduces to: where Ei(x) ≡ − ∞ −x dt e −t /t is the exponential integral.Figure 1a plots P s as a function of Q o for two values of β.For sufficiently high value of β, P s exhibits two increasing branches connected by a decreasing branch.Let β* be the critical value of β above which these three branches exist.β* is found from the condition dP s /dQ o = d 2 P s /dQ 2 o = 0 to be β*=2.989.
The frequency equation is taken from Eq. ( 46) with ∂G s /∂p s = 0.It reads: where (a,x) ≡ ∞ x z a−1 e −t dt is the incomplete Gamma function.From this expression, it is straightforward to obtain the Hopf bifurcation curve ε H , an example of which is illustrated in Fig. 1b for β=5.5.The vertical dotted lines correspond to the two values In fact, it is easy to show analytically from P s = Q o µ and from Eq. (47) (in the limit → 0) that the boundaries of the Hopf bifurcation curve are bounded by Q o *.This is true for an arbitrary growth rate G of the form (34) and it shows an important property: a necessary condition for a Hopf bifurcation to exist is to have β > β* and to choose a value of Q o for which the corresponding P s (Q o ) lies on a decreasing branch dP s /dQ o <0.However, in contrast to what was suggested in Wylie et al. (1999), this not a sufficient condition: ε must be sufficiently large (ε > ε H ) for an instability to occur.
For the two other growth models (SBG and MBG), expressions for P s (Q o ) and ε H are not analytically available.However, a numerical approach is straightforward.Figures 2  and 3 show the stationary overpressure and an example of  a Hopf bifurcation curve for the two growth models, respectively.The choice of parameter values is discussed in the next section.For this choice, the critical values of β* are 2.60 for the SBG model and 3.75 for the MBG model.Although the details are quantitatively different from the WVW case, the topology of the curves is identical.In summary, the stationary state is stable when β < β*.For β > β*, there exists values of Q o for which P s decreases as a function of Q o .The stationary state is stable when Q o is chosen in a range for which dP s /dQ o > 0. Otherwise, it is stable when ε is sufficiently small.In all case, loss of stability derives from a Hopf bifurcation.

Numerical approach
In order to investigate the nature of the attractor in the unstable regime, numerical solutions of the system (28, 29, 31) were generated.The overpressure Eq. ( 29) is solved by a second-order finite difference scheme in the time-like variable τ of step size h: Table 1.Parameter values used in the calculations of Sect.4.2.These are typical of the Montserrat system (Wylie et al., 1999).The value of β is estimated from the temperature and from the empirical relation ( 19) (Shaw, 1972).The value of Henry constant K H , the volatile diffusion constant D and the initial bubble radius R(0) are taken from Proussevitch et al. (1993).The bubble number density N is inferred from Proussevitch et al. (1993).

Parameter
Units Value where P n and Q n denote respectively the overpressure and the flow rate at τ = nh with n an integer.The advection equation (Eq.28) is solved using an upstream explicit scheme.
The step sizes h in z-and τ -space are chosen to be equal, so as to minimize numerical dispersion.Let c n i and p n i denote the concentration and pressure fields at position z=ih and τ =nh where i is an integer.This algorithm gives 1/Q n+1 is chosen in a self-consistent manner (via a Newton-Raphson algorithm) so as to be consistent with the discretized version of the integral constraint (31): where w i are weight constants associated with the particular numerical integration method used (with h small enough, a trapezoidal integration method was found to be sufficient).With the field c updated, the full pressure field can be finally found by the numerical integration of Eq. ( 30).The algorithm is fast, simple and convergent.

Numerical solutions
In this section, we have left the parameters Q o and ε free but have fixed the other parameter values according to Table 1.  1.The initial condition has been chosen as Q(0)=1 with an initial pressure and concentration profiles that have the steady state form for this value of Q(0).
With these parameter values, the time scale t is 1 h for the WVW model and 1.58 h for the SBG and the MBG models.The scaled parameters are α=29.1,γ =14.2, δ=0.3 and A=0.046.In the growth rate expressions, we have also used the approximation of Eq. ( 34) in order to verify the bifurcation analysis performed in Sect.3. Figure 4a illustrates for the MBG model an overpressure and flow rate time series in a regime where dP s /dQ o < 0 but for a case where the steady state is stable (ε < ε H ). Figure 4b illustrates the corresponding time series when ε is just above ε H .In agreement with the analysis, the steady state is indeed unstable in the latter case.The solution evolves towards a stable limit cycle.Figure 4c illustrates a limit cycle obtained for a larger value of ε > ε H .This illustrates that the pressure cycle approaches a saw-tooth time-behavior and that the flow rate exhibits bursting.These features are typical and are similarly found in the two other growth models, WVW and SBG (not illustrated).
The cycle period T (with respect to the time variable t) is a highly sensitive function of Q o and ε.This is illustrated in Fig. 5 for the MBG growth model where each curve shows T as a function of ε for various values of Q o .The curves T (ε) are similar for the two other growth models.In general, T decreases as Q o increases for a fixed value of ε.But for fixed Q o , the period is an increasing function of ε.The lower limit of each curve corresponds to the period at the Hopf bifurcation point 2π/(Q o ) where is the Hopf frequency found from Eq. (47).At the other limit (for large ε), it is  1.
easy to show that the period takes the asymptotic value (see Appendix B): where P s (Q) is the steady state overpressure curve and the labels A, B, C, D correspond to the boundary of the limit cycle in {P , Q} space, as shown, for instance, in Fig. 1a.The fact that that the period curves obtained from the numerical solutions go smoothly to its analytical asymptotic expression is a further indication of the adequacy of the numerical algorithm.

Comparison with observations: the Soufrière system
In this section, we present numerical results for the case where the full pressure expression is kept in the growth rate, so that the approximation (34) is relaxed.The general features of the solution are qualitatively identical as those described above.We will use this approach to see how the model applies to actual observations.For concreteness, we will consider the tilt angle measurements performed on the shallow Soufrière system in August 1997 (top of Fig. 5 in Voight et al., 1998).When the steady state is unstable, the asymptotic state of our model is a limit cycle with a welldefined period and amplitude.However, the observed periods and oscillation amplitudes are not quite constant.This could be due to a combination of transient effects and/or systematic slow or random variations in the system parameters which are not described by our simple model.Nevertheless, by requesting that the limit cycle solution is close to a typical measured cycle, our model can be used to obtain reasonable estimates of two of the parameter values that are not easily accessible: the time scale t and the dimensionless magma chamber elastic parameter ε.
For the purpose of illustration, we selected the particularly well-defined cycle from 2 August to 3 August, which has a period T obs =18 h and is characterized by a dimensionless asymmetry parameter a obs ≡ (t max − t min )/T obs =0.6 where t max and t min are the times at the cycle maximum or minimum, respectively.Both T obs and a obs do not depend on the specific (possibly non linear) conversion between tilt angle measurements and overpressure values.We only need to assume that extrema in tilt angles occur at the same time as extrema in overpressure.From the limit cycle solutions P (t), we have adjusted t and ε by minimizing the residual relative error e where T calc is the dimensionless computed period and a calc the calculated asymmetry parameter.We have taken (Wylie et al., 1999) Q o =10 m 3 /s.For the MBG model, the initial pressure profile was taken as the steady state one.Table 2 gives the values of t and ε that minimize e for the three growth models.As an example, Fig. 6 illustrates the overpressure and flow rate cycles for the MBG model for these  2).The origin of the time axis is arbitrary as the transients are not illustrated.
values of t and ε.The curves are qualitatively similar for the two other growth models.Nevertheless, the minimal residual error is about seven times larger for the SBG growth model than for the two other models.Taking the estimates (Wylie et al., 1999) Y =3×10 10 Pa and σ =0.2 with L, and r from Table 1, we obtain the product V η o from Eq. ( 22).This estimate is reported in Table 2.For comparison, the estimate V η o =2.34×10 5 km 3 Pa-s is inferred from Wylie et al. (1999), with a time scale arbitrarily chosen as t=1 h.Table 2 also gives the ratio (in s −1 ) of the overpressure drop over a cycle, P , by the viscosity η o .For η o =10 6 Pa-s (Voight et al., 1999), the overpressure drop is of the order of the MPa, consistent with the observations.The table also gives the flow rate amplitude Q in m 3 /s.Notwithstanding the difference in the values of t (which implies a difference in the pressure and flow rate scales), the amplitude of the variations P and Q are comparable for all three bubble growth rate models.
Assuming a linear scaling between the overpressure P at the base of the conduit and tilt angle measurements (Voight et al., 1999), the following chi-square can be estimated.Here θ i are measured tilt angles (from the 2-3 August cycle in Fig. 5 of Voight et al., 1998), B and C are two constants defining the assumed linear relation between the calculated (dimensionless) overpressure P i and the tilt angle.σ i (estimated to be 10% of θ i ) is the error on the tilt angle measurement and f the number of degrees of freedom.
We have selected the parameters t and ε from Table 2.The values of χ are reported in Table 2.These values suggest that the SBG model does not provide a good fit, whereas both the MBG and WVW growth models lead to a marginally good fit (with a slightly better fit in favor of the MBG model).Nevertheless, on physical grounds, one expects the bubble growth expression to have a more complex form than what Eq. ( 11) implies.Notwithstanding the large number of simplifying assumptions, the application of this model to data allows one to constraint the values of the time scale t and of the elastic response parameter ε, which are not easily available otherwise.

Conclusion
In this contribution, we have extended the model of Wylie et al. (1999) describing a simple mechanism for the generation of an oscillatory behavior in a shallow magmatic system.In this model, magma flow is coupled to the magma water content through the explicit volatile dependence of the melt viscosity.We have extended the original model by relaxing three of its simplifying assumptions: use is made of an arbitrary volatile bubble growth rate law expression, we do not assume that the volatile concentration profile has a steady-state form and the spatial integral of the viscosity in Eq. ( 31) is properly evaluated.We have also performed a linear stability analysis of the system and have established that the oscillatory behavior occurs via a Hopf bifurcation.The range of parameter values for which oscillatory behavior exists have been found for three different bubble growth models: the linear growth model of Wylie et al., the single bubble diffusion-limited growth model and an approximate many-bubble growth model that has been recently published (L'Heureux, 2007).We have found that this modeling approach can generate oscillatory magma dynamics with a period that is compatible with the observations.The quality of the fit with the observations can be used to constraint some system parameters that are not easily measured.Since we assumed that only variations in the volatile concentration affect the melt viscosity, the present model is more suitable to shallow magma reservoirs.But other factors (such as temperature variations and changes in magma chemical composition due to mineral growth) may also play an important role in changing the viscosity.In addition, crystal precipitation is expected to modify the porosity of the fluid system and to change its flow characteristics.Some previously published models (Barmin et al., 2002;Melnik and Sparks, 1999;Melnik and Sparks, 2005) have considered these effects but not with a many-bubble growth model.Such considerations are the object of our continuing research.It will also be interesting to investigate how random parameter fluctuations influence the oscillatory dynamics of a volcanic system (Wiesenfeld, 1985).This is part of our planned work.
I. L'Heureux: Effect of volatile bubble growth rate on periodic dynamics of shallow volcanic systems from Q o , the pressure cycle is described by a slow motion along the ascending branches of steady state diagrams such as those of Figs.1a, 2a or 3a.The overpressure slow dynamics is described by Eq. ( 29): Thus it is seen that the pressure can not assume values on the decreasing branch of the steady state diagram(where dP s /dQ o < 0).Rather, the motions along the two ascending branches are connected by horizontal lines joining B to C and D to A in Fig. 1a, resulting in very fast variations in Q(τ ).To calculate the cycle period T , we rewrite Eq. (B1) as εdP /dt = Q o − Q where we used dτ /dt = Q (Eq.27).
The period is then: where we neglected the time used to run along the horizontal paths BC and DA.Integrating by parts with respect to Q with P = P s = Qµ s finally gives Eq. ( 54).

Fig. 1 .
Fig. 1.(a) Steady state overpressure P s as a function of the flow rate Q o at the base of the conduit for the WVW bubble growth model, for two values of the viscosity response coefficient β.As discussed in Appendix B, the path ABCD describes the oscillatory cycle in (P , Q) space obtained in the limit where ε → ∞.(b) Steady state stability phase diagram in parameter space (Q o , ε) for β=5.5 for the WVW bubble growth model.S=stable node; SF=stable focus; U=unstable focus.The Hopf bifurcation line (continuous curve) resides within the region where the steady state P s (Q o ) exhibits two extrema (dashed vertical lines).

Fig. 2 .
Fig. 2. (a) Steady state overpressure P s as a function of the flow rate Q o at the base of the conduit for the SBG model with β=5.5.(b) Steady state stability phase diagram in parameter space (Q o , ε) for β=5.5 for the SBG bubble growth model.S=stable node; SF=stable focus; U=unstable focus.The Hopf bifurcation line (continuous curve) resides within the region where the steady state P s (Q o ) exhibits two extrema (dashed vertical lines).

Fig. 3 .
Fig. 3. (a) Steady state overpressure P s as a function of the flow rate Q o at the base of the conduit for the MBG model with β=5.5.(b) Steady state stability phase diagram in parameter space (Q o , ε) for β=5.5 for the MBG model.S=stable node; SF=stable focus; U=unstable focus.The Hopf bifurcation line (continuous curve) resides within the region where the steady state P s (Q o ) exhibits two extrema (dashed vertical lines).

Fig. 5 .
Fig. 5. Cycle period T as a function of ε for the MBG model with β=5.5 and for various values of Q o , as indicated on the curves.Parameters values as in Table1.

Fig. 6 .
Fig. 6.Limit cycle P (t), Q(t) for the MBG model with β=5.5 and for values of the time scale t and ε that minimize the relative error e (Table2).The origin of the time axis is arbitrary as the transients are not illustrated.

Table 2 .
Values of the parameter ε and the time scale t leading to a minimum in the error function e of Eq. (55) for various bubble growth models.Here, T obs =18 h and a obs =0.6.The dimensionless reservoir flow rate Q o is calculated from Eq. (22) with the knowledge of t and assuming Q o =10 m 3 /s.Also, T calc = tT calc .The product V η o is estimated from Eq. (22) with the knowledge of t and ε, and assuming Y =3×10 10 Pa and σ =0.2.P is the amplitude of the overpressure drop over one cycle.Q is the corresponding flow rate amplitude.χisdefined in Eq. (56).
a Wiley-Voight-Whitehead; b Single-Bubble Growth; c Multiple-Bubble Growth