L Forestier-Coste

S Mancini

A Burgisser

F James

Numerical resolution of a mono-disperse model of bubble growth in magmas

Keywords: coupled ode system and diffusion equation, mass preserving scheme, bubble growth AMS Class.: 86A04, 8608, 65L20, 65M06

Growth of gas bubbles in magmas may be modeled by a system of differential equations that account for the evolution of bubble radius and internal pressure and that are coupled with an advection-diffusion equation defining the gas flux going from magma to bubble. This system of equations is characterized by two relaxation parameters linked to the viscosity of the magma and to the diffusivity of the dissolved gas, respectively. Here, we propose a numerical scheme preserving, by construction, the total mass of water of the system. We also study the asymptotic behavior of the system of equations by letting the relaxation parameters vary from 0 to ∞, and show the numerical convergence of the solutions obtained by means of the general numerical scheme to the simplified asymptotic limits. Finally, we validate and compare our numerical results with those obtained in experiments.

Introduction

All volcanic eruptions involve a decompression of the magma during its ascent from the Earth's crust to the surface. This decompression causes the volatiles dissolved into the magma to come out of solution as gas bubbles. The way these bubbles are growing, whether they coalesce with one another or travel faster than or with the magma, are all conditioning the way the 1 volcanic eruption will unfold. Bubbles that remain trapped with the magma they originally grew from will accumulate gas pressure until failure of the magma releases it suddenly to produce an explosive eruption. Such scenario is most likely when the magma is highly viscous and prevents bubble motion. This situation is propitious to modeling because bubbles can be considered as immobile with respect to the magma and the resulting spherical geometry allows one to reduce bubble growth to a system of differential equations describing the evolution of pressure and gas mass in a bubble coupled with an advection-diffusion equation describing the drainage of the dissolved gas towards the bubble. A further assumption is that bubbles are exclusively made of water vapor, which can be justified by the fact that water is, by far, the most abundant volatile species in such viscous magmas.

Since the seminal work done in [START_REF] Scriven | On the dynamics of phase growth[END_REF] several numerical schemes that solve such system of differential equations have been proposed in the context of visco-elastic fluids (see [START_REF] Arefmanesh | Diffusion-induced growth of a gas bubble in a viscoelastic fluid[END_REF][START_REF] Vrentas | Slow bubble growth and dissolution in a viscoelastic fluid[END_REF][START_REF] Welch | Direct simulation of vapor bubble growth[END_REF][START_REF] Beechem | Bubble growth mechanism in carbon foams[END_REF][START_REF] Rosebrock | Study of the growth and motion of graphitic foam bubbles[END_REF]). Application to gas bubble in magmas is slightly more recent (see [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF][START_REF] Proussevitch | Dynamics and energetics of bubble growth in magmas: Analytical formulation and numerical modeling[END_REF][START_REF] Toramaru | Numerical study of nucleation and growth of bubbles in viscous magmas[END_REF][START_REF] Blower | Coupling of viscous and diffusive controls on bubble growth during explosive volcanic eruptions[END_REF][START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF][START_REF] Chouet | Dynamics of diffusive bubble growth and pressure recovery in a bubbly rhyolitic melt embedded in an elastic solid[END_REF]). All these schemes have in common a discretization of the advection-diffusion equation that is not conservative by construction with respect to the diffused species. In fact, they involve user-defined discretization parameters that have to be empirically adjusted to ensure sufficient convergence and/or accuracy of the scheme.

Developing alternate, robust schemes would allow including the dynamics of bubble growth into more sophisticated model that take into account, for instance, that bubble have different sizes, or that, if magma viscosity is low enough, bubble may rise with respect to the magma. In this paper we present a new numerical scheme, in which the flux in the advection-diffusion equation is computed in order to conserve the total water mass in the bubble-magma system at a discrete level, and this despite the mesh discretization we apply. Moreover, under the assumption of constant in time diffusion coefficient, we give some explicit solutions of the proposed model when the viscosity or the diffusion are very large (infinity) or very small (zero), we shall call these asymptotics limit regimes; we also numerically verify the convergence of the proposed scheme towards these limit regimes.

The present work is developed as follows. In section 2, we recall the differential equations describing the respective evolution of bubble radius and mass, together with the advection-diffusion equation describing the behavior of the water concentration in the magma. Following [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF][START_REF] Proussevitch | Dynamics and energetics of bubble growth in magmas: Analytical formulation and numerical modeling[END_REF] and [START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF] we write the problem in dimensionless form, introducing two relaxation parameters Θ V and Θ D . Section 3, is devoted to the numerical approximation of the model. The main novelty is the discretization of the advection-diffusion equation, see section 3.2, in which we explain how to compute the mesh and flux at each iteration in such a way that the total mass is conserved. In section 4 we deal with the asymptotics of the dimensionless problem, when the ratio between the relaxation parameters varies from 0 to ∞. Three main regimes are underlined: viscous, diffusive, and equilibrium. For each limit, we also propose a way to discretize it. Numerical results, convergence of the solution towards the simplified asymptotic limits, comparisons with experiments and with the code of reference [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF], are discussed in section 5. Finally, in section 6, we summarize our study and suggest possible extensions of the modeling of bubbles growth in magma.

The model

We are interested in the modeling of bubble growth in a highly viscous, crystal-free magma. This has two main consequences on the model. The first one is that we assume that bubbles do not interact with each other, in particular there are no coalescence effects. This is strongly limitative for the simulation of a magmatic conduit, but the presence or absence of coalescence can be controlled in laboratory experiments, see section 5.2. The other point is that, due to the high viscosity of the magma, bubbles travel along with the same velocity as the melt. In other words, they can be considered as immobile with respect to the melt.

At this stage, we can consider that a bubble can be described with two parameters, its volume V and its gaseous mass M . In this section, we denote with a hat the dimensional variables. Taking into account that the bubble is made only of water in gaseous form, we can write the perfect gases law inside the bubble in order to relate the gas pressure P to the gas density ρ:

GT ρ = P M w , (1)
with M w the molar mass for water, G the perfect gas constant and T the gas temperature. Next, following [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF][START_REF] Proussevitch | Dynamics and energetics of bubble growth in magmas: Analytical formulation and numerical modeling[END_REF] and [START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF] we assume that the bubble is spherical, with radius R, so that V = 4π R3 /3, and we set for future convenience M = ρ R3 , so that the bubble mass is 4π M /3. Thus we can choose the radius R = R(t) and the variable M = M (t), proportional to the mass, to describe the evolution of the bubble, and seek for a system of differential equations for these variables. Notice that in [START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF] an equation on the pressure P (t) is given, we choose here to track the bubble mass because it leads to a better handling of mass conservation at the numerical level. Such a model gives a description of the growth of a single bubble, or for a population of identical, non-interacting bubbles: this is the so-called mono-disperse case.

Basic equations

Two main physical processes drive bubble growth, both originating from magma decompression caused by magma ascent towards the surface. On the one hand, the gas trapped into the bubble is expanding; on the other hand, the water dissolved in the magma is diffusing and eventually is vaporized in the bubble, so that the water concentration profile in the melt has to be considered as well.

The equations describing the time evolution for the bubble radius and pressure and for the water concentration have been described in the literature several times. Therefore, we shall not reproduce here this derivation and, for example, we refer the reader to [START_REF] Arefmanesh | Diffusion-induced growth of a gas bubble in a viscoelastic fluid[END_REF][START_REF] Heureux | A new model of volatile bubble growth in a magmatic system: Isobaric case[END_REF][START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF][START_REF] Proussevitch | Dynamics and energetics of bubble growth in magmas: Analytical formulation and numerical modeling[END_REF][START_REF] Scriven | On the dynamics of phase growth[END_REF]. We recall briefly the origin of each equation, and their coupling.

The growth of the bubble in the magma involves viscous effects, surface tension effects and the ambient pressure in the magma Pa ( t), which is a given function of time. We define the decompression rate ∆P for a given ambient pressure function Pa ( t) as

∆P = P i t * ,
where P i is the initial ambient pressure and the time t * is such that Pa (t * ) = 0.

From the momentum conservation of Navier-Stokes equation, neglecting the inertial terms and considering the incompressibility of the melt, one obtains

P ( t) -Pa ( t) = 2σ R( t) + 4 Ṙ( t) R( t) ηeff , ( 2
)
where σ is the surface tension, assumed constant in this paper, and ηeff is the effective viscosity. In the following the viscosity may be constant or vary in space as:

ηeff (r) = η( R) 1 - η( Ŝ) η( R) α + R3 η( R) Ŝ R dη( r ) dr 1 r3 dr , (3)
with r the radial distance from the bubble radius outwards, η given by [START_REF] Hess | Viscosities of hydrous leucogranitic melts: a non-Arrhenian model[END_REF], and Ŝ the radius of the so called influence region. This region is the magma volume surrounding the bubble from which gas diffuse into the bubble. It is quite natural for a spherical bubble to consider this influence region as a sphere, centered at the bubble center and with radius Ŝ. The evolution of the radius Ŝ is obtained by assuming that the volume of the influence region is constant in time (see [START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF]), so that

Ŝ( t) = Ŝ3 0 + R( t) 3 1/3 , ( 4
)
where Ŝ0 is a constant representing the radius of the influence region when the bubble has a null radius. Finally, in equation ( 3), α = α( t) is the magma porosity, or void fraction. Porosity is actually a macroscopic notion, but as far as it can be given a meaning for a single bubble, we choose, following [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF][START_REF] Proussevitch | Dynamics and energetics of bubble growth in magmas: Analytical formulation and numerical modeling[END_REF][START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF], to define it as

α( t) = R( t) 3 Ŝ( t) 3 .
We rewrite equation ( 2) as a differential equation on the radius R:

Ṙ( t) = R( t) 4η ef f P ( t) -Pa ( t) - 2σ R( t) . (5)
To obtain the internal pressure P ( t), or equivalently by (1) the gaseous water mass M ( t), we need to write the mass conservation of water. First we consider the volatile mass balance at the bubble-magma interface, which reads in spherical geometry 4π 3

d d t ρ R3 = 4π ρm F R( t) , (6)
where ρ m is the magma density and F R( t) represents the water flux from the magma into the bubble at the interface. Now, we introduce the concentration of water in the melt, which is a function C = C(r, t) of the time t and on the radial distance from the boundary of the bubble, r ∈ [ R( t), Ŝ( t)]. With this notation, the flux F = F (r, t) in equation ( 6) is given by

F R( t) = D R2 ∂C ∂r r= R( t) , (7)
where the diffusion coefficient D = D(r) may depend on the distance r via the concentration C, see for example [START_REF] Proussevitch | Dynamics and energetics of bubble growth in magmas: Analytical formulation and numerical modeling[END_REF][START_REF] Proussevitch | Bubbledrive-1: A numerical model of volcanic eruption mechanisms driven by disequilibrium magma degassing[END_REF], but is assumed to be constant in all the numerical simulations we have performed.

The definition of the influence region implies that the total water mass inside it, that is the sum of the water mass in the bubble and of the water dissolved in its influence region, must remain constant in time, see [START_REF] Welch | Direct simulation of vapor bubble growth[END_REF]. Assuming that for a bubble of radius zero the water concentration in magma is a constant C 0 , this may be expressed in the following form: 4π ρ( t) 3 R3 ( t) + 4π ρm

Ŝ( t) R( t) r2 C(r, t) dr = 4π ρm 3 Ŝ3 0 C 0 . (8)
Within the influence region, that is for radii r ∈] R( t), Ŝ( t)[, the water concentration is assumed to follow an advection-diffusion equation

∂C ∂ t + v m ∂C ∂r = 1 r2 ∂ ∂r r2 D ∂C ∂r ,
where v m is the radial velocity in the melt and is obtained by solving the continuity equation in radial form and considering the incompressibility of the melt (see [START_REF] Scriven | On the dynamics of phase growth[END_REF]):

v m = Ṙ R2
r2 . The advection-diffusion finally reads

∂C ∂ t + Ṙ R2 r2 ∂C ∂r = 1 r2 ∂ ∂r r2 D ∂C ∂r , (9)
and has to be complemented with boundary conditions. For r = R, that is at the bubble-melt interface, the pressure has to be in equilibrium with water concentration, following Henry's law

C( R, t) = K H P , (10)
where K H is the Henry constant. The other boundary condition is given at the external interface of the influence region, and follows from the global mass balance [START_REF] Proussevitch | Dynamics and energetics of bubble growth in magmas: Analytical formulation and numerical modeling[END_REF]. Indeed, stating that the time derivative of (8) has to be zero, a straightforward computation taking into account ( 9) and [START_REF] Rosebrock | Study of the growth and motion of graphitic foam bubbles[END_REF] shows that the water flux on the boundary of the influence region, r = Ŝ, is null:

∂C ∂r r= Ŝ = 0, (11)
Summarizing, we are lead to the following system of differential equations

Ṙ = R 4η ef f P -Pa - 2σ R , (12)
Ṁ = 3ρ m D R2 ∂C ∂r r= R ,
where C solves the advection-diffusion equation [START_REF] Proussevitch | Bubbledrive-1: A numerical model of volcanic eruption mechanisms driven by disequilibrium magma degassing[END_REF], with boundary conditions [START_REF] Toramaru | Numerical study of nucleation and growth of bubbles in viscous magmas[END_REF] and [START_REF] Blower | Coupling of viscous and diffusive controls on bubble growth during explosive volcanic eruptions[END_REF].

Dimensionless problem

In the preceding model, the physical parameters involved may vary for several orders of magnitude and in a very intricate manner. Table 1 recalls their meaning and presents a sample of these values, which come from the experimental results quoted below. The behavior of the solutions to the model can vary drastically with these values, from one experimental situation to another, and computational times may be very long. The aim of this section is to provide a dimensionless set of equations, in order to identify several specific regimes and to eventually give analytical solutions to each regime.

Following [START_REF] Proussevitch | Dynamics and energetics of bubble growth in magmas: Analytical formulation and numerical modeling[END_REF][START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF][START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF], a set of five dimensions gives a physically relevant scaling, namely a bubble radius, a pressure, a density, a viscosity coefficient and a diffusion coefficient. These characteristic dimensions are chosen here as the corresponding initial values: the initial bubble radius, R i , the initial gas density ρ i , the initial ambient pressure P i = P a (t = 0) and the diffusion and viscosity coefficients, D i and η i , leading to the following scalings on variables and parameters:

η ef f = ηeff η i D = D D i ρ m = ρm ρ i Σ = 2σ R i P i R = R R i S = Ŝ R i r = r R i Ṙ = ṘP i R i ∆P P = P P i P a = Pa P i ρ = ρ ρ i t = t∆P P i
The time evolution equation of the radius of the influence region (4) becomes

S 3 (t) = S 3 0 + R 3 (t). (13)
A convenient model for Pa (t), quite simple and compatible with experimental conditions, is a constant linear decompression function, which in dimensionless variables reads:

P a = 1 -t (14)
This ambient pressure function is the one considered in [START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF], and we shall apply it in the numerical simulations.

In this paper, we assume the temperature T to be constant, then the perfect gas law reads M = ρR 3 and straightforward computations show that the mass conservation equation ( 8) becomes

M + 3ρ m S(t) R(t) r 2 C dr = S 3 0 C 0 ρ m , (15)
and that the system of differential equations on mass and radius [START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF] rewrites

Ṙ = R Θ V η ef f P -P a - Σ R , (16)
Ṁ = 3ρ m Θ D r 2 D ∂C ∂r r=R , (17)
where we have introduced the relaxation parameters, see [START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF]:

Θ V = 4η i ∆P P 2 i , Θ D = R 2 i ∆P D i P i . (18)
Finally, the water concentration in the melt satisfies the dimensionless advectiondiffusion equation in the influence region

∂ t C + ṘR 2 r 2 ∂ r C = 1 Θ D 1 r 2 ∂ r r 2 D∂ r C (19)
with the following boundary conditions:

C(R, t) = C H √ P , ∂C ∂r r=S = 0, ( 20
)
where C H is given by C H = K H √ P i . The dimensionless parameters Θ D and Θ V defined by ( 18) are driving the behavior of all the equations involved, and will be referred to as the viscosity Θ V and diffusion Θ D relaxation parameters. Their values may vary of several order of magnitude and define several specific regimes which are analyzed in section 4 below.

At this stage, it is worth to focus on initial conditions. From the above adimensionalization, the initial radius is R = 1, and the corresponding initial pressure is assumed to be at equilibrium with the ambient pressure and the surface tension: P (0) = (P i + Σ)/Θ V . In experiments, after the first decompression jump ensuring bubble nucleation, we actually wait for the bubble mass and radius to reach an equilibrium state, so that for instance bubble radius stops to evolve. The initial water concentration in the influence region is C(0, r) for r ∈ [R(0), S(0)]. We choose a constant initial concentration C(0, r) = C i . The behavior of the solutions drastically depends on the relationships between the initial concentration and the initial bubble pressure. A particular case is given by

C(0, r) ≡ C i = C H P (0), r ∈ [R(0), S(0)]. (21)
These "well-prepared" initial data correspond physically to some equilibrium between the inner bubble pressure and the water concentration. The solution behaves nicely as expected: the radius increases, the inner pressure decreases. On the other hand, when (21) is not satisfied, one can observe some jumps at small times (depending on the value of Θ V ), where bubble pressure and water concentration try to reach equilibrium. If C i > C H P (0), there is a water excess in the melt, so that the radius may decrease or the pressure increase very fast before reaching a smooth regime. On the contrary, no bubble should exist when the water is not sufficient, C i < C H P (0), and the solution may exhibit almost discontinuous behavior, or even not exist.

Numerical approximation

In this section we consider the numerical approximation of the model ( 16)-( 17)-( 19), on R(t) and M (t) and C(r, t), together with the boundary conditions (20) and the external assumption ( 13)-( 14). We propose a numerical scheme for the advection-diffusion equation which conserves exactly the water mass by construction. This is a delicate point of the discretization; the flux at the bubble border has to be carefully computed because the magnitude of the relaxation parameters Θ V and Θ D may differs of several orders. We shall first present the straightforward discretization of the system of differential equations ( 16)-( 17) and then present the more delicate discretization of the advection-diffusion equation, (19).

The differential system

We describe here the basic elements of the numerical scheme for the differential system ( 16)- [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF]. Let us first define, for n ∈ N, the time t n+1 and the ambient pressure P n+1 a = P a (t n+1 ), at the iteration step n + 1, respectively by:

t n+1 = t n + ∆t n , P n+1 a = 1 -t n+1 ,
where the time step ∆t n is computed at each iteration and must satisfy some stability conditions which will specified later on, see section 3.3. The numerical results, see section 5, will be plotted in term of the ambient pressure P a , which may be considered as a time variable. We choose a semi-implicit scheme for the discretization of ( 16), in the sense that the discrete bubble radius R n = R(t n ) is treated implicitly, whereas the pressure P n = P (t n ) remains explicit. Thus the evolution of the discrete radius is given by

R n+1 = R n -∆t n Σ Θ V η n ef f 1 -∆t n (P n -P n a ) Θ V η n ef f -1 . ( 22
)
Next we discretize the equation for the mass balance at the magma-bubble interface, [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF], by a semi-explicit scheme. Defining the discrete bubble mass by M n = M (t n ), we recall that the pressure is given for all n by P n (R n ) 3 = M n . We denote by F n 0 = F (R, t n ) the discrete flux at the interface r = R, and we set

M n+1 = M n + 3ρ m ∆t n Θ D F n+1 0 . (23)
Finally, the discrete radius of the influence region, S n = S(t n ) is given by

S n = (R n ) 3 + S 3 0 1/3 .
We turn now to the definition of the discrete flux F n+1 0 , which follows from the discretization of the advection-diffusion equation (19).

The advection-diffusion discretization

The advection-diffusion equation for the water concentration C(t, r) is solved by splitting the equation between the advection step and the diffusion step. The advection step consists in discretizing the following transport equation:

r 2 ∂ t C + ṘR 2 ∂ r C = 0. ( 24
)
We choose to solve it by a Lagrangian method, namely, a set of mesh points at time t n being given, r n i , 0 ≤ i ≤ N , we compute a new mesh at time t n+1 by solving explicitly the equation of characteristics

r 2 i dr i dt = R 2 Ṙ, (25)
which integrates in

r n+1 i = R n+1 3 -(R n ) 3 + (r n i ) 3 1/3 . (26)
The above relation defines the mesh for all n ≥ 1 as soon as the initial discretization r 0 i , 0 ≤ i ≤ N , is fixed. Inspection of (26) shows that, if r 0 0 = R 0 and r 0 N = S 0 are known, we have for all n ≥ 1: r n 0 = R n and r n N = S n , so that any choice R 0 < r 0 1 < • • • < r 0 N -1 < S 0 is relevant, in particular the uniform grid defined by r 0 i = (i/N )S 3 0 1/3 . Finally, we denote by ∆r n i the (non uniform) space discretization step, ∆r n i = r n i+1 -r n i . The diffusion step consists in discretizing the equation

r 2 d t C = 1 Θ D ∂ r r 2 D∂ r C .
Following a standard finite volume strategy, we integrate the above equation on the mesh [t n , t n+1 [×]r n+1 i , r n+1 i+1 [, looking for piecewise constant solutions C n i on the mesh. We obtain, for i = 1, ..., N -1:

C n+1 i = C n i + ∆t n Θ D 3 F n+1 i+1 -F n+1 i r 3 i+1 -r 3 i , (27)
where F n+1 i stands for the discrete flux between cells i -1 and i. As an approximation of r 2 D∂ r C for r = r n+1 i , we choose a centered finite difference:

F n+1 i = D n i 2 r n+1 i 2 r n+1 i+1 -r n+1 i-1 C n i -C n i-1 . (28)
The boundary conditions (20) become

C n+1 0 = C H √ P n , F n+1 N = 0. ( 29
)
We are now in position to close equation ( 23) by setting the value of the discrete flux at the bubble-magma interface, F n+1 0 . The trick here is to obtain a discrete analogue of the mass preservation [START_REF] Jaupart | Gas content, eruption rate and instabilities of eruption regime in silicic volcanoes[END_REF]. Defining the discrete total water mass M n by

M n = M n + ρ m N -1 i=0 C n i r 3 i+1 -r 3 i , (30)
we have the following result.

Proposition 1 Let F n+1 0 be given by

F n+1 0 = F n+1 1 - Θ D ∆t n r 3 1 -r 3 0 3 (C n+1 0 -C n 0 ). ( 31
)
Then the numerical scheme (22)-(31) conserves the discrete total water mass, that is M n = M 0 , ∀n ∈ N.

Proof:The total water mass at time t n+1 is equal to:

M n+1 = M n+1 + ρ m N -1 i=0 C n+1 i r 3 i+1 -r 3 i .
Applying (23), splitting the sum for i = 0 and i = 1...N -1 and replacing (27) in the sum, we obtain:

M n+1 = M n + 3ρ m ∆t n Θ D F n+1 0 + ρ m C n+1 0 (r 3 1 -r 3 0 )+ +ρ m N -1 i=1 C n i (r 3 i+1 -r 3 i ) + 3ρ m ∆t n Θ D N -1 i=1 F n+1 i+1 -F n+1 i .
Hence, simplifying the last sum we get:

M n+1 = M n + 3ρ m ∆t n Θ D F n+1 0 + ρ m C n+1 0 (r 3 1 -r 3 0 )+ +ρ m N -1 i=1 C n i (r 3 i+1 -r 3 i ) + 3ρ m ∆t n Θ D (F n+1 N -F n+1 1
).

Recalling that F n N = 0 for all n, and splitting, in the definition M n , the sum for i = 0 and i = 1...N -1, we must have that:

3ρ m ∆t n Θ D F n+1 0 + ρ m C n+1 0 (r 3 1 -r 3 0 ) -3ρ m ∆t n Θ D F n+1 1 = ρ m C n+1 0 (r 3 1 -r 3 0 ),
which is verified since F n+1 0 is defined by (31). Remark Notice that Proposition 1 holds true for any choice of the discrete flux in (28) for 1 ≤ i ≤ N -1.

Stability conditions

In this section, we describe how to compute for each n ≥ 1 a time step ∆t n ensuring some stability conditions on the numerical approximations. The idea is to compute a bound for ∆t n for each numerical approximation (radius, mass, concentration), and then to take as ∆t n the minimum of these three stability conditions.

We cannot obtain a completely satisfactory stability condition ensuring the bubble mass positivity in (23). However, a partial condition is given by asking that the discrete bubble pressure P n remains larger than the ambient pressure P n a at each iteration n. This leads to

∆t n < Θ V |(P n a -P n )(R 3 ) n -3ρ m (r 3 1 -r 3 0 )C n 0 | |3ρ m F n 1 | (32)
This condition yields reasonable time steps in almost all cases, namely ∆t n ∼ 10 -9 instead of 10 -12 when both Θ V and Θ D are small. This condition is not sufficient to avoid oscillations in the solution, in particular when the initial conditions are not well prepared in the sense given above, see equation (21). In some rare case, the oscillations blow up, but this is consistent with the physical incompatibilities concerning (21). The scheme is, nevertheless, robust in the sense that these oscillations, which appear at the beginning of the computation, tend to disappear when time increases.

We look now for stability conditions for ( 22) and ( 27). In particular we search a bound of the time step ∆t n such that each solution is positive.

Proposition 2 Assume that for n ∈ N ∆t n < min R n Θ V η n ef f Σ , Θ V η n ef f |P n -P n a | . ( 33
)
Then the solution to the numerical scheme (22) is positive, i.e. R n > 0.

Proof: Let us assume that at the iteration n all the variables are positive, then at the iteration n + 1, using (22), ∆t n must be such that:

R n -∆t n Σ Θ V η n ef f 1 -∆t n P n -P n a Θ V η n ef f -1 > 0.
We have two possibilities. The first one is when P n -P n a ≤ 0. Then we have:

1 -∆t n P n -P n a Θ V η n ef f > 0.
Hence,

R n -∆t n Σ Θ V η n ef f > 0 which implies, ∆t n < Θ V η n ef f R n Σ .
The second one is: P n -P n a > 0. Then ∆t n must be the positive solution of a second order equation with a positive dominant coefficient:

Σ(P n - P n a )/(Θ V η n ef f ) 2 . Hence R n+1 is positive, if ∆t n is external to the roots: Θ V η n ef f P n -P n a , R n Θ V η n ef f Σ .
Finally, we remark that the choice (33) verifies both conditions.

Proposition 3 Assume that, for n ∈ N,

∆t n < Θ D 6 min i (r 3 i+1 -r 3 i ) r n+1 i+1 -r n+1 i-1 D n+1 i (r n+1 i ) 2 . ( 34
)
Then the solution to the numerical scheme (27) is positive, i.e. C n i > 0 for all i = 1, . . . , N -1.

Proof: We first remark that the Dirichlet condition on the boundary r = R implies that C n 0 > 0 for all n ∈ N. Assuming that at the n th iteration C n i is positive for all i = 0...N -1, we want ∆t n to verify C n+1 i > 0. Thus from (27), it must be, for i = 1...N -1:

C n i + 3∆t n Θ D (r 3 i+1 -r 3 i ) F n+1 i+1 -F n+1 i > 0.
i is given by (28), collecting the terms with respect to C n i-1 , C n i and C n i+1 and considering that C n i are positive for all i, we get a sufficient condition for the positivity of C n+1 i in the form, for i = 1...N -2:

1 - 6∆t n Θ D (r 3 i+1 -r 3 i ) D n+1 i+1 (r n+1 i+1 ) 2 r n+1 i+2 -r n+1 i + D n+1 i (r n+1 i ) 2 r n+1 i+1 -r n+1 i-1 > 0. Since r n+1 i+2 -r n+1 i D n+1 i+1 (r n+1 i+1
) 2 > 0, the time step ∆t n given by (34) verifies the above condition. If now i = N -1, recalling that F n N = 0 for all n, we obtain that ∆t n must satisfy:

∆t n < Θ D 6 (r 3 N -r 3 N -1 ) r n+1 N -r n+1 N -2 D n+1 N -1 (r n+1 N -1
) 2 which concludes the proof.

Limit cases

As mentioned in section 2, the system of equation ( 16)-( 17)-(19) has two relaxation times, Θ V and Θ D defined by [START_REF] Hess | Viscosities of hydrous leucogranitic melts: a non-Arrhenian model[END_REF], which may differ by several order of magnitude, depending on the values of, for instance, diffusivity or viscosity. In many experiments Θ V and/or Θ D are very small, of the order of 10 -7 . The time steps ∆t n depending on these values, the computational time needed to reach a porosity close to 1 is very large. The study of the limit cases, such as when Θ V and Θ D tends to ∞ or to 0, is thus an attractive alternative to solving the full system because it leads to simplified models with smaller simulation times. In particular, we will classify the different limits in regimes of bubble growth by considering the ratio Θ V /Θ D . Following [START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF], we define a viscous regime when the ratio is very small, see 4.1, an equilibrium regime when the ratio is of order 1, see 4.2, and a diffusive regime when the ratio is large, see 4.3. At the end of each section, we will also summarize when necessary the numerical scheme corresponding to the simplified cases. As we are mainly interested in the behavior of the bubble physical dimensions (pressure P and radius R), we shall only describe how to compute these two quantities. In particular, we recall that, when comparing with experiments, we consider the porosity: α = R 3 /S 3 , and we note that all the simulations performed in this paper tend to compute until a porosity α as close as possible to 1, even if it is not realistic from a physical viewpoint: the crossover value above which our model of bubbles in a melt is no longer valid is around 0.7.

We list all the possible limits in table 2, which references each simplification to compute the bubble radius R and mass M (or pressure P ). We note that these simplifications are relevant only under the assumption that temperature is constant, so that the perfect gas law reads M = P R 3 . We give when possible the expression of the simplified solutions for variable viscosities and diffusivities. In table 2 we also give the orders of magnitude delimiting each regime.

0 10 -5 10 1 ∞ Θ V 4.1 Viscous regime : Θ V /Θ D << 1
We first consider the case when the viscous relaxation parameter is smaller than the diffusion one. There are three possibilities: Θ V tends to zero and Θ D is of order 1 or tends to infinity; and Θ V is of order 1 and Θ D tends to infinity.

Θ

V ∼ 1 and Θ D → ∞
Since Θ D is very large, equation (19) reads r 2 d t C(r, t) = 0, which yields to:

S(t) R(t) r 2 C(r, t) dr = S(0) R(0) r 2 C 0 dr. ( 35
)
Concerning the bubble mass evolution, since Θ D >> 1 the water mass inside the bubble is constant and equation ( 17) reads Ṁ = 0, which is equivalent to:

M = M (0). ( 36
)
In fact, a large diffusive relaxation parameter may be physically given by a very small value for the diffusivity in magma; hence there will not be diffusion of water from the magma into the bubble, and the bubble water mass will not change.

As Θ V ∼ 1, no simplification is possible for the equation giving the evolution of the bubble radius [START_REF] Gardner | Bubble growth in highly viscous silicate melts during continuous decompression from high pressure[END_REF]. Still, recalling that M (0) = M = P R 3 and that the ambient pressure is given by ( 14), we can write a differential equation only depending on R:

Ṙ = 1 Θ V η ef f M (0) R 2 -R P a -Σ , (37)
or equivalently on P :

Ṗ = - 3P Θ V η ef f P -P a -Σ P M (0) 1/3
. This result is identical to that obtained in [START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF]. Equation (37) can be easily solved with an implicit scheme.

4.1.2 Θ V → 0 and Θ D → ∞
When Θ D is large, the simplifications ( 35) and ( 36) are always true. In particular, water mass in the bubble is constant: M = M (0). Regarding equation ( 16), since Θ V is very small, multiplying by Θ V , recalling that M (0) = P R 3 , and simplifying allows us to obtain a third order equation on R: R 3 P a + ΣR 2 -P (0) = 0, which admits an unique real solution given by the explicit relation

R = c 6a + 2b 2 3ca - b 3a , ( 38
)
where a = P a /P (0), b = Σ/P (0) and c = (108a

2 -8b 3 +12 √ 3 √ 27a 2 -4b 3 a) 1/3 . 4.1.3 Θ V → 0 and Θ D ∼ 1
If Θ D ∼ 1, the simplifications of section 4.1.1 no longer hold. There is, for instance, no possible simplification for equations ( 17) and ( 19). Nevertheless, since Θ V << 1, we have:

P = P a + Σ R (39)
which links the bubble pressure P to the radius R. On the other hand, considering that M = P R 3 and equation ( 17), we obtain the following differential equation for R:

Ṙ = 3ρ m Θ D r 2 D ∂C ∂r r=R + R 3 2ΣR + 3R 2 (1 -t) -1 , ( 40
)
where the water concentration is obtained solving the advection-diffusion equation ( 19). Numerically, we compute M n+1 and C n+1 i using respectively equations ( 23) and ( 27), then we apply relation M = P R 3 to compute P and finally from (39) we obtain R.

Equilibrium regime:

Θ V /Θ D ∼ 1
In this section we deal with those regimes in which the relaxation parameters Θ V and Θ D are of the same order of magnitude. More precisely, when both Θ V and Θ D tend to zero, or to ∞, since otherwise no simplification is possible.

4.2.1 Θ V → 0 and Θ D → 0
This is an interesting situation because on the one hand computational time is very long and on the other hand it corresponds to the so-called equilibrium growth, which is a common situation in natural magmas: the bubble is always at its maximum possible radius. First notice that, since Θ V << 1, following the discussion of section 4.1.3, we have the simplification (39).

Next, let us consider Θ D << 1, then multiplying the water concentration equation ( 19) by Θ D and simplifying, we have:

1 r 2 ∂ ∂r r 2 D ∂C ∂r = 0,
for which the solution, taking into account the boundary conditions (20), reads:

C(r, t) = C(R, t) = C H √ P , ∀ r ∈ [R(t), S(t)], (41)
When Θ D tends to zero, equation ( 17) is no longer valid to compute the water mass variation inside the bubble. Therefore, we consider the total mass conservation equation [START_REF] Jaupart | Gas content, eruption rate and instabilities of eruption regime in silicic volcanoes[END_REF]. Replacing C by ( 41) in [START_REF] Jaupart | Gas content, eruption rate and instabilities of eruption regime in silicic volcanoes[END_REF] and recalling [START_REF] Chouet | Dynamics of diffusive bubble growth and pressure recovery in a bubbly rhyolitic melt embedded in an elastic solid[END_REF], we obtain

M = ρ m S 3 0 (C 0 -C H √ P ), (42)
Now since M = P R 3 , equation ( 42) is a second order equation in X = √ P , which turns out to have one positive solution, namely

P = -Y + Y (Y + 4R 3 C 0 /C H ) 2R 3 2
, where

Y = C H ρ m S 3 0 . ( 43
)
Hence R and P (or M ) are uniquely defined by combining (39) with ( 43) or (42) (which we have used in our numerical resolution).

Equation ( 42) allows one to establish an expression for the porosity α = α(t). Recalling that the porosity is defined by α = R 3 /S 3 , we have

R 3 = S 3 0 α 1 -α .
Replacing M by P R 3 in (42) gives readily

P α 1 -α = ρ m C 0 -C H √ P ,
which in turns leads to

α = β P + β , where β = ρ m C 0 -C H √ P . (44)
Expressing (44) in dimensional gives:

α = γ M w P + γ , where γ = RT ρm K H P0 -P , ( 45
)
where P0 is the pressure of a bubble of radius zero. Equation ( 45) is equivalent to the most commonly used formula to calculate porosity in the equilibrium regime (e.g. [START_REF] Jaupart | Gas content, eruption rate and instabilities of eruption regime in silicic volcanoes[END_REF]). We underline, however, that in the common formula P = Pa , whereas the pressure in (45) accounts for surface tension because P = Pa + 2σ/ R.

4.2.2 Θ V → ∞ and Θ D → ∞
This situation could result from a very large viscosity, which yields Θ V >> 1, combined with a very small diffusivity, which yields Θ D >> 1. Following experimental evidence (see [START_REF] Gardner | Bubble growth in highly viscous silicate melts during continuous decompression from high pressure[END_REF]), we can imagine that the physical system is "fixed" or "frozen".

On the one hand, as discussed in section 4.1.1, the water mass in the bubble is constant, see equation (36). On the other hand, since Θ V >> 1, from equation ( 16) we also obtain Ṙ = 0, that is:

R(t) = R(0) = 1. ( 46
)
Since both the mass M and the radius R are constant, the pressure P is explicitly determined by M = P R 3 , and no numerical scheme is needed.

Diffusive regime :

Θ V /Θ D >> 1
In this last section, we treat regimes which have the viscous relaxation parameter larger than the diffusion one. We have to differentiate three cases: when Θ D is small and Θ V is of order one or tends to infinity, and when Θ D is of order one and Θ V tends to infinity.

4.3.1 Θ V → ∞ and Θ D ∼ 1
As shown before, when Θ V → ∞, we obtain equation (46) and the bubble radius is constant in time. Recalling that M = P R 3 , we have Ṁ = Ṗ , hence from [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF] we get the following differential equation on P :

Ṗ = 3ρ m Θ D r 2 D ∂C ∂r r=R , (47)
with C solution of (19). Numerically, the radius R is constant, and the pressure P is computed using the numerical scheme of section 3.

4.3.2 Θ V → ∞ and Θ D → 0
Considering the discussion of sections 4.2.2 and 4.2.1, both simplifications (46) and (43) hold. From (46), the bubble radius is constant, R = 1, so that the bubble pressure is also constant, and is explicitly obtained by simplifying equation (43):

P = -Y + Y (Y + 4C 0 /C H ) 2 2 , Y = C H ρ m S 3 0 . (48) 4.3.3 Θ V ∼ 1 and Θ D → 0 Since Θ V ∼ 1,
there is no possible simplification for the radius equation ( 16), but, from section 4.2.1, the bubble pressure is computed by (43). In this case, the radius R must be computed using the numerical scheme of section 3, while the pressure P is explicitly given by (43).

Numerical results

In this section we compare first the numerical results obtained using the general scheme of section 3 to those obtained with the numerical approximation of the simplified schemes of section 4. This is followed by a comparison between the behavior of our numerical results and experimental data described in [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF]. Finally, we consider the reference code described in [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF] and compare its numerical results with those obtained using the general scheme. Le us first discuss the dependence of our results on the number of discretization points N with respect to the radial variable r. Various numerical tests show that a small number of points is sufficient in order to well capture the behavior of the discrete flux F n 0 on the bubble surface. The relative errors for the bubble radius R, the bubble pressure P , and the porosity α with respect to the reference ones with N = 2500, are of order 10 -3 for N = 50 and of order 10 -4 for N = 250, respectively. Therefore we choose N = 50 in all the following computations.

In figure 1 we show the evolution of the concentration function C(t, r) computed solving the general scheme with Θ V = 0.000236 and Θ D = 5.28929. Bubble size evolution is sketched as grey circles of increasing radius R. We clearly see the mesh refinement near the bubble wall (the grey circle portion) when the concentration function becomes stiffer.

Numerical convergence

We show now through a few selected examples the numerical convergence of the global numerical scheme defined in section 3 towards the simplified limit cases discussed in section 4. Convergence is determined by fixing either Θ V or Θ D and varying the other one. In figure 2 we show the convergence for the bubble radius towards selected limit cases. On the left, we plot the bubble radius evolution with respect to the ambient pressure P a which linearly decreases in time, fixing Θ D = 0.1 and varying Θ V from 10 -5 to 10 1 , together with the radius computed as explained in 4.1.3 or just defined as the constant 1, as justified in 4.3.1. On the right, we plot the bubble radius evolution with respect to time, fixing Θ V = 0.1 and varying Θ D from 10 -3 to 10 3 , together with the radius obtained as explained in 4.1.1 and 4.3.3. In both cases we can observe the transition of the general unsimplified regime (the middle case in table 2) from one simple growth regime to the next.

Experimental data vs. numeric

Controlled decompression experiments on high temperature magmas are able to produce gas bubbles. By varying the end pressure, data on bubble size and porosity have been retrieved for different initial conditions such as magma chemical composition, temperature, decompression rate, etc. In this work we illustrate how comparisons between such experimental data and model outputs can be carried out. Let us first present here the experimental framework used in [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF]. Samples of viscous magma are saturated in water and maintained under pressure for about 5 days for the water to be homogeneously dissolved into the magma. Then an instantaneous decompression gives rise to bubble nucleation. After waiting for a few minutes that these initial, small bubble reach their equilibrium sizes, a linear decompression is applied until a final pressure where samples are quenched by a sudden cooling to ambient temperature. The cold samples are then sliced and analyzed to obtain bubble sizes and porosity. One experiment quenched just after the sudden decompression that nucleates the bubbles gives the initial conditions for our model.

The physical values used and measured during the experiments are the following: the initial radius R i = 17.5 • 10 -6 , the diffusion coefficient D i = 5.79 • 10 -12 , the initial concentration C 0 = 3.44 • 10 -2 , the initial pressure P i = 10 8 , the surface tension σ = 0.1, the viscosity η i = 5.9 • 10 -4 , the magma density ρ m = 2400, the gas porosity density ρ i being calculated by the perfect gas law (1), the temperature T = 1098, and the decompression rate ∆P = 10 -5 . We tested two different experimental series. In the first series, bubbles growth was only due to gas expansion and water diffusion form the magma. In the second series, bubbles growth was also due to coalescence processes. The porosity evolution of both series is comparable but the evolution of bubble size differ.

In figure 3 we show the evolution of the porosity α and of the radius R with respect to the ambient pressure P a . On the left graph are represented three numerical results for different viscosity calculation and initial porosities along with the experimental results obtained in [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF] (triangles). The run represented by the grey line had a constant viscosity η, whereas the black line and the dashed line had variable viscosity η ef f computed applying the formula (3). We remark that considering a variable viscosity η ef f instead of a constant one has an impact on the numerical result only when the ambient pressure becomes very small because the grey and the black curves diverge only when P a < 0.1. The dashed line is a numerical result computed starting from the porosity measured on the experiment quenched just after the sudden decompression (porosity α = 0.056). The other two runs use the equilibrium state in Section 4.2.1, which predicts a porosity of 0.0779 instead of 0.056. We note that the experimental points fit better the dashed line at high P a and are closer to the two other curves at low P a . This leads us to conclude that during the first phase of the experiments the time between nucleation and the beginning of the decompression was not enough to reach the equilibrium. The right graph shows the numerical radii R as a function of the ambient pressure P a for the same three numerical runs. Experimental results are now represented by squares centered on the median value of the experimental radii and a standard deviation representing the spread of measured bubble radii. We note that the three numerical runs are very similar, regardless of initial conditions, and that the fit between experiments and numerical results is worse for radius than porosity. The larger misalignment of experimental radii compared to that of porosity is explained by the fact that each experimental point in figure 3

Discussion of a particular experiment

In figure 4 we compare the experimental results obtained in [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF] when bubble coalescence occurs to our numerical results with variable viscosity for three different initial porosities. The left graph, which displays the evolution of porosity with ambient pressure, indicates that this experimental series is best reproduced numerically by starting from a initial porosity of 0.035. This is lower than the best fit value of the other series. A tentative explanation of this situation is an even shorter delay between bubble nucleation and the start of the decompression. We plotted on the right graph, which shows the evolution of R with ambient pressure, distribution histograms of measured bubble radii. The three computed bubble radii R all fit the experimental measurements within a standard deviation, but, considering the large spread of bubble sizes, it is not possible to choose which numerical results has the better fit. As discussed in [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF], the poly-disperse nature of bubble growth in the experiments was caused by bubble coalescence. Hence, one should consider poly-disperse modeling of the bubble population in order to produce more accurate results that resolve the spread in bubble sizes.

Tasting various decompression rates

The behavior of the proposed numerical scheme when considering linear decompression has been discussed in figure 3. Here we present some numerical test ensuring the robustness of the numerical scheme when considering nonlinear decompression of the ambient pressure P a (t). We have performed one test that considers decompression by pressure jumps, and one test that considers an increasing ambient pressure, i.e. ∆P < 0. In figure 5, we show the behavior of the porosity (on the left) and radius (on the right) for the run with pressure jumps. The physical values we have considered are: R 0 = 4.4 • 10 -6 , α i = 0.041, D = 2 • 10 -12 , C i = 0.036, P i = 7 • 10 7 , σ = 0.1, ρ m = 2354, T = 998, a variable viscosity η ef f and ∆P = 2.5 • 10 5 for the linear decompression; whereas, for the nonlinear decompression, we have a series of instantaneous pressure jumps every 20 seconds. In the increasing ambient pressure test, we have first linearly decompressed the ambient pressure from 1 to 0 and then linearly compressed it back, with the same decompression rate (in absolute value). There is a very good agreement between the numerical results of the decompression phase and the compression one. For instance the relative error between the initial radius and the one at the end of the decompression/compression cycle is of order 10 -3 .

Comparison with existing code

Lastly, we compare our numerical code with the one of reference described in [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF]. The main differences between the two schemes are: in [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF] temperature and diffusion are not constant, as it is the case in our study; vitrification and crystallization are taken into account and we do not consider them; the viscosity function is not the same in the two schemes. We find that computational times of the [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF] code are much smaller than those of our code. In figure 6 we show the radius evolution with respect to the ambient pressure and computed with the two codes starting from the same initial data. The small difference at the beginning of the computation (when the ambient pressure is 1) is due to the lack of control of the time discretization in the code of [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF], so that the results miss some initial time points to better capture the pressure jump. The difference at final time (when ambient pressure is smaller than 0.1) may be due to the difference in the viscosity formula, or to crystallization or vitrification. Overall, both numerical results are in good agreement. Further comparisons with other experimental results are under investigation and will be exposed in a future work.

Figure 6: Bubble radius as a function of the ambient pressure P a : PSA-line is the numerical result using the reference code [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF].

Conclusions

We have proposed and applied a numerical scheme for the approximation and simulation of the solution of a non-linear system of differential equations coupled with an advection diffusion equation, previously proposed in the volcanology literature (see for example [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF], [START_REF] Proussevitch | Dynamics and energetics of bubble growth in magmas: Analytical formulation and numerical modeling[END_REF] and [START_REF] Lensky | Bubble growth during decompression of magma : experimental and theorical investigation[END_REF]). Our goals were to give a discretization of the system that is conservative by construction and to study the asymptotic limits when the two relaxation parameters Θ V and Θ D tend to 0 or ∞.

In the recent years (see for example [START_REF] Chouet | Dynamics of diffusive bubble growth and pressure recovery in a bubbly rhyolitic melt embedded in an elastic solid[END_REF] and the reference therein), the numerical approach to solve the model governed by equations ( 16)-( 17)-( 19) is based on the one proposed in [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF]: the transport term in the advectiondiffusion equation is simplified by the means of a change of variable at the continuum level, leading to an heat equation with non standard diffusion term. Nevertheless, with the method proposed in [START_REF] Proussevitch | Dynamics of diffusive bubble growth in magma: isothermal case[END_REF], a large number N of discretization points in the radial direction is greatly reduced by the means of a variable mesh size. This size is controlled by an empirically defined parameter that ensures the conservation of water mass and that precisely captures the behavior of the flux on the bubble border. With our approach, a small number of points, N = 50, also guarantees precise results, but the mesh size is automatically defined. In fact, the discrete flux on the bubble border is defined in such a way that the numerical scheme preserves water mass (see Proposition 1). There is thus no need to adjust an empirical parameter to ensure scheme accuracy. This advantage is balanced by more strict stability conditions that yield small time steps and long computation times.

Concerning the asymptotic behavior of the coupled system of equations, we have analytically deduced simplified models in three regimes: viscous, when Θ V /Θ D << 1, diffusive, when Θ V /Θ D >> 1, and in equilibrium, when Θ V /Θ D ∼ 1. In particular, when both Θ V and Θ D tend to 0, we retrieve the equilibrium state of the coupled system of equations. We have numerically shown the convergence of the scheme towards the solutions of the three regimes when varying the relaxation parameters. We also determined numerically the boundaries between the various regimes.

We compared our numerical results with data obtained from decompression experiments of natural magmas. We performed numerical tests that verify the behavior of the code when considering non-linear decompressions with pressure jumps and decompression/compression cycles. This validation of the code gives also a feedback on the quality of experimental results. In particular, we show that, unlike originally assumed by the authors, decompressions in [START_REF] Burgisser | Experimental constraints on degassing and permeability in volcanic conduit flow[END_REF] started while bubbles were still growing, i.e. equilibrium was not reached. Finally, it appears that the simplified mono-disperse framework is not accurate enough to capture the spreading bubble size distributions produced by coalescence. We thus infer that an extension of the physical model to include a poly-disperse description for the bubbles population is a worthy pursuit.

Figure 1 :

 1 Figure 1: Bubble growth and gas concentration function in the influence region.

Figure 2 :

 2 Figure 2: Bubble radius evolution and convergence. Left: convergence towards the limit cases 4.1.3 and 4.3.1. Right: convergence towards the limit case 4.1.1 and 4.3.3.

Figure 3 :

 3 Figure 3: Porosity α (on the left) and radius R (on the right) with respect to ambient pressure P a : effect of a variable viscosity η ef f .

Figure 4 :

 4 Figure 4: Bubble radius R with respect to ambient pressure P a : monodisperse vs. poly-disperse.

Figure 5 :

 5 Figure 5: Evolution of the porosity α (on the left) and the radius R (on the right) with respect to a nonlinear decompression P a .

Table 1

 1

		: Physics constants values
	D	diffusivity	10 -12 (m 2 • s -1 )
	ηeff viscosity	10 4 (P a • s)
	M w molar mass for water	18.10 -3 (kg • mol -1 )
	G	perfect gases constant	8.3144 (J • mol -1 • K -1 )
	T	temperature	1098.15 (K)
	σ	surface tension	0.1 (J • m -2 )
	ρ m	magma density	2154 (kg • m -3 )
	K H Henry constant	3.44.10 -6 (kg -1/2 • m 1/2 )
	S 0	influence radius for R = 0	6.204.10 -5 (m)
	∆P decompression rate	10 5 (P a • s -1 )
	P i	initial ambient pressure	10 8 (P a)
	C 0	water concentration for R = 0 4.21 (wt.%)

Table 2 :

 2 Limit cases ∞

		4.1.2	4.1.1	4.2.2
		(38), (36)	(37), (36)	(46), (36)
	10 3
	Θ D	4.1.3		4.3.1
		(39), (40)	(16), (17)	(46),(47)
	10 -3
		4.2.1	4.3.3	4.3.2
		(39), (43)	(43), (16)	(46), (48)

  is a full decompression run starting from P a = 1. As a result, bubble nucleation dynamics occurring during the initial decompression step is only approximately similar from one experiment to the next.

	1							8
			Initial conditions
	0.8		α i =0.0779, η=1 α i =0.0779, η=variable		6
	0.6		α i =0.056, η=variable
	α						R	4
	0.4
	0.2							2
	0							0
	0	0.2	0.4	P a	0.6	0.8	1	0	0.2	0.4	P a	0.6	0.8	1

Acknowledgment

This work was partially founded by the ERC-starting grant DEMONS (n. 202844) under the European FP7.