2. Basic equations

Let us consider a large volume of a gas-saturated magmatic melt subject to rapid decompression from a certain initial pressure to a final pressure (assumed fixed). Let us assume that the decompression rate is so high that it can be considered instantaneous. Note that the solution obtained below can easily be generalised to the case of finite-rate decompression in the spirit of Chernov et al. (Reference Chernov, Davydov and Pil’nik2023b ), where the bubble growth was studied without taking into account barodiffusion.

Suppose a gas bubble is formed and starts to grow in the volume of the melt as a result of decompression. As shown earlier (Chernov et al. Reference Chernov, Pil’nik, Davydov, Ermanyuk and Pakhomov2018), the bubble growth process can be divided into three stages. The first one is inertial. At this stage, the gas pressure in the bubble is almost equal to the initial pressure, and the time dependence of the bubble radius is exponential. The second stage is transient, characterised by a drop of gas pressure in the bubble to the pressure in the surrounding liquid. This leads to emergence of diffusion and, as result, to an increase of the mass of gas in the bubble. Finally, the third stage is asymptotically diffusive. At this stage, the gas pressure in the bubble is nearly equal to the pressure of the surrounding liquid, and the bubble growth occurs solely as result of diffusion.

Although the growth of gas bubbles in supersaturated solutions is well studied, the specific features of this process in magmatic melts should be carefully addressed. In particular, magmatic melts have a large mass fraction (up to $7\,\%$ ) of dissolved volatile components, mainly water, which has a strong impact on the viscosity of magma. The latter is often approximated by an Arrhenius-type dependence (Proussevitch & Sahagian Reference Proussevitch and Sahagian1998):

(2.1) \begin{equation} \eta (C) = \eta ^*\exp {\big\{\big(E^*_{\eta }(1-k_{\eta }\, C)/\mathcal{R} T\big)\big\}}, \end{equation}

where $\eta$ is the dynamic viscosity of the magmatic melt; $C$ is the mass fraction of gas dissolved in magma; $T$ is the magma temperature; $\mathcal{R}$ is the universal gas constant; and $\eta ^*$ , $E^*_{\eta }$ and $k_{\eta }$ are empirical constants. Note that, if necessary, a more precise empirical dependence $\eta (C)$ can be used in the model formulated below.

Since the viscosity of magma can increase by orders of magnitude during degassing, the transitional stage of bubble growth (in contrast with bubble growth in ordinary liquids) can last quite a long time, of the order of $10^2 {-} 10^3$ seconds. At this stage the effect of barodiffusion, which is of central interest in this paper, may appear and significantly affect the entire dynamic pattern of the process. Indeed, the formation of a diffusion boundary layer around the growing bubble entails the emergence of a large pressure gradient in the liquid and, as result, a barodiffusion flow comparable to diffusion. Let us consider this issue in more detail.

To describe the dynamics of bubble growth, we write down the Navier–Stokes equation taking into account the spherical symmetry of the problem:

(2.2) \begin{equation} \rho _l\left (\dfrac {\partial v}{\partial t} + v\dfrac {\partial v}{\partial r}\right ) = -\dfrac {\partial p_l}{\partial r} + \eta \left ( \dfrac {1}{r^2} \dfrac {\partial }{\partial r} \left ( r^2 \dfrac {\partial v}{\partial r}\right ) - \dfrac {2 v}{r^2} \right ) + 2 \dfrac {\partial \eta }{\partial r}\dfrac {\partial v}{\partial r}, \end{equation}

where $t$ is time; $r$ is the radial coordinate with origin in the centre of the bubble; $v$ is the radial velocity of the liquid; $p$ is pressure; $\rho$ is density; and $\eta$ is the dynamic viscosity of the liquid. Hereinafter, the subscripts ‘ $l$ ’ and ‘ $g$ ’ denote the liquid and gas phases, respectively. The liquid phase is assumed to be incompressible.

The dependence of the velocity of the liquid $v$ on the radial coordinate $r$ is found from the continuity equation as follows:

(2.3) \begin{equation} v = \dot {R} R^2/r^2, \end{equation}

with $R$ the bubble radius and $\dot {R} ={\rm d}R/{\rm d}t$ the bubble growth rate. For an incompressible liquid, the latter quantity coincides with the velocity of the liquid at the bubble boundary. Substituting (2.3) into (2.2), we get

(2.4) \begin{equation} \rho _l \left ( \dfrac {\ddot {R}R^2}{r^2} + \dfrac {2\dot {R}^2 R}{r^2} - \dfrac {2\dot {R}^2 R^4}{r^5} \right ) = -\dfrac {\partial p_l}{\partial r} - 4 \dfrac {\dot {R} R^2}{r^3}\dfrac {\partial \eta }{\partial r}. \end{equation}

To eliminate the unknown pressure gradient on the right-hand side of (2.4), we integrate this equation from $R$ to infinity:

(2.5) \begin{equation} \rho _l \left ( \ddot {R}R + \dfrac 32\dot {R}^2 \right ) = p_l(R) - p_l(\infty ) + 4 \eta (R)\dfrac {\dot {R}}{R} - 12 \dot {R}R^2 \int _R^{\infty } \dfrac {\eta (r)}{r^4} {\rm d}r. \end{equation}

Note that under changing external conditions, the pressure away from the bubble depends on time. In the idealised case of instantaneous decompression, it is equal to the final pressure.

The condition of equality of forces acting on the surface of the bubble from the gas and liquid sides is written as follows:

(2.6) \begin{align} p_g - \left \lbrace p_l(R) - 2 \eta (R) \left ( \dfrac {\partial v}{\partial r} \right )_{r=R} \right \rbrace = \dfrac {2\sigma }{R}, \end{align}

where $\sigma$ is the surface tension. Note that this condition is different from the classic condition for liquids with negligibly low viscosity by reason of the presence of the viscous term proportional to the product of velocity gradient and viscosity. At high value of viscosity, i.e. in the case considered in the present work, the role of the viscous terms in the stress tensor is crucially important for the bubble dynamics. Importantly, owing to the strong viscosity gradient, the bubble appears to be surrounded by a highly viscous ‘shell’. The effect of this shell should always be kept in mind as it defines the dynamic background of the system superposed by the effects of barodiffusion, which is of central interest in this paper.

Expressing the dependence of $p_l(R)$ from the last equation, we obtain

(2.7) \begin{equation} \rho _l \left ( R\ddot {R} + \dfrac {3}{2}\dot {R}^2 \right ) = p_g - p_l(\infty ) - \dfrac {2 \sigma }{R} - 4\hat {\eta } \dfrac {\dot {R}}{R}; \quad \hat {\eta } = 3 R^3 \int _R^{\infty } \dfrac {\eta (r)}{r^4} {\rm d}r. \end{equation}

Equation (2.7) is an analogue of the Rayleigh–Plesset equation with the only difference being that, instead of the usual viscosity, the integral viscosity $\hat {\eta }$ appears in its right-hand part. The non-uniformity of viscosity is due to formation of a diffusion boundary layer around the bubble. Since the viscosity of the magmatic melt strongly depends on the mass fraction of dissolved volatile components, the viscosity gradient around the bubble can take large values, which must be taken into account. It is easy to see that in the case of uniform viscosity, equation (2.7) reduces to the classical Rayleigh–Plesset equation.

The dynamics of the mass fraction of the gas dissolved in the melt $C(t,r)$ is described by the transfer equation, which, taking into account the effect of barodiffusion, is written as follows (unlike barodiffusion, thermodiffusion is not taken into account here, since the temperature gradient in the liquid turns out to be quite small; see below):

(2.8) \begin{equation} \dfrac {\partial C}{\partial t} + v \dfrac {\partial C}{\partial r} = \dfrac {1}{r^2} \, \dfrac {\partial }{\partial r} \left (r^2 \left \lbrace D \dfrac {\partial C}{\partial r} + \dfrac {D_p}{p} \dfrac {\partial p}{\partial r}\right \rbrace \right )\!, \end{equation}

where $D$ is the diffusion coefficient and $D_p$ is the barodiffusion coefficient. Note that the convective term appearing in (2.8) is due to the radial motion of the liquid.

In this paper we assume for simplicity that both the diffusion and barodiffusion coefficients are constant. The diffusion coefficient depends on the viscosity of the liquid (Landau & Lifshitz Reference Landau and Lifshitz1987), and, therefore, in the specific case of a magmatic melt, it depends on the mass fraction of dissolved gas. Physically, the barodiffusion coefficient and the diffusion coefficient are expected to have similar dependencies on the viscosity of the melt. Therefore, the barodiffusive ratio $k_p = D_p/D$ is expected to be nearly constant. Taking into account dependencies of the diffusion and barodiffusion coefficients on the mass fraction of dissolved gas may affect the growth rate of a bubble; however, it is unlikely to produce a strong impact on the relative contribution of barodiffusion to the gas flux into the bubble. The latter issue is addressed in our work.

At the initial moment of time, the melt is assumed to be saturated with gas to an equilibrium state: $(C)_{t=0} = C_i$ , where $C_i = C_s(p_i)$ . Here $C_s$ is the mass fraction of gas dissolved in the saturation state in magma, which, according to Henry’s law, is a function of pressure $p$ . For water dissolved in magma, this dependence has the following form (Stolper Reference Stolper1982): $C_s(p) = \sqrt {K_{{H}}\, p}$ , where $K_{ {H}}$ is Henry’s constant. Hereinafter, the subscripts ‘ $i$ ’ and ‘ $f$ ’ denote the initial and final states. The condition of local thermodynamic equilibrium is fulfilled at the bubble boundary: $(C)_{r=R} = C_s(p_g)$ . Away from the bubble, the mass fraction of the gas dissolved in magma is equal to the initial one: $(C)_{r\to \infty } = C_i$ . The gas in the bubble is considered ideal: $p_g =\rho _g\bar {\mathcal{R}} T$ , where $\bar {\mathcal{R}} =\mathcal{R}/\mu _g$ is the reduced gas constant with $\mu _g$ being the molar mass of the gas.

The formulated boundary-value problem is complemented by the equation of material balance:

(2.9) \begin{equation} \dfrac {{\rm d} m_g}{{\rm d} t} = 4\pi R^2 \rho _l \left ( D \dfrac {\partial C}{\partial r} + \dfrac {D_p}{p} \dfrac {\partial p}{\partial r} \right )_{r = R}\!, \end{equation}

where $m_g = (4\pi /3 ) \rho _g R^3$ is the mass of gas in the bubble.

Regarding the correctness of using a single-bubble growth model as a proxy for the process of ‘collective’ bubble growth during degassing, we believe that a single-bubble model remains reasonably accurate as long as the thickness of the diffusion boundary layer around the bubble is much less than the thickness of liquid domains separating individual bubbles. At large oversaturation this condition may be fulfilled even at relatively large volume fraction of gas since the thickness of the diffusion boundary layer is much less than the bubble size (Chernov et al. Reference Chernov, Kedrinsky and Pil’nik2014, Reference Chernov, Davydov and Pil’nik2023b ).

In conclusion, we emphasise that, in contrast to gases (Elperin, Kleeorin & Rogachevskii Reference Elperin, Kleeorin and Rogachevskii1997; Ke & He Reference Ke and He2012; Zhdanov Reference Zhdanov2017), there are currently no models that would reliably determine the value of the barodiffusion coefficient for liquids, especially for highly viscous ones, and the only way to find this coefficient remains experimental. It is only known that in liquids, as in gases, in the process of barodiffusion, heavy molecules tend to move to the area of increased pressure, and light molecules tend to move to the area of reduced pressure.

3. Analytical solution of the problem

The problem is solved using the following non-dimensional variables: $\tau = t/t_0$ ; $\bar {r} = r/R_0$ ; $\bar {R}= R/R_0$ ; $\bar {p} = (p - p_{\!f})/\Delta p_i$ ; and $\bar {C} = (C - C_f)/\Delta C_i$ , where $t_0$ is the characteristic time of the process; $R_0$ is the characteristic size; $\Delta p_i= p_i - p_{\!f}$ ; $\Delta C_i = C_i - C_f$ ; and $C_f = C_s(p_{\!f})$ . The characteristic size $R_0$ is chosen in such a way that the characteristic times of the diffusion $t_0^D = R_0^2/D$ and dynamic $t_0^{\eta } = \eta _i/\Delta p_i$ processes are equal, i.e. $t_0^D = t_0^{\eta } \equiv t_0$ . Then $R_0 = \sqrt {D\eta _i/\Delta p_i}$ , where $\eta _i = \eta (C_i)$ . We also introduce the characteristic speed of the process $v_0 = R_0/t_0 = \sqrt {D\Delta p_i/\eta _i}$ . Note that the above-defined choice of $R_0$ is the only reasonable one, as it reduces the number of free parameters of the problem and simplifies the parametric analysis of the entire process under consideration. For convenience, we introduce the variable $\chi = \bar {r}/\bar {R}$ : in terms of this variable the interphase boundary becomes stationary.

The modified Rayleigh equation (2.7) in non-dimensional variables is written as follows:

(3.1) \begin{equation} \boldsymbol {P\!r}_{\!\boldsymbol{D}}^{-1} \big ( \bar {R} \bar {R}'' + (3/2) \bar {R}'^2 \big ) = \bar {p}_g - \bar {p}_l(\infty ) - \bar {R}_{cr}/\bar {R} - 4 \hat {\bar {\eta }}\, \bar {R}'/\bar {R}, \end{equation}

where $\boldsymbol{P\!r}_{\!\boldsymbol{D}} = \nu _i/D$ is the Prandtl diffusion number; $\nu _i = \mu _i/\rho _l$ is kinematic viscosity of the melt; $\hat {\bar {\eta }} = 3 \int _1^{\infty } \bar {\eta }(\zeta )/\zeta ^4 {\rm d} \zeta$ is non-dimensional integral viscosity; $\bar {\eta } = \eta /\eta _i$ ; $\bar {R}_{cr} =R_{cr}/R_0$ ; $R_{cr} = 2\sigma /\Delta p_i$ is the critical radius of the bubble; $\bar {R}' = {\rm d}\bar {R}/{\rm d}\tau$ ; and $\bar {R}'' = {\rm d}^{2}\bar {R}/{{\rm d}\tau ^2}$ . Since the viscosity of magmatic melts is quite high (the characteristic values of the Prandtl diffusion number are of the order of $10^{12}$ and higher), the inertial terms in (2.1) can be neglected, which also applies to the equations for pressure and pressure gradient in liquid (see below). Also, if we consider the growth of bubbles whose sizes are much larger than the critical one, the Laplace pressure can be neglected in (2.1). However, if necessary, all these terms can be taken into account.

The pressure gradient in liquid in dimensionless variables in the absence of inertial terms will take the following form (see (2.4)):

(3.2) \begin{equation} \dfrac {\partial \bar {p}_l}{\partial \chi } = -4 \dfrac {\bar {R}'}{\bar {R}} \dfrac {1}{\chi ^3} \dfrac {\partial \bar {\eta }}{\partial \chi }. \end{equation}

The viscosity gradient of the liquid in (2.2) can be expressed in terms of the gradient of the mass fraction of the gas dissolved in the liquid:

(3.3) \begin{equation} \dfrac {\partial \bar {\eta }}{\partial \chi } = \left (\dfrac {{\rm d} \bar {\eta }}{{\rm d} \bar {C}}\right ) \boldsymbol{\cdot }\left (\dfrac {\partial \bar {C}}{\partial \chi }\right )\!. \end{equation}

Note that the pressure gradient in the liquid remains always positive; moreover, for some regimes in the absence of barodiffusion the pressure near the bubble may be negative. In general, owing to specific properties of magmatic melts, the pressure gradient in magma is due to the gradient of the mass fraction of the dissolved gas, which is precisely a distinctive feature of magma as compared with ordinary liquids. As for the flux due to barodiffusion, it is essential only in the liquid domain where the pressure gradient is high. The latter is directly related to the gradient of the mass fraction of dissolved gas. Thus, barodiffusive flux is essential only in the thin diffusive boundary layer formed around the bubble and only at the transient stage of the process (see below). The thickness of this layer is small compared with the radius of the bubble at the considered time scales of the process. Away from the bubble, both the pressure gradient and the flux due to barodiffusion are negligibly small.

If dependence (2.1) is used for the viscosity of the magmatic melt, which in dimensionless form has the form $\bar {\eta } = \exp \{ \varkappa _i\, (1 - \bar {C}) \}$ , then ${\rm d}{\bar {\eta }}/{{\rm d}\bar {C}}= - \varkappa _i\bar {\eta }$ , where $\varkappa _i = \varkappa _0 \Delta C_i$ ; and $\varkappa _0 = k_{\eta } E_{\eta }^* / (\mathcal{R} T)$ is a coefficient depending only on the properties of the magmatic melt. As can be seen from (2.3) and (2.2), the barodiffusion flow, like the diffusion flow, is directed towards the bubble, since a zone of low pressure is formed near the bubble. Molecules of water, as is well known, are lighter than molecules of silicon dioxide, the key constituent of the ‘dry’ magmatic melt.

The pressure distribution in the liquid around the bubble can be found by integrating (2.2):

(3.4) \begin{equation} \bar {p}_l(\chi ) = \bar {p}_l(\infty ) - 4 \dfrac {\bar {R}'}{\bar {R}} \left (\dfrac {\bar {\eta }(\chi )}{\chi ^3} - 3\int _{\chi }^{\infty } \dfrac {\bar {\eta }(\zeta )}{\zeta ^4} {\rm d}{\zeta } \right )\!. \end{equation}

For pressure and pressure gradient at the interface, we have the following expressions:

(3.5) \begin{align}& \big(\partial \bar {p}_l / \partial \chi \big)_{\varGamma } = - 4 \big(\bar {R}'/\bar {R}\big) \big( {\rm d}\bar {\eta }/ {\rm d} \bar {C}\big)_{\varGamma } \big(\partial \bar {C}/\partial \chi \big)_{\varGamma }, \end{align}

(3.6) \begin{align}&\qquad\,\, \bar {p}_{l \varGamma } = \bar {p}_l(\infty ) - 4 \big(\bar {R}'/\bar {R}\big) \big(\bar {\eta }_{\varGamma } - \hat {\bar {\eta }}\big)\!. \end{align}

Hereinafter, the subscript ‘ $\varGamma$ ’ denotes the value of a quantity at the interphase boundary, i.e. at $\chi = 1$ .

To find a solution to the diffusion problem (2.8), (2.9), we use an approach similar to that of Chernov & Pil’nik (Reference Chernov and Pil’nik2015, Reference Chernov and Pil’nik2018) and Chernov et al. (Reference Chernov, Pil’nik, Vladyko and Lezhnin2020) that proved to be efficient for a description of the growth of new phase germs in various physical systems in a metastable state. Using variables ( $\tau$ , $\chi$ ), the problem to be solved takes the following form:

(3.7) \begin{align}& \dfrac {\partial \bar {C}}{\partial \tau } = \dfrac {\bar {R}'}{\bar {R}} \left ( -\dfrac {1}{\chi ^2} + \chi \right ) \dfrac {\partial \bar {C}}{\partial \chi } + \dfrac {1}{\bar {R}^2} \dfrac {1}{\chi ^2} \dfrac {\partial }{\partial \chi } \left ( \chi ^2 \alpha \, \dfrac {\partial \bar {C}}{\partial \chi } \right )\!; \end{align}

(3.8) \begin{align}&\quad (\bar {C})_{\tau = 0} = 1; \quad (\bar {C})_{\chi = 1} = \bar {C}_{\varGamma }; \quad (\bar {C})_{\chi \to \infty } = 1, \\[12pt] \nonumber \end{align}

where $\alpha (\tau ,\chi ) = 1-4 k_p \chi ^{-3} (\bar {R}'/\bar {R} ) (\bar {p}_{l}(\tau ,\chi ) + \psi ^{-1} )^{-1} {\rm d} \bar {\eta }(\tau ,\chi )/{\rm d}C$ ; ${\rm d}\bar {\eta }/{\rm d}C = (\Delta C_i )^{-1} {\rm d}\bar {\eta }/{\rm d}\bar {C}$ ; $\bar {C}_{\varGamma }(\tau ) = \bar {C}_s (\bar {p}_g(\tau ) ) = {(\sqrt {\psi \bar {p}_g (\tau )+ 1} -1)}/{(\sqrt {\psi + 1} - 1)}$ ; $\psi = \Delta p_i/ p_{\!f}$ is relative decompression value (at $\psi \gg 1$ : $\bar {C}_{\varGamma }(\tau ) \approx \sqrt {\bar {p}_g (\tau )}$ ); and $k_p = D_p/D$ is the barodiffusion ratio. Note that if barodiffusion is neglected, then $\alpha = 1$ .

As is shown in Chernov et al. (Reference Chernov, Pil’nik, Davydov, Ermanyuk and Pakhomov2018, Reference Chernov, Pil’nik, Vladyko and Lezhnin2020), in terms of $\tau$ , $\chi =\bar {r}/\bar {R}$ a quasi-stationary state is quickly reached in the system. Therefore, when integrating (3.7) one can neglect the time derivative $\partial \bar {C}/ \partial \tau$ , which appears on the left-hand side of (3.7). Within the framework of this assumption, an analytical solution can be found that agrees well with the numerical one in a wide range of governing parameters of the process, and describes all its stages: from inertial to diffusion. This solution has the form

(3.9) \begin{equation} \bar {C}(\tau , \chi ) = 1 - \left \lbrace 1 - \bar {C}_{\varGamma }(\tau ) \right \rbrace \, \dfrac {\boldsymbol{I}\,(\tau , \chi , \beta )}{\boldsymbol{I}\,(\tau ,1, \beta )}, \end{equation}

where $\boldsymbol{I}\,(\tau ,\chi , \beta ) = \int _{\chi }^{\infty } ({1}/{(\alpha \zeta ^2)}) \exp \lbrace - ({\beta }/{\alpha }) (({1}/{\zeta }) + ({\zeta ^2}/{2))}\rbrace {\rm d} {\zeta }$ is the integral function of coordinate, time and time-dependent functions $\alpha (\tau )$ and $\beta (\tau ) = \bar {R} \bar {R}'$ , which include the sought-for dependence $\bar {R}(\tau )$ . Taking into account (3.9), the latter can be found from the equation of material balance (2.9) and expressed in terms of the following integro-differential equation:

(3.10) \begin{equation} \boldsymbol{I}^*\,(\tau ) \, \dfrac {{\rm d}}{{\rm d} \tau } \big \lbrace \big(1 + \psi \bar {p}_g \big) \bar {R}^3 \big\rbrace = 3\, \varepsilon \bar {R} \left (1 - \bar {C}_{\varGamma }\right )\!, \end{equation}

where $\boldsymbol{I}^*\,(\tau , \beta ) = \int _1^{\infty } ({1}/{(\alpha \zeta ^2)}) \exp \lbrace - ({\beta }/{\alpha }) (({1}/{\zeta }) + ({\zeta ^2}/{2})) + (3/2)({\beta }/{(\alpha _{\varGamma })}) \rbrace {\rm d} {\zeta }$ ; $\varepsilon ={(\rho _l \Delta C_i)}/{(\rho _{g}(p_{\!f}))}$ is the dimensionless phase transformation criterion characterising the initial supersaturation of the medium caused by decompression, which is essentially an analogue of the Jacob number appearing in the problem of the growth of a vapour bubble in an overheated liquid (Chernov et al. Reference Chernov, Pil’nik, Vladyko and Lezhnin2020); and

(3.11) \begin{equation} \alpha _{\varGamma } (\tau ) = 1+j_p/j = 1 - 4 k_p \big(\bar {R}'/\bar {R}\big) \big(\bar {p}_{l \varGamma } + \psi ^{-1}\big)^{-1} \big({\rm d}\bar {\eta }/{{\rm d}C}\big)_{\varGamma }\!. \end{equation}

Note that in the latter formula, the second term on the right-hand side is the ratio of the density of the gas flow $j_p$ on the interfacial surface due to barodiffusion to the density of the flow $j$ due to molecular diffusion. For large supersaturation ( $\varepsilon \gg 1$ ), which corresponds to large values of $\beta$ typical for real volcanic processes (Chernov et al. Reference Chernov, Pil’nik, Davydov, Ermanyuk and Pakhomov2018), the integral expression $\boldsymbol{I}^*$ has a simple asymptotic form: $\boldsymbol { I}^*\approx \sqrt {\pi /(6\alpha _{\varGamma } \beta )} - 4/(9\beta )$ .

Regarding the integral viscosity, let us note that since the integral function in the expression for $\hat {\bar {\eta }}$ (see (2.7)) at large $\varkappa _i$ decreases very quickly with the coordinate $\chi$ , finding $\hat {\bar {\eta }}$ can be greatly simplified. Representing the viscosity of the melt as an exponential function of the mass fraction of the dissolved gas, decomposing the exponent into the Taylor series in the vicinity of $\chi = 1$ so that $\bar {\eta } \approx \bar {\eta }_{\varGamma }\, \exp \lbrace \displaystyle ({({\rm d}\ln \bar {\eta })}/{({\rm d}\bar {C})})_{\varGamma } ({(\partial \bar {C})}/ {(\partial \chi)} )_{\varGamma } (\chi -1 ) \rbrace$ , integrating and taking into account (3.9), we get

(3.12) \begin{align} \hat {\bar {\eta }} \approx \bar {\eta }_{\varGamma } \dfrac {3\, \alpha _{\varGamma }\, \boldsymbol { I}^*}{\left (-\dfrac {{\rm d}\ln \bar {\eta }}{{\rm d}\bar {C}}\right )_{\varGamma } \left (1 - \bar {C}_{\varGamma }\right ) }. \end{align}

If (2.1) is used for the viscosity of the magmatic melt, then ${({\rm d}\ln \bar {\eta })}/{({\rm d}\bar {C})} = -\varkappa _i$ .

Note that the obtained solution can also be used under non-stationary external conditions, since the pressure $\bar {p_l}(\infty )$ on the right-hand side of (2.1) is generally a function of time (while in the particular case of instantaneous decompression $\bar {p}_l(\infty ) = 0$ ).

Let us consider some limiting cases. At short times ( $\tau \ll 1$ ), when the pressure in the bubble is almost equal to the initial one, in the case of instantaneous decompression, (2.7) can be integrated. As a result, we get

(3.13) \begin{equation} \bar {R}(\tau ) = \bar {R}_{cr}(1 + \delta \, e^{\tau /4}), \end{equation}

where $\delta =R_i/R_{cr} - 1$ in which $R_i$ is the initial (supercritical) radius, which, generally speaking, can be chosen with some arbitrariness provided that it somewhat exceeds the critical value.

At long times ( $\tau \gg 1$ ), the gas pressure in the bubble asymptotically tends to the pressure of the ambient liquid (i.e. away from the bubble): $\bar {p}_g\to 0$ , and function $\beta$ tends to a constant value $\beta _0$ , depending only on the degree of the metastability $\varepsilon$ , while the concentration field of the gas dissolved in the liquid tends to its stationary state, which indicates the self-similarity of this stage of the process. Here the self-similar variable is $\chi$ . This is the diffusion stage where the effects of inertia and barodiffusion are no longer significant (see expression for $\alpha _{\varGamma }$ (3.11)), and bubble growth occurs solely as a result of molecular diffusion. At this stage the distribution of the mass fraction of the dissolved gas and the law of bubble growth can be found explicitly:

(3.14) \begin{equation} \bar {C}(\chi ) = 1 - \boldsymbol{I}\,(\chi , \beta _0)/\boldsymbol{I}\,(1, \beta _0); \quad \bar {R}(\tau ) = \sqrt {\bar {R}_*^2 + 2\beta _0 \left (\tau - \tau _*\right )}, \end{equation}

where $\tau _*$ and $\bar {R_*} = \bar {R}(\tau _*)$ are the characteristic time of establishment of the diffusion stage and the radius of the bubble at this moment of time, respectively. The coefficient $\beta _0$ is found from the integral equation: $\beta _0\, \boldsymbol{I}^*\,(\beta _0) = \varepsilon$ . For large supersaturation ( $\varepsilon \gg 1$ ), there is a good asymptotic decomposition for $\beta _0$ , namely $\beta _0\approx (6/\pi )\, (\varepsilon +4/9)^2$ . It should be noted that (3.14) can be obtained by solving a purely diffusion problem, without taking into account inertial effects and effects associated with barodiffusion (Scriven Reference Scriven1959). Once again, we draw attention to the fact that the diffusion stage of growth for bubbles in magmatic melts is not necessarily achievable within a reasonable span of time. This is defined by conditions under which a particular natural process takes place.

The obtained asymptotic solutions (3.13) and (3.14) allow us to estimate the contribution of barodiffusion to the short-time and long-time dynamics of bubble growth:

(3.15) \begin{align}&\quad\,\,\, j_p/j \approx k_p \varkappa _0 \psi \quad {\text{at}} \quad \tau \lesssim 1; \end{align}

(3.16) \begin{align}& j_p/j\approx 2 k_p \varkappa _0 \psi \bar {\eta }_f \tau ^{-1} \quad {\text{at}} \quad \tau \gg 1. \end{align}

As can be seen from the above formulas, the effect of barodiffusion is significant only at the initial and transitional stages. The impact of this effect decreases with time, and can be safely ignored at $\tau \gg 2 k_p \varkappa _0 \psi \bar {\eta }_f$ .

4. Thermal problem

Let us evaluate the possible influence of thermal effects on the dynamics of magma degassing. At first glance, this effect is deemed to be significant, since a lot of energy is spent on the desorption of gas dissolved in the magmatic melt that may affect, in particular, the conditions at the interphase boundary. When the entire excess mass of gas formed as a result of decompression is released from the magma, the melt temperature decreases by $T_i\, {\boldsymbol {Ph}}$ , where ${\boldsymbol {Ph}} = r_d \Delta C_i/(c_l T_i)$ is the criterion of phase transformation (analogous to the Stefan number), with $r_d$ the specific heat of desorption and $c$ the specific heat capacity. For typical values of the mass fraction of dissolved gas ( ${\sim}5\,\%$ ) an overall temperature decrease of the order tens of kelvins may be expected. To describe in detail the dynamics of the magma cooling in the vicinity of a growing single bubble, we need to solve the following thermal problem.

The temperature field in a magmatic melt is described by the thermal conductivity equation:

(4.1) \begin{align} \rho _l c_l \left ( \dfrac {\partial T}{\partial t} + v \dfrac {\partial T}{\partial r} \right ) = \dfrac {1}{r^2}\dfrac {\partial }{\partial r} \left (\lambda _l r^2 \dfrac {\partial T}{\partial r} \right )\!, \end{align}

where $\lambda$ is the coefficient of thermal conductivity. We assume that at the initial moment of time, the temperature of the magma is homogeneous $ (T )_{t=0} = T_i$ and exceeds the liquidus temperature so that the cooling due to gas desorption does not lead to partial crystallisation. The thermal balance condition is fulfilled at the bubble boundary: $-\lambda _l ({\partial T}/{\partial r} )_{r=R}= - r_d \rho _l D ({\partial C}{\partial r} )_{r=R}$ . At infinity, the temperature field remains undisturbed, i.e. $ (T )_{r\to \infty }\to T_i$ . Being formulated in non-dimensional variables, the thermal problem takes the following form:

(4.2) \begin{equation} \dfrac {\partial \varTheta }{\partial \tau } = \dfrac {{\bar {R}'}}{\bar {R}} \left ( -\dfrac {1}{\chi ^2} + \chi \right ) \dfrac {\partial \varTheta }{\partial \chi } + {\boldsymbol {Le}}^{-1} \dfrac {1}{\bar {R}^2}\dfrac {1}{\chi ^2}\dfrac {\partial }{\partial \chi } \left (\chi ^2\dfrac {\partial \varTheta }{\partial \chi }\right )\!; \nonumber\\[3pt] \left (\varTheta \right )_{\tau = 0} = 1; \nonumber\\[3pt] \left (\varTheta \right )_{\chi = 1} = \varTheta _{\varGamma }; \quad -\left ({\partial \varTheta }/{\partial \chi }\right )_{\chi = 1} = -{\boldsymbol {Ph}} \, {\boldsymbol {Le}} \, \left ({\partial \bar {C}}/{\partial \chi }\right )_{\chi = 1}; \quad \left (\varTheta \right )_{\chi \to \infty } = 1. \end{equation}

Here $\varTheta = T/T_i$ is dimensionless temperature; ${\boldsymbol {Le}}= D/a_l$ is the Lewis number; and $a = \lambda / (\rho c )$ is the thermal conductivity coefficient.

By analogy with the diffusion problem, the solution of the thermal problem can be obtained in the following form:

(4.3) \begin{equation} \varTheta (\tau , \chi , \beta ) = 1 - \left \lbrace 1 - \varTheta _{\varGamma }(\tau ) \right \rbrace \, \dfrac {\boldsymbol{I}\,(\tau , \chi , \beta {\boldsymbol {Le}})}{\boldsymbol{I}\,(\tau , 1, \beta {\boldsymbol {Le}})}, \end{equation}

where $\varTheta _{\varGamma }$ is the temperature at the gas–liquid interface. It should be noted here that, in contrast with the diffusion problem, in the solution of the thermal problem (4.3) in the integral function $\boldsymbol{I}$ the argument $\beta$ is replaced by $\beta {\boldsymbol {Le}}$ . Substituting the temperature profile (4.3) into the heat balance equation, we get

(4.4) \begin{equation} \varTheta _{\varGamma }(\tau ) = 1 - {\boldsymbol {Ph}} \, {\boldsymbol {Le}} \, \dfrac {\boldsymbol{I}^*\,(\tau , \beta {\boldsymbol {Le}})}{\boldsymbol{I}^*\,(\tau , \beta )}. \end{equation}

Formula (4.4) determines the temperature of the melt at the gas–liquid interface where the effect of cooling is most pronounced. As can be seen from (4.3) and (4.4), at the diffusion stage of the process the temperature profile around the bubble becomes stationary and the resulting solution is self-similar.

5. Numerical results and discussion

This section presents the results of numerical simulation of the process of gas bubble growth in a rhyolitic magmatic melt, which is subject to rapid decompression at the initial moment of time. This particular type of magma is chosen because the viscous effects, and hence the effects associated with barodiffusion, are expected to be quite pronounced. The calculations were carried out based on the analytical solution obtained above, and also via numerical simulations. Initially, the melt is assumed to be gas-saturated. It should be noted that our scope is to clearly demonstrate the effect of barodiffusion under specific conditions rather than to describe all possible scenarios of bubble growth covering a wide range of governing parameters and initial conditions. The calculation results are presented in dimensionless form, which greatly simplifies interpretation since all spatial and temporal scales of the process are directly related to the corresponding reference values, and the overall number of free parameters is decreased.

All physicochemical properties of magma used in calculations (table 1) are taken from Proussevitch & Sahagian (Reference Proussevitch and Sahagian1998), which presents a fairly complete summary of the essential data available in the literature. Formula (2.1) is used to describe the dependence of the magma viscosity on the mass fraction of the dissolved gas (water). The choice of specific value of the barodiffusion coefficient remains problematic due to the lack of any reliable information in the literature. In this regard, the barodiffusion coefficient in this work is set as a free parameter. By varying this parameter, we can say under which conditions barodiffusion has (or has not) a significant effect on the overall dynamics of the bubble growth. Naturally, the choice of this coefficient should be varied within physically reasonable limits. We expect that its value should not differ much from the estimated value for gases. Accordingly, the value of the barodiffusion ratio $k_p$ in the calculations was chosen to be $10^{-3}$ and $2{\times }10^{-3}$ . Table 2 presents the values of initial parameters, characteristic values and dimensionless criteria at which the calculations were performed.

Table 1.

Some properties of rhyolite magmatic melt.

Table 2.

Initial parameters, characteristic values and dimensionless criteria of the process used in numerical simulations.

In practice, the key parameters of the problem vary within a wide range. The initial pressure in the magma chamber may reach 2000 atm, and then reduce to normal atmospheric pressure during decompression. The value of initial decompression defines the initial oversaturation of the medium. The decompression may be either very fast (e.g. due to rarefaction wave propagating along the magma conduit at the initial stage of an explosive-type eruption) or quite slow (e.g. due to ascent of highly viscous magma in the case of an extrusive eruption). It is impossible to model all possible scenarios within a single paper. Therefore, we restricted the scope of the present research to demonstration of existence of regimes where barodiffusion plays essential role in the dynamics of bubble growth. The choice of initial conditions in our calculations (table 2) falls into the range of immediate interest for realistic volcano eruptions. This choice is motivated by our intention to consider two cases characterised by small (Set 1) and large (Set 2) values of non-dimensional criterion of phase transition $\varepsilon$ introduced in our paper. Note that with our choice this quantity varies by roughly an order of magnitude to ensure consistency of conclusions.

Figure 1.

Dependencies of pressure in liquid $\bar {p}_l$ (a) and mass fraction of dissolved gas $\bar {C}$ (b) on radial coordinate $\chi$ at time $\tau = 30$ without taking into account barodiffusion, $k_p = 0$ (1), and taking into account barodiffusion, $k_p = 10^{-3}$ (2), Set 2.

Figure 1 shows the liquid pressure $\bar {p}_l$ and the mass fraction of dissolved gas $\bar {C}$ as functions of the radial coordinate $\chi$ (Set 2) at the moment in time ( $\tau = 30$ ) when the barodiffusion flux reaches its maximum value. At this time, the difference between the solutions obtained with and without taking into account barodiffusion is maximum. As can be seen in figure 1, barodiffusion leads to extremely sharp variation of $\bar {p}_l(\chi )$ and $\bar {C}(\chi )$ in the immediate vicinity of the bubble surface. This leads to a consistent increase in the bubble growth rate as compared with the case where barodiffusion is not taken into account. The difference between the dependencies $\bar {p}_l(\chi )$ and $\bar {C}(\chi )$ corresponding to $k_p = 10^{-3}$ and $k_p = 2{\times }10^{-3}$ is almost imperceptible; for this reason the latter case is not shown in figure 1.

The effect of pressure decrease in the liquid near the bubble can be tentatively explained as follows. In the presence of a viscosity gradient, the bubble growth rate is lower than it would be in the case of uniform viscosity under the same conditions (as discussed in our previous paper (Chernov et al. Reference Chernov, Pil’nik, Davydov, Ermanyuk and Pakhomov2018)). Since the bubble in the presence of the viscosity gradient does not grow to the size it would reach in the case of uniform viscosity, tensile stresses arise in the liquid, being manifested by a pressure decrease in the diffusion boundary layer around the bubble, where the viscosity gradient is significant. Barodiffusion somewhat smooths out this effect.

In figure 2(a,b), the time dependencies of the bubble radius and gas pressure in the bubble are presented, calculated within the framework of (3.10), both with and without barodiffusion taken into account. The major part of the calculated results are presented in logarithmic scale, since the characteristic spatial and temporal scales of the various stages of the process differ by orders of magnitude. The same dependencies are shown in figure 2(c,d) for a narrow time interval as illustration of the dynamics observed at the transient stage. One can clearly see all three stages of the process under consideration: (i) inertial (exponential), (ii) transient and (iii) diffusive. Naturally, the transition stage is of central physical interest, representing also a significant difficulty in modelling. At this stage the effects associated with inertia and barodiffusion manifest themselves. Note that the duration of this stage can be long and, depending on the type of magma and the conditions of the process, can reach hundreds or even thousands of seconds. This is due to the high viscosity of magmas. For ordinary liquids, the transition stage is usually short, and, therefore, it is not taken into account.

Figure 2.

Dependencies of the bubble radius $\bar {R}$ (solid lines 1–3) and gas pressure in the bubble $\bar {p}_g$ (solid lines 4–6) on time $\tau$ : excluding barodiffusion, $k_p = 0$ (1, 4) and taking into account barodiffusion, $k_p = 10^{-3}$ (2, 5) and $k_p = 2\times 10^{-3}$ (3, 6). The dashed lines 7 illustrate the asymptotic solutions of (3.13) and (3.14). (a,b) Set 1; (c,d) Set 2.

In figure 2 it can be seen that without barodiffusion the gas pressure in the bubble drops monotonically. By contrast, if barodiffusion is taken into account, the gas pressure time history exhibits non-monotonic behaviour, with development of a secondary pressure peak in a certain time interval. This effect becomes more pronounced for higher values of $\varepsilon$ and $k_p$ . At this stage of the bubble growth the pressure gradient in the liquid reaches very high values, intensifying the barodiffusion flow. Under certain conditions, this flow can surpass the diffusion flow, which is indirectly manifested in the observed pressure dynamics. At the subsequent stage, the overall dynamics is controlled by diffusion and the effect of barodiffusion ceases to manifest itself. Therefore, regardless of whether barodiffusion is taken into account or not, for sufficiently large values of time the dependence of the bubble radius on time asymptotically tends to (3.14).

It is worth emphasising that, in contrast to ordinary liquids, the bubble growth in a magmatic melt is characterised by formation of a low-pressure region near the bubble surface (see (3.4)). The higher the supersaturation (in other words, the greater the value of $\varepsilon$ ), the lower is the pressure in the vicinity of the bubble boundary. In the absence of barodiffusion, the pressure near the bubble might become negative for certain regimes. In other words, tensile stresses might arise in a narrow liquid layer around the bubble, the thickness of which is comparable to the thickness of the diffusion boundary layer. That would lead to intense cavitation of the liquid, with formation of a cloud of small bubbles around the growing central bubble. Such a scenario is physically unlikely, although it follows from the obtained solution. An alternative scenario is expected if we take into account the effects associated with barodiffusion. Indeed, the higher the values of $\varepsilon$ , the greater is the pressure gradient (and the flow due to barodiffusion) in the liquid during bubble growth. When the pressure near the bubble reaches very low values, the barodiffusion flow exceeds the diffusion flow, thereby hindering the development of a negative pressure region in the liquid and suppressing cavitation. Note, however, that in the present work we use a linear dependence of the density of the barodiffusion flow on the pressure gradient. Whether or not this dependence is correct for magmatic melts under extreme conditions remains an open issue, as well as the question regarding the reference value of the barodiffusion coefficient.

Figure 3.

Dependence of the coefficient $\beta$ on time $\tau$ excluding barodiffusion, $k_p = 0$ (1) and taking into account barodiffusion, $k_p = 10^{-3}$ (2), $ k_p = 2\times 10^{-3}$ (3); the dashed line (4) indicates the asymptotic $\beta =\beta _0$ . (a) Set 1; (b) Set 2.

The effect of barodiffusion on the overall dynamics of the bubble growth can also be clearly illustrated by the dependence of function $\beta$ on time $\tau$ as shown in figure 3. In fact, function $\beta (\tau )$ determines the growth rate of the bubble. It has a clearly defined long-time asymptotic value $\beta _0 = \textrm {const.}$ ; the approach of $\beta (\tau )$ to $\beta _0$ identifies the transition to the diffusion stage of the bubble growth. As can be seen in figure 3, the effect of barodiffusion leads to faster growth of $\beta$ as compared with the case without barodiffusion. This difference manifests itself at the initial and transient stages, being more pronounced at higher values of $\varepsilon$ and $k_p$ . At large $\tau$ , when the barodiffusion flow becomes small, function $\beta (\tau )$ asymptotically tends to $\beta _0$ . However, it is important that in the case with barodiffusion $\beta (\tau )$ may significantly overshoot $\beta _0$ at the transient stage (see the data for Set 2 in figure 3).

The time history of function $j_p/j$ , the density of the barodiffusion flow $j_p(\tau )$ at the gas–liquid interface normalised by the density of the diffusion flow $j(\tau )$ , is shown in figure 4. Initially the ratio $j_p/j$ is nearly constant. It then increases with time at the transient stage, reaches a maximum value and asymptotically drops to zero value at the diffusion stage. One can clearly observe the existence of time intervals where both flows are comparable or where the flow due to barodiffusion significantly exceeds the diffusion flow. The existence of such intervals and relevant degassing regimes depend on particular conditions: the relevant value of dimensionless criterion $\varepsilon$ is the key factor here.

Figure 4.

Time history of the ratio $j_p/j$ of the density of the barodiffusion flow $j_p(\tau )$ at the gas–liquid interface to the density of the diffusion flow $j(\tau )$ : $k_p = 10^{-3}$ (1); $ k_p = 2\times 10^{-3}$ (2); the dashed lines illustrate the asymptotic solutions (3.15) and (3.16). (a) Set 1; (b) Set 2.

Note that the obtained asymptotic solutions (3.15) and (3.16) describe the initial and final stages of the process quite well and can serve as a simple estimate of the magnitude of the initial decompression and the characteristic time at which the influence of barodiffusion is significant. So, if $\psi \sim (k_p\varkappa _0 )^{-1}$ , then the flow due to barodiffusion at $\tau \lesssim 2 k_p\varkappa _0 \psi \bar {\eta }_f$ must be taken into account (see (3.15) and (3.16)). Although this is an upper estimate, it remains reasonably accurate when the value of $\varepsilon$ is not too large and the peak value of $j_p/j$ may significantly exceed the initial value.

Now let us discuss the thermal aspects of the problem. As mentioned above, during the growth of the bubble, a temperature boundary layer is formed near the bubble, where the melt is cooled due to absorption of latent heat associated with the phase transition during desorption. For the case of a significant decrease of the melt temperature, this would affect the diffusion problem, in particular, the boundary conditions at the gas–liquid interface. This assumption is intuitively justified owing to the large mass fraction of gas dissolved in magmas and significant energy required for its release. However, the calculations based on (4.3) and (4.4) show that the melt around the bubble cools only slightly: the temperature of the melt near the gas–liquid interface is lower than the temperature at infinity by a little more than 1 K. The reason for this lies in the fact that the thickness of the temperature boundary layer is significantly greater than the thickness of the diffusion layer (the Lewis number for magmatic melts is small: ${\boldsymbol {Le}}\sim 10^{-4}$ ). Therefore, with a good accuracy, the process can be considered isothermal. However, this is valid only for the problem of a single-bubble growth. In the case of growth of a bubble ensemble, with overlap of the temperature boundary layers of neighbouring bubbles at the late stage of the process, the melt will cool significantly.

Verification of the obtained analytical solution was carried out by numerical simulation. Computations were carried out using HPC computing resources of the Lavrentyev Institute of Hydrodynamics SB RAS. Equation (2.1) was integrated using the fourth-order Runge–Kutta–Merson scheme with automatic time step correction, and (3.7) and (3.10) (taking into account the initial and boundary conditions (3.8)) were integrated using an implicit scheme of the first order in time and the second order in spatial variable. The calculated results show good agreement with analytical solutions (3.9) and (3.10) and limiting cases (3.13) and (3.14). The difference between numerical and analytical solutions was less than $2\,\%$ .

Finally, we stress that whatever is the value of the barodiffusion coefficient, there exist relevant regimes and transient time intervals where the effects due to barodiffusion play a significant or even leading role. This is important for adequate description of processes physically related to gas bubble growth in highly viscous liquids (e.g. magma degassing during its ascent through a volcano channel). Indeed, at the transient stage of such processes, i.e. at a time scale much less than or comparable with the time scale of reaching the diffusion stage, the effect of barodiffusion fully manifests itself. Neglecting barodiffusion at this stage can lead to a significant error in describing the evolution of a single bubble in particular and the entire degassing process as a whole.