The dynamics of dual-magma-chamber system during volcanic eruptions inferred from physical modeling

A magma plumbing system with dual magma chambers beneath active volcanoes is commonly observed through petrological and geophysical measurements. This paper developed a physical model for the dynamics of a dual-magma-chamber system during volcanic eruptions. The model consists of the plumbing system where two elastically deformable magma chambers are connected in series with non-deformable conduits. Based on this model, we obtained an analytical solution that describes temporal changes in pressures at the two chambers accompanied by the eruption. The analytical solution showed that the feature of the chamber pressure changes is mainly controlled by two non-dimensional numbers C′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C'$$\end{document} and Ω′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega '$$\end{document}. Here, C′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C'$$\end{document} is the ratio of the parameter controlling the magnitude of pressure change in the shallower chamber to that in the deeper chamber, and Ω′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega '$$\end{document} is the ratio of conduit’s conductivity (inverse of resistivity to magma flow) between the shallower chamber and the surface to that between the chambers. For smaller C′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C'$$\end{document} and Ω′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega '$$\end{document}, the shallower chamber’s pressure is kept constant during the decrease in the pressure at the deeper chamber in the initial phase of the eruption. This corresponds to a deformation pattern commonly observed in some eruptions, in which the deflation of the deeper chamber was predominant. The estimation of C′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C'$$\end{document} and Ω′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega '$$\end{document} based on the parameters related to magma properties and geometries of the chambers and the conduits revealed that the smaller C′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C'$$\end{document} and Ω′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega '$$\end{document} conditions are satisfied under realistic magmatic and geological parameters. This indicates that the magma dynamics in the dual-chamber system generally cause the dominance of the deeper chamber’s deflation.


Introduction
A magma plumbing system beneath an active volcano shows a wide variation in size, geometry, and magma properties inside it. This variation causes complex features of geodetic and petrological observations during volcanic eruptions (e.g., Dvorak and Dzurisin 1997;Cashman et al. 2017;Edmonds et al. 2019). To retrieve information about the eruption dynamics from the observations, it is necessary to understand the relationship between the magma plumbing system and observed signals. The simplest magma plumbing system is composed of a single and elastically deformable magma chamber connected to the surface by a volcanic conduit, in which inflation or deflation of the magma chamber is caused by a balance between a magma influx from depth to the chamber and a magma outflux from the chamber to the conduit. In this case, the inflation/ deflation of the magma chamber, which can be detected by the geodetic observations of the surface deformation, is straightforwardly linked with the magma activity in the plumbing system (e.g., Mogi 1958;Dzurisin 2003;Segall 2013). When additional effects such as magma lateral migration from the magma chamber and viscoelastic deformation of the chamber are considered, the relationship between the chamber deformation and the magma activity becomes more complex (e.g., Ueda et al. 2005;Segall 2016; Anderson et al. 2019).
One of the complex magma plumbing systems observed in the volcanoes with silicic eruptions is the system composed of dual magma chambers. In some petrological studies, the dual-chamber system has been proposed to explain the magma mixing process: two types of magma with lower and higher degrees of fractionation (i.e., lower and higher SiO 2 ) are stored in the deeper and shallower chambers, respectively, and the magma injection from the deeper to shallower chambers leads to the magma mixing (Sparks et al. 1977). The magma chamber's depth is deduced from the magma storage pressure that is estimated from thermodynamic constraints for petrological data. In geodetic analyses for volcanic activities such as based on GPS, tiltmeter, strainmeter, and InSAR, multiple deformation sources are sometimes detected to explain the surface deformation patterns (e.g., Kohno et al. 2008;Foroozan et al. 2011;Reverso et al. 2014). In this case, the chamber depth is directly estimated from the depth of the deformation source.
In the dual-magma-chamber system's observations, there is a concern that the petrological imaging of the chamber system does not always agree with the geodetic imaging. In the observations during some well-monitored eruptions, there is a common feature that the dual chambers at the shallower and the deeper depths were inferred from the petrological observations, whereas the geodetic observations mainly detected the deformation of the deeper chamber. At Mount Usu, Hokkaido, Japan, two magma chambers with storage pressures (depths) of 200-250 MPa (8-10 km) and 100-150 MPa (4-6 km) were inferred from petrological experiments (Tomiya et al. 2010). On the other hand, geodetic observations during the first few days of the 2000 eruption report a deflation of a deeper magma chamber at a depth of 10 km (Murakami et al. 2001). In the 2011 eruption at Kirishima Shinmoe-dake volcano, Kyushu, Japan, petrological observations revealed that basaltic andesite magma and silicic andesite magma were stored in deeper and shallower magma chambers with pressures (depths) of 250 MPa (10 km) and 125 MPa (5 km), respectively (Suzuki et al. 2013). The erupted magma originates from the mixing of these magmas. Precise geodetic observations based on borehole-type tiltmeter and GPS measurements detected a deflation of a spherical source at a depth of 9.8 km that was accompanied with three sub-Plinian eruptions with the duration of 2-3 h in each event and lava extrusion with two days duration, whereas a shallower source was not detected (Ueda et al. 2013). At Sakurajima volcano, Kyushu, Japan, petrological observations based on melt inclusion data revealed that magma storage depth immediately before the 1914 Plinian eruption is 0.9-3.2 km, and a few data show a storage depth of 4.1 km (Araya et al. 2019). A geodetic deformation source accompanied by recent Vulcanian explosions is estimated to be located at a depth of 4 km beneath the crater (Iguchi et al. 2013), implying a shallower magma chamber. In contrast, a precise leveling survey revealed that a sudden deflation source during the 1914 eruption was deeper, located at a depth of 10 km beneath Aira caldera (Mogi 1958). These observations at three volcanoes show a common feature: a deeper chamber deflation was predominant during the eruption in the dual-magma-chamber system.
The dual-chamber system dynamics have been investigated in some previous studies, mainly focusing on the magma plumbing system's long-term evolution. At Soufrière Hills volcano on Montserrat, West Indies, a magma transfer between the dual chambers connected by a vertical conduit has been imaged by inverting magma efflux data and surface deformation data by GPS during eruptive and inter-eruptive periods for about 12 years (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) (Elsworth et al. 2008;Foroozan et al. 2011). These results also showed that the eruptive episodes deplete the deeper chamber only (Elsworth et al. 2008), whereas whether the shallower chamber's deflation is apparent or not depends on model assumptions of the plumbing system and the GPS data quality (Foroozan et al. 2011). Melnik and Costa (2014) developed a numerical model for the plumbing system composed of two magma chambers connected by a dyke and a shallower conduit with dyke-to-cylindrical shape by referring to Soufrière Hills volcano. The model showed that connectivity between the chambers significantly influences the dynamics, such as the timing of chamber pressure change and the period of cyclic behavior of discharge rate corresponding to 2-3 years eruption cyclicity. At Grímsvötn Volcano, Iceland, posteruptive displacements measured by dense GPS network after the 1998, 2004, 2011 eruptions were explained by an analytical model that has two connected magma chambers (Reverso et al. 2014). More recently, Bato et al. (2018) applied a sequential data assimilation technique (Ensemble Kalman Filter) to the GPS data after the 2011 eruption based on the two-chamber model and revealed the temporal evolution of the basal magma inflow and a possible interplay between the activities of Grímsvötn volcano and its adjacent volcano (Bárðarbunga).
Recently, the short-term geodetic deformations accompanied with the eruptions were successfully measured by high-resolution observations such as based on tiltmeter and strainmeter (e.g., Bonaccorso et al. 2013;Iguchi et al. 2013;Ueda et al. 2013). As shown in the previous studies for the long-term geodetic data, an appropriate physical model is also necessary to observe the short-term phenomena. In this study, we try to develop a physical model for the dual-magma-chamber system dynamics during the eruption. Using this model, we obtain an analytical solution for temporal changes in pressures at the two chambers accompanied by the eruption, which describes the relationship between a pattern of the chamber pressure changes and magmatic and geological conditions. We discuss whether a deformation pattern in which the pressure decrease in a deeper magma chamber is predominant, commonly observed in some eruptions, can occur under realistic conditions.

A model for dual-magma-chamber system
We consider a magma plumbing system in which two elastically deformable magma chambers are connected in series with non-deformable conduits (Fig. 1). A pressure change in each magma chamber is controlled by the balance between magma influx from the deeper conduit and outflux to the shallower conduit. In this case, the temporal changes in the chamber pressures are formulated as and where (1) dp p is the chamber pressure, t is time, Q in and Q out are the magma influx and outflux, respectively (in kg s −1 ), ρ is the magma density, V is the chamber volume, κ m , is the magma compressibility in the chamber, and κ ch is the compressibility of the magma chamber. The number subscripts "1" and "2" correspond to the deeper and shallower magma chambers, respectively ( Fig. 1). Here, C 1 and C 2 are the parameters controlling a variability of the chamber pressure: a larger C 1 or C 2 leads to a larger magnitude of the chamber pressure change for given magma influx and outflux. The magma flow in the conduit is modeled as a Poiseuille flow. To set the conduit geometry as general as possible, we consider a dyke-like conduit with an elliptical cross section (Costa et al. 2007, Fig. 1). Under the assumption that a temporal change in the magma density in the conduit is negligible, the mass conservation in the conduit between the two chambers is expressed as where w is the ascent velocity, a and b are the semi-major and semi-minor axes of the ellipse, respectively, and the subscript "12" represents the conduit between the two chambers ( Fig. 1). In the momentum conservation, we assume a linear pressure gradient and that the inertia term is negligible. In this case, we obtain: where L is the conduit length, g is the gravity acceleration, and η is the magma viscosity. Similarly, the mass and momentum conservations in the conduit between the shallower chamber and the surface are approximated as and where p s is the pressure at the surface, and " 2s " represents the conduit between the shallower chamber and the surface (Fig. 1). From Eqs. (4)-(7), we obtain and (3) where The parameters 12 and 2s correspond to "conduit conductivity" (Slezin 2003): a larger 12 or 2s leads to a higher magma flux for a given pressure gradient in the conduit due to lower viscous resistance. Here, the chamber pressure is expressed using the overpressure ( p ) with respect to lithostatic pressure as where ρ c is the crust density ( ∼ 2500 kg m −3 ). From Eqs. (1), (2), (8), (9), and (11) and under the assumption that p s ≪ p 2 , we obtain and Eqs. (12) and (13)  and � 2s = ρ 2s πa 3 2s b 3 2s 4η 2s L 2s (a 2 2s + b 2 2s ) . (11) It is noted that when ρ 12 and ρ 2s in Eq. (17) are approximated as (1 − φ 12 )ρ lc and (1 − φ 2s )ρ lc , respectively, where φ is the magma porosity in the conduit and ρ lc is the liquid-crystals density which is close to ρ c , ρ ′ corresponds to the ratio of the porosity in the shallower conduit to that in the deeper conduit, φ 2s /φ 12 . Here, the non-dimensional numbers C ′ and ′ are the ratio of C 2 to C 1 and that of 2s to 12 , respectively, expressed as and where κ ′ is the ratio of the sum of the magma and the chamber compressibilities in the shallower chamber to that in the deeper chamber: and r is the ratio of the semi-major to semi-minor axes of the ellipse in the conduit cross section: From the simultaneous differential equations (15) and (16) and under the assumption that ρ ′ , L ′ , C ′ , ′ , and Q ′ in1 are constant during the pressure changes, we obtain an analytical solution for p ′ 1 and p ′ 2 as and where (21) When we use the boundary conditions that p ′ 1 = p ′ 1,0 and p ′ 2 = p ′ 2,0 for t ′ = 0 , A + and Aare expressed as

Results
We investigate temporal changes in the chamber pressures during an eruption using the analytical solution (Eqs. (23) and (24)). For simplicity, the magma influx to the deeper chamber Q in1 is set to be 0 under the assumption that it is substantially smaller than the magma flux in the conduits during the eruption, and p ′ 1,0 and p ′ 2,0 are set to be 0. We first show the typical results for temporal changes in p ′ 1 and p ′ 2 for ρ ′ L ′ = 3 (Fig. 2). This value of ρ ′ L ′ is obtained under realistic conditions in which, for example, ρ 12 /ρ c = 0.9 , ρ 2s /ρ c = 0.7 , and L 2s /L 12 = 1 . These conditions correspond to the case that the average magma porosities in the deeper and the shallower conduits ( φ 12 and φ 2s ) are 0.1 and 0.3, respectively, and the deeper magma chamber is two times deeper than the shallower chamber. The temporal changes in p ′ 1 and p ′ 2 strongly depend on the non-dimensional parameters C ′ and ′ . When C ′ = ′ = 1 , both p ′ 1 and p ′ 2 decrease simultaneously (Fig. 2e). For larger C ′ and ′ , the decrease in p ′ 2 is followed by the decrease in p ′ 1 (Fig. 2c). On the other hand, for smaller C ′ and ′ , the decrease in p ′ 1 is followed by the decrease in p ′ 2 (Fig. 2g). In this case, p ′ 2 is kept constant during the decrease in p ′ 1 in the initial phase. This feature reflects the fact that a larger C 1 (i.e., higher variability of the pressure in the deeper chamber) or a larger 12 (i.e., higher conductivity of the deeper conduit) leads to an efficient pressure decrease in the deeper chamber. For larger C ′ and smaller ′ , p ′ 2 increases during the decrease in p ′ 1 (Fig. 2i). The deformation pattern observed during some eruptions, in which the deeper magma chamber's deflation was predominant, corresponds to the case of smaller C ′ and ′ (Fig. 2g). The dependence of the features of the temporal changes in the chamber pressures on the parameters can be evaluated by checking the magnitude of p ′ 2 when p ′ 1 reaches a critical value (Fig. 3). To investigate the features of the pressure changes during the initial phase of the eruption, we set the critical value of p ′ 1 as −0.5 or −1 . Figure 3 shows the contour map for p ′ 2 at the critical value of p ′ 1 in the parameter space of C ′ and ′ for ρ ′ L ′ ranging from 1 to 10. The range of ρ ′ L ′ was set under the assumptions that φ 12 = 0.1 , φ 2s = 0.2-0.5, and L ′ = 0.5-2. The region of negative p ′ 2 smaller than the critical value of p ′ 1 corresponds to the case in which the decrease in p ′ 2 is followed by the decrease in p ′ 1 (i.e., Fig. 2c). On the other hand, the region of p ′ 2 ∼ 0 corresponds to the case in which the decrease in p ′ 1 is followed by the decrease in p ′ 2 (i.e., Fig. 2g). The region of positive p ′ 2 at p ′ 1 = −0.5 (Fig. 3a, c, and e) reflects the case in which p ′ 2 increases during the decrease in p ′ 1 (i.e., Fig. 2i). Figure 3 indicates that the typical patterns of the temporal changes in the overpressures shown in Fig. 2 exist for a wide range of ρ ′ L ′ . From the quantitative point of view, the boundary of the region depends on ρ ′ L ′ : As ρ ′ L ′ increases, the region of negative p ′ 2 expands to the smaller ′ range, whereas the region of positive p ′ 2 contracts. In the region of smaller C ′ and ′ , we can specify a condition for p ′ 2 to be close to 0 during the decrease in p ′ 1 . As an example, we now specify the condition of the ranges of C ′ and ′ for p ′ 2 to be in the range of −0.1 to 0.1 both at p ′ 1 = −0.5 and −1 (referred to as "smallp ′ 2 condition"). When ρ ′ L ′ = 1 , the small-p ′ 2 condition is satisfied for C ′ < 0.200 and � ′ < 0.493 (Fig. 3a, b). As ρ ′ L ′ increases, the range of C ′ is unchanged, whereas that of ′ decreases: � ′ < 0.048 at ρ ′ L ′ = 10 (Fig. 3e, f ). Therefore, under a wide range of ρ ′ L ′ , the small-p ′ 2 condition is satisfied when C ′ < 0.200 and � ′ < 0.048.

Discussion
We have specified the ranges of C ′ and ′ for the smallp ′ 2 condition in the previous section. This condition is interpreted as the case in which the deflation of the deeper magma chamber is predominant. Here, we discuss whether the small-p ′ 2 condition is satisfied under realistic magma properties and geological conditions. Among the parameters included in C ′ and ′ (Eqs. (19) and (20)), V 2 /V 1 , L ′ , b 2s /b 12 , r 2s , and r 12 are related with the geometry of the magma plumbing system such as the chamber volume and the conduit length and depth, and ρ 2 /ρ 1 , ρ 2s /ρ 12 , and η 2s /η 12 are determined by the Fig. 4 Magma compressibilities in the deeper ( κ m1 ; solid curves) and the shallower ( κ m2 ; dashed curves) chambers (a), κ ′ for S c1 = S c2 = 1 (b), and κ ′ with varying S c1 and S c2 (c), as a function of H 2 O. The value of CO 2 is varied from 500 to 2000 ppm in a and b and set as 500 ppm in c. See text for details magma properties such as the density and the viscosity. The parameter κ ′ , which is defined as the ratio of the sum of the magma compressibility κ m and the chamber compressibility κ ch (Eq. (21)), is a complex function of the magmatic and geological parameters: κ m is expressed by a function of the magma density ρ as (∂ρ/∂p)/ρ , and κ ch is controlled by the chamber shape and the rigidity of the surrounding rocks.
To obtain the relationship between κ ′ and the magmatic and geological parameters, we first calculate the magma compressibility by formulating the magma density as the density of the gas-liquid-crystals mixture, in which the gas density follows the equation of state for ideal gas and the liquid-crystals density is constant ( = 2500 kg m −3 ). The equilibrium volatile exsolution is assumed to follow the H 2 O-CO 2 solubility model by Newman and Lowenstern (2002). We also assume that the magma temperature is 1000 o C , the melt-crystals compressibility is 10 −10 Pa −1 , and the pressure and the crystal content in the deeper and the shallower chambers are 250 MPa and 0 vol.%, and 125 MPa and 30 vol.%, respectively. Under these formulations and assumptions, we can calculate the magma compressibilities in the deeper and the shallower chambers ( κ m1 and κ m2 ) as a function of initial H 2 O and CO 2 contents (Fig. 4a). Furthermore, the compressibility of the magma chamber is expressed by (3S c )/(4µ) where µ is the rigidity of the surrounding rocks set to be 10 GPa, and S c is a non-dimensional parameter used for evaluating the effect of the shape change in the magma chamber on the compressibility. The value of S c is equal to 1 when the magma chamber is purely spherical, and it increases as the chamber changes to a penny-shaped sill. When the ratio of the longer to shorter axes of the sill is 10, S c reaches up to about 10 (Anderson and Segall 2011). We defined S c for the deeper and the shallower chambers as S c1 and S c2 , respectively.
In the case that both the deeper and the shallower chambers are spherical (i.e., S c1 = S c2 = 1 ), the variation of κ ′ ranges from 1 to about 20 (Fig. 4b), which originates from the effect of volatile exsolution on the magma compressibility (Fig. 4a). For higher H 2 O, the volatiles are exsolved and the gas phase is generated, leading a drastic increase in the magma compressibility. On the other hand, for lower H 2 O, the volatiles are not exsolved and the magma compressibility is equal to the melt-crystals compressibility. For given volatile contents, κ m2 is larger than κ m1 because lower pressure and higher crystallinity lead to more efficient volatile exsolution in the shallower chamber. Especially, κ ′ drastically increases under the conditions of larger κ m2 due to volatile exsolution and constant κ m1 without volatile exsolution (i.e., κ m1 is equal to the melt-crystals compressibility), as shown in the case of H 2 O = 4 wt % and CO 2 = 500 ppm (Fig. 4b). Figure 4c shows the effects of the change in the combination of the chamber shape on κ ′ , in which S c1 and/or S c2 are set as 10 to take into account the penny-shaped chamber. Compared to the case of S c1 = S c2 = 1 , κ ′ is larger for S c1 = 1 and S c2 = 10 , whereas it is smaller for S c1 = 10 and S c2 = 1 . When S c1 = S c2 = 10 , the increase in the magnitude of κ ch in κ ′ weakens the effect of the contrast between κ m1 and κ m2 , leading to a smaller κ ′ than the case of S c1 = S c2 = 1 . Nevertheless, κ ′ is larger than 1 except for the case of extremely small H 2 O. When S c1 = 1 and S c2 = 10 , which corresponds to the combination of the Fig. 5 a C ′ as a function of V 2 /V 1 for ρ 2 /ρ 1 = 0.8 with varying κ ′ . Color map and labeled contours for κ ′ are shown. b ′ as a function of η 2s /η 12 for r 2s = 1 , L ′ = 1 , and ρ 2s /ρ 12 = 0.8 with varying r 12 and b 2s /b 12 under the constraint that the cross-sectional area of the conduit is constant. The ranges of C ′ and ′ in which −0.1 < �p ′ 2 < 0.1 both at p ′ 1 = −0.5 and −1 for ρ ′ L ′ = 10 (Fig. 3c) are shown by red dashed lines and arrows Kozono Earth, Planets and Space (2021) 73:103 deeper spherical chamber and the shallower pennyshaped chamber, κ ′ reaches up to 20. Now we investigate the assemblage of the magmatic and geological parameters which satisfies the small-p ′ 2 condition (i.e., C ′ < 0.2 and � ′ < 0.048 ). Figure 5a shows the dependence of C ′ on the ratio of the chamber volume ( V 2 /V 1 ) and κ ′ . If the volumes of the two chambers are comparable (i.e., V 2 /V 1 ∼ 1 ), a larger κ ′ due to the effect of the compressibility contrast between the chambers plays a key role in satisfying the condition of C ′ < 0.2 . Under an extreme condition that κ ′ = 1 , in which the two chambers with the same shape are filled with incompressible magma, a larger V 2 /V 1 up to about 10 is needed to satisfy the condition of C ′ < 0.2 . On the other hand, when κ ′ reaches up to 20 because of larger magma compressibility due to efficient volatile exsolution and larger chamber compressibility due to penny shape in the shallower chamber (Fig. 3b, c), the condition of C ′ < 0.2 is satisfied even in the case of V 2 /V 1 < 1 . Figure 5b shows the dependence of ′ on the ratio of the magma viscosity in the conduit ( η 2s /η 12 ) and conduit geometry. In the case that the conduit is entirely cylindrical and keeping the same size ( r 12 = r 2s = 1 and b 2s /b 12 = 1 ), the condition of � ′ < 0.048 is satisfied when η 2s /η 12 is larger than about 20. Here, we consider a combination of the deeper dyke-like conduit and the shallower cylindrical conduit ( r 12 > 1 and r 2s = 1 ) under the constraint that the cross-sectional area of the conduit is constant (i.e., a 12 b 12 = a 2s b 2s ). In this case, a larger b 2s /b 12 is set for a larger r 12 , which corresponds to a thinner dyke with a larger aspect ratio. Because the thinning of the dyke leads to the decrease in the conductivity of the deeper conduit and the increase in ′ , the minimum value of η 2s /η 12 satisfying the condition of � ′ < 0.048 is larger than the case of the entirely cylindrical conduit, and it reaches 10 3 for r 12 = 100 and b 2s /b 12 = 10 . Generally, the magma viscosity increases during magma ascent because of crystallization and the evolution of melt composition such as the decrease in H 2 O dissolved in the melt. Furthermore, the basaltic magma injection from the deeper to the shallower chambers has been inferred from many petrological studies, implying a lower magma viscosity in the deeper conduit. Therefore, η 2s /η 12 is estimated to be much greater than 1, and the condition of � ′ < 0.048 is easily satisfied. These results indicate that the deformation pattern characterized by the dominance of the deeper chamber's deflation can occur under realistic magmatic and geological conditions. Finally, we evaluate the timescale and the magnitude of the volume change of the magma chamber analyzed in this study. The non-dimensional timescale for p ′ 1 to reach −1 is of the order of 1 (i.e., t ′ ∼ 1 ; Fig. 2). Therefore, the dimensional timescale is approximated as t = 1/(C 1 � 12 ) (Eq. (14)), which is a function of magmatic and geological parameters included in C 1 and 12 ( a 12 , b 12 , L 12 , V 1 , κ m1 , κ ch1 , η 12 , ρ 1 , and ρ 12 ; see Eqs. (3) and (10)). Under the assumptions that a 12 ∼ b 12 ∼ 10 m, L 12 ∼ 5000 m, V 1 ∼ 10 8 -10 10 m 3 , κ m1 ∼ κ ch1 ∼ 10 −10 Pa −1 , η 12 ∼ 10 3 -10 5 Pa s −1 , and ρ 1 ∼ ρ 12 , t is estimated as about 10-10 5 s. From the definition of the chamber compressibility, the volume change of the deeper magma chamber ( δV 1 ) is related with the pressure change ( δp 1 ) as δV 1 = κ ch1 V 1 δp 1 . Because the dimensional magnitude of the pressure change in the deeper chamber ( |δp 1 | ) at p ′ 1 = −1 is equal to |p| = (ρ c − ρ 12 )gL 12 (Eq. (14)), the magnitude of the volume change ( |δV 1 | ) is expressed as κ ch1 V 1 (ρ c − ρ 12 )gL 12 . By setting the values of κ ch1 , V 1 , and L 12 same as those for the estimation of t and assuming that ρ 12 ∼ 2000-2250 kg m −3 (i.e., the porosity is 0.1-0.2), |δV 1 | is estimated as about 10 5 -10 7 m 3 . The estimated values of t and |δV 1 | cover the duration of the crustal deformation and the magnitude of the change in the chamber volume observed during the 2000 Usu and the 2011 Shinmoe-dake eruptions (Murakami et al. 2001;Ueda et al. 2013). Although we recognize that both t and |δV 1 | strongly depend on unconstrained parameters such as the chamber volume, a rough agreement between the analytical results and the observations implies that the analytical solution has the potential to be combined with the observation data during the eruptions.

Conclusions
We developed the physical model of the co-eruptive magma dynamics in the plumbing system, where the two elastically deformable magma chambers are connected in series with the conduits. Based on the physical model, we obtained the analytical solution for the temporal changes in the two chambers' pressures accompanied by the eruption. The solution showed that the feature of the pressure changes is mainly controlled by the non-dimensional numbers C ′ and ′ : C ′ is the ratio of the parameter controlling the magnitude of the pressure change in the shallower chamber to that in the deeper chamber, and ′ is the ratio of the conduit's conductivity between the shallower chamber and the surface to that between the chambers. We found that for smaller C ′ and ′ , the shallower chamber's pressure is kept constant during the decrease in the pressure at the deeper chamber in the initial phase of the eruption. This corresponds to the deformation pattern commonly observed in some eruptions, in which the deflation of the deeper magma chamber was predominant. The estimation of C ′ and ′ based on the parameters related to the magma properties and the geometries of the chambers and the conduits revealed that the smaller C ′ and ′ conditions are satisfied under realistic magmatic and geological