Analysis of stochastic model for nonlinear volcanic dynamics

Motivated by important geophysical applications we consider a dynamic model of the magma-plug system previously derived by Iverson et al. (2006) under the influence of stochastic forcing. Due to strong nonlinearity of the friction force for a solid plug along its margins, the initial deterministic system exhibits impulsive oscillations. Two types of dynamic behavior of the system under the influence of the parametric stochastic forcing have been found: random trajectories are scattered on both sides of the deterministic cycle or grouped on its internal side only. It is shown that dispersions are highly inhomogeneous along cycles in the presence of noises. The effects of noise-induced shifts, pressure stabilization and localization of random trajectories have been revealed by increasing the noise intensity. The plug velocity, pressure and displacement are highly dependent of noise intensity as well. These new stochastic phenomena are related to the nonlinear peculiarities of the deterministic phase portrait. It is demonstrated that the repetitive stick–slip motions of the magma-plug system in the case of stochastic forcing can be connected with drumbeat earthquakes.


Introduction
It is well-known that the behavior of volcanic systems is enormously complex so that a lot of non-linear feedbacks lead to multiple states even during a single eruption (Tanaka et al., 2014). Without better modeling forecastings of these dynamic processes the highly important questions where, when and how volcanic eruptions occur will remain substantially empirical. Nowadays, an elaboration of the adequate mathematical models for volcanic dynamics is a challenging problem (Melnik and Sparks, 1999;Barmin et al., 2002;Nakanishi and Koyaguchi, 2008;Costa et al., 2012).
Many uncertainties in physical parameters of volcanic dynamics (Woo, 2000) lead to a conclusion that like the climate systems (see, among others, Saltzman, 2002;Alexandrov et al., 2014) volcanoes, representing stochastic and chaotic systems, need to Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | be described in terms of probabilities (Sparks, 2003;Bebbington and Marzocchi, 2011). Stochastic approaches and mathematical formalisms can be found in (Gardiner, 2009).
Some of silicic volcanoes analyzed in detail over the last few decades represent complex periodic systems (Denlinger and Hoblitt, 1999;Michaut et al., 2013). The dome-building eruption of Mount St Helen (MSH) during 2004 and 2005 has represented a nearequilibrium cyclic system with the solid plug uplift caused by magma ascent from below with a nearly steady rate of roughly 1-2 m 3 s −1 . This eruption was accompanied by drumbeat earthquakes recurred every 1-2 min with magnitudes < 2 and focal depth < 1 km (Iverson et al., 2006;Moore et al., 2008;Matoza and Chouet, 2010). A stick-slip mechanism explains a cyclic behavior of such earthquakes as a consequence of stick-slip motion of a plug pushed by compressible magma (Denlinger and Hoblitt, 1999). A new dynamic approach based on this mechanism (connecting the MSH behavior with a damped oscillator) was suggested by Iverson et al. (2006). We use this model to demonstrate unusual dynamic behavior of similar volcanic systems under the influence of parametric noises.
In order to explain interactions between solid-state extrusion and persistent drumbeat earthquakes at MSH, Iverson et al. (2006) developed a model based on recurrent stick-slip motions of the solid plug along its margins with the friction force F (Fig. 1). Let us briefly discuss the main principles of this dynamic model. The magma influx comes to the base of an eruptive conduit from below with a nearly steady rate Q. A solid dacite plug of solidified magma blocks the conduit from above so that its lower boundary is mobile due to the effects of pressure and basal accretion with mass rate ρB from below (ρ and B stand for the magma bulk density and the volumetric rate of magma crystallization). The total plug mass Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | m changes with time because the difference in mass rates ρB and ρ r E as m = m 0 + κt (here ρ r is the plug bulk density, E is the volumetric rate of surface erosion, m 0 is the initial plug mass, t is the process time, and κ = ρB − ρ r E is assumed constant). The horizontal cross-sectional area A, the magma compressibility α 1 and the conduit wall compliance α 2 are estimated by Iverson et al. (2006). The dynamic process of plug extrusion is controlled by the plug weight mg and the friction force F dependent of the plug velocity u (g is the acceleration due to gravity) whereas the conduit volume V is governed by the law of mass conservation. A three-parametric differential model connecting independent variables u, p and V (p is the pressure) was derived by Iverson et al. (2006). Below we use this model to demonstrate some new special aspects of non-linear dynamics of volcanic systems under the influence of stochastic noises.

The model and its deterministic behavior
The following system of reduced governing equations based on the laws of conservation of the solid plug linear momentum, solid plug mass and conduit fluid mass was derived and discussed in detail by Iverson et al. (2006). These equations can be written in the form of where (Iverson et al., 2006). Here, p 0 and ρ 0 represent the static equilibrium pressure and magma density. The key aspects of MSH friction force measured in experiments (Moore et al., 2008) can be described by a function (Iverson et al., 2006) where sgn(u) is the sign of u, F 0 is the friction force at static equilibrium, c ≪ 1 is a rateweakening parameter and u ref is a reference value of u (Iverson et al., 2006). Equation (4) includes the main physical aspects of the process: the friction force at u = 0 abruptly changes its sign due to the fact that the gravity force, which shifts the plug in downward direction, is opposite to the friction force. However, an abrupt behavior of expression (Eq. 4) is not good physical approximation of the friction force. Therefore, let us model this force by the close continuous function where In present paper, we focus on the autonomous case, when κ = 0. In Fig. 2, u-p-projections of phase trajectories of the system (Eqs. 1-3) for the fixed initial value V 0 = 6.32 × 10 5 m 3 are plotted by the thin black lines. These trajectories tend to the close curve (thick black line) of the cycle. Time series of this cycle are presented in Fig. 3. A vertical left part of this cycle in Fig. 2 corresponds to the slow movement, and the rest arc part of the cycle reflects fast movement. Slow dynamics switches to fast one at the corner point C in the case of movement in a clockwise direction along the thick black curve. The Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | stability of this cycle is highly non-uniform: the vertical part is extremely stable, but arc curve possesses a neutral stability.
Essential details of the phase portrait are shown in Fig. 2b by an enlarged fragment of Fig. 2a. As one can see, there exists a pseudo-separatrix (dashed red line) which divides two types of dynamics. If the initial state lies to the left of this red curve, then the trajectory quickly verges towards the vertical part of the cycle (arrows pointing to the left). If the initial state lies to the right of this pseudo-separatrix, then the phase trajectory goes away from the cycle (arrows pointing to the right), and only after a long excursion, the trajectory approaches to the vertical part of the cycle. Physically it means that small deviations in u at sufficiently large p may redeploy the dynamic system through its pseudo-separatrix.
This feature of the deterministic phase portrait playing an important role in understanding of stochastic phenomena will be discussed below.
Note, that this cycle is not a limit cycle in the classical mathematical sense. Indeed, for different values of V 0 , the non-linear system (Eqs. 1-3) exhibits different closed curves. The cycles and time series for various values of V 0 are compared in Fig. 4. Note that an increase of V 0 implies an increase of both amplitude and period of oscillations.

The role of stochastic forcing
In order to study possible deviations of the friction force from expression (Eq. 5) let us consider the parametric random disturbances. Such disturbances simulate the influence of different physical processes and phenomena leading to variations in the friction force behavior (e.g. the effects of frictional melting, temperature-dependent friction, and so on).
If the noise intensity is small enough, such bundle has a small dispersion and is localized near the deterministic cycle (green lines in Fig. 5a). As the noise intensity increases, along with the natural increase of dispersion, the following unexpected phenomenon is observed: the bundle's right side of stochastic trajectories is shifted inside the deterministic cycle (blue and red lines in Fig. 5a). Some details of the corresponding probabilistic distributions are presented in Figs. 5a,b. The probability density functions of u coordinates of intersection points of the random trajectories with the line p = 1.2935×10 7 Pa are plotted for three values of the F -noise intensity in Fig. 5b whereas the probability density functions of time intervals between successive intersections are shown in Fig. 5c. As one can see, with increasing noise, both the amplitude and period of stochastic oscillations decrease.
This stochastic phenomenon can be explained by the phase portrait peculiarities of initial deterministic system (see Fig. 2b) near the upper part of vertical fragment of the cycle. In the deterministic case, the phase trajectory slowly moves along the vertical part of the cycle up to the point C. At the point C, this trajectory abruptly changes the direction, and begins to move along the arc part quickly. Under the stochastic disturbances, random trajectories deviate from this vertical part of the deterministic cycle. As a result of this deviation, the random trajectory can cross the red pseudo-separatrix, and then it falls within the region of large arc-form excursions. In this case, the random trajectory turns right before the point C. The more noise, the earlier this turn. Such stochastic deformation of the random flow results in a decrease of u-and p-oscillation amplitudes and of the period.
Under the further increasing of noise intensity, random states of the system (Eqs. 2, 3 and 6) are localized and leave the interior of deterministic cycle. This noise-induced shift is demonstrated in Fig. 6. Here, an essential decrease of the dispersion of p coordinate is observed. In other words, p stabilizes near its certain value with increasing the noise intensity.
The dynamics of plug displacement is shown in Fig. 7. If the noise intensity is large enough so that the system leaves its cycle, the plug displacement increases with noise. If Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the system is within its cycle, the displacement is also within the corresponding deterministic stepwise curve (black line in Fig. 7).
In order to study an influence of possible changes in magma influx, let us consider a role of Q-noise: Q → Q(1 + δξ(t)), where ξ(t) is a Gaussian white noise, and δ is a Q-noise intensity. In this case, a non-linear dynamic system consists of Eq. (1) as well as of the following stochastic equations For weak noise, stochastic trajectories are localized near the deterministic cycle (Fig. 8a). As noise intensity increases, the dispersion of random trajectories increases as well (Figs. 8b,c). It can be seen that the dispersion is extremely non-uniform along the cycle. For the vertical part, a dispersion is small even for large noise, and the random trajectories do not differ from the deterministic cycle. Along the arc part, the dispersion of the random trajectories increases. The probability density functions of u coordinates of intersection points of the random trajectories with the line p = 1.2935 × 10 7 Pa are plotted for three values of the Q-noise intensity in Fig. 9a. Panel b of this figure shows the probability density functions of time intervals between successive intersections. As one can see, the dispersions grow and the mean values are practically unchangeable with increasing noise.
As one can see, the volcanic model under consideration demonstrates quite different qualitative and quantitative response on the random perturbations of different parameters. The system is extremely sensitive to F -noise so that even a weak F -noise implies a crucial deformation of the oscillatory behavior. An increase of F -noise leads to a decrease of the period and amplitude of oscillations. Note that the system is also sensitive to Q-noise so that a large noise intensity implies a dispersion increase of the arc-part of stochastic oscillations.

Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 4 Conclusions
The phase portrait of deterministic system contains a point of unstable equilibrium and a pseudo-separatrix, which subdivides the system into different dynamic areas (point C and dashed red line in Fig. 2). In addition, if point (u, p) of the phase plane lies below this pseudo-separatrix (Fig. 2b), the system quickly reaches its equilibrium state (thick vertical line in Fig. 2). If, however, the phase point (u, p) is above this pseudo-separatrix, the system tends to its equilibrium state in course of long time interval. An important point of the deterministic behavior is that a certain constant value of the plug velocity u establishes at different pressures p in the state of equilibrium (thick vertical line in Fig. 2b). This equilibrium state also persists at different conduit volumes V 0 (Fig. 4). A time of system transition to its equilibrium state (vertical line in Fig. 4a) therewith increases with increasing V 0 . In additiion, more broad conduits might have a rather big variation in u and p than more narrow ones during the deterministic process of volcanic plug evolution. By this is meant that a time required to attain the equilibrium state increases with increasing the conduit's cross-sectional area.
In order to analyze a role of variations of two main parameters of the plug motion (friction force F and magma influx Q), two types of noises have been introduced in the model equations: F -noise and Q-noise. It was shown that these noises lead to different evolutionary types of the dynamic system. Let us summarize the main aspects of those behavior. In the first place, random trajectories are scattered either on both sides of the deterministic cycle (Q-noise) or on its internal side (F -noise). Dispersions corresponding to random trajectories of both noises upon that grow with increasing the noise intensities. As this takes place, one can see that both dispersions are highly inhomogeneous along cycles (Figs. 5 and 8). Note that these dispersions are small enough in the vertical parts of corresponding cycles for F -and Q-noises.
An important point is that, an increase in dispersion occurs in the vicinity of pseudoseparatrix under the influence of F -noise. This is due to the fact that phase trajectories intersect the pseudo-separatrix and the phase points undergo transitions across it under Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the action of F -noise. Let us especially emphasize that F -noise phase trajectories leave the corresponding deterministic cycle and form some stochastic bundle shifted into the cycle's interior. As this takes place, the bundle's dispersion increases while the period and amplitude of oscillations decrease with increasing the F -noise intensity (Fig. 5). By this is meant that the presence of F -noise reduces possible variations in the plug velocity u and pressure p and decreases a time required to attain the equilibrium state (thick vertical line in Fig. 5a). It is significant that the effect of pressure stabilization near a certain value (dependent of the noise intensity) occurs with a rise in the F -noise intensity. The random trajectories therewith leave the corresponding cycle and localize in the vicinity of this value (Fig. 6). A dynamic behavior of the plug displacement is dependent of whether the dynamic system is within or beyond its phase cycle. In the former case, the plug displacement oscillates within the bounds of the corresponding deterministic stepwise curve (Fig. 7). In the latter case when the noise intensity is sufficiently large, the plug displacement increases drastically.
It is known that the eruption of Mount St. Helens was accompanied by rather regular repetitive long-period (or drumbeat) earthquakes over a long time. Moreover, such drumbeat events were more random from time to time. In addition, subevents in the form of randomly occurring series of smaller seismic events (produced by a separate random process) have been imposed upon these long-period events (Matoza and Chouet, 2010). The present study demonstrates that repetitive stick-slip motions of the plug representing stochastic oscillations can be connected with these drumbeat earthquakes. The calculated period between drumbeats (see Figs. 5 and 9) is in agreement with experimental data (30-300 s, Iverson et al., 2006). The physical reason is that such earthquakes observed at shallow depths (< 1 km at MSH) can be caused by the stick-slip motions of the magma-plug system under the influence of noises where the driving force acting on a compliant crustal body is large enough (the force drop responsible for this kind of seismicity can be estimated from our calculations as ∆F = ∆pA ∼ 6 × 10 7 kg m s −2 , where ∆p ∼ 2 × 10 3 Pa).