Two-phase magma flow with phase exchange. Part II. 1.5D numerical simulations of a volcanic conduit

Received: 10 October 2023 Revised: 14 May 2024 Accepted: 2 July 2024 DOI: 10.1111/sapm.12747 ORIGINAL ARTICLE Two-phase magma flow with phase exchange: Part II. 1.5D numerical simulations of a volcanic conduit Alain Burgisser1Marielle Collombet1Gladys Narbona-Reina2 Didier Bresch3 1Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTerre, Grenoble, France 2Dpto. Matemática Aplicada I E.T.S Arquitectura, Universidad de Sevilla, Sevilla, Spain 3LAMA UMR5127 CNRS, Université Savoie Mont Blanc, Le Bourget du lac, France Correspondence Alain Burgisser, Univ. Savoie Mont Blanc, CNRS,ISTerre,LeBourgetduLac, France. Email: [email protected] Funding information Agence Nationale pour la Recherche, Grant/Award Numbers: ANR-19-CE31-0007, ANR-18-CE40-0027, ANR-23EXMA-004, ANR-23-CE40-0014-01; MICIU/AEI/10.13039/501100011033 and European Union NextGenerationEU/ PRTR, Grant/Award Numbers: REC-B-22376-1, REC-B-22376-2; MCIN/AEI/10.13039/501100011033 and ERDF, Grant/Award Number: PID2022-137637NB-C22; University Savoie Mont Blanc; CNRS Abstract In a review paper in this same volume, we present the state of the art on modeling of compressible viscous flows ranging from single-phase to two-phase systems. It focuses on mathematical properties related to weak stability because they are important for numerical resolution and on the homogenization process that leads from a microscopic description of two separate phases to an averaged two-phase model. This review serves as the foundation for Parts I and II, which present averaged two-phase models with phase exchange applicable to magma flow during volcanic eruptions. Part I establishes a two-phase transient conduit flow model ensuring: (1) mass and volatile species conservation, (2) disequilibrium degassing considering both viscous relaxation and volatile diffusion, and (3) dissipation of total energy. The relaxation limit of this model is then used to obtain a drift-flux system amenable to simplification. Here, in This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made. © 2024 The Author(s). Studies in Applied Mathematics published by Wiley Periodicals LLC. Stud Appl Math. 2024;153:e12747. wileyonlinelibrary.com/journal/sapm 1of54 https://doi.org/10.1111/sapm.12747 2of54 BURGISSER et al. Part II, we summarize this model and propose a 1.5D simplification of it that alleviates three issues causing difficulties in its numerical implementation. We compare our model outputs to those of another steady-state, equilibrium degassing, isothermal model under conditions typical of an effusive eruption at an andesitic volcano. Perfect equilibrium degassing is unreachable with a realistic water diffusion coefficient because conduit extremities always contain melt supersaturated with water. Such supersaturation has minor consequences on mass discharge rate. In contrast, releasing the isothermal assumption reduces significantly mass discharge rate by cooling due to gas expansion, which in turn increases liquid viscosity. We propose a simplified system using Darcy’s law and omitting several processes such as shear heating and liquid inertia. This minimal system is not dissipative but approximates the steadystate mass discharge rate of the full system within 10%. A regime diagram valid under a limited set of conditions indicates when this minimal system captures the ascent dynamics of effusive eruptions. Interestingly, the two novel aspects of the full model, diffusive degassing and heat balance, cannot be neglected. In some cases with high diffusion coefficients, a shallow region where porosity and velocities tend toward zero develops initially, possibly blocking an eventual steady state. This local porosity loss also occurs when a steady-state solution is subjected to a change in shallow permeability. The resulting shallow porosity loss features many characteristics of a plug developing prior to a Vulcanian eruption. KEYWORDS numerical simulations, two-phase models, volcanic eruptions 1 INTRODUCTION The diversity of volcanic eruptive regimes, from gas-rich, particle-laden jets to slow effusion of gas-poor lavas, is a direct consequence of the multiphase nature of magmas, which are composed of silicate liquid, crystals, and a compressible fluid that can be either a gas or a supercritical fluid. 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 3of54 This diversity, which has not yet been rationalized in a single physical framework, led early modeling attempts to focus on end-member cases. Wilson,1for instance, described explosive, gas-rich eruptions with the energy equation of an expanding gas to relate gas pressure to ejection velocity. The propelling volcanic gas, however, is made of volatile species that, at depth, are either dissolved into the silicate liquid or present as individual gas bubbles. The eruptive regime depends greatly on how these volatiles come out of solution during the ascent toward the surface to form a continuous gas phase. The changes in pressure and temperature due to the ascent deeply affect the relative proportions of silicate liquid, crystals, and gas, and a thorough review of the associated physical processes can be found in Gonnermann and Manga.2 Quickly,3gas-based models were complemented by models of volcanic conduit flow based on the mass and momentum conservations of a single-phase magma where gas bubbles, silicate liquid, and crystals are represented by a Navier–Stokes equation with Boussinesq approximation along a one-dimensional cylindrical conduit assuming a parabolic velocity distribution across the conduit (often referred to as the 1.5D approximation) (e.g., Refs. 4–8). The single-phase assumption was motivated by the high viscosity of most silicate liquids, which restricts the motion of crystals and bubbles relative to the liquid. While this relative immobility holds for individual gas bubbles, it does not accurately represent the observed tendencies of bubbles to coalesce and form an interconnected network permeable to the gas they contain. Considering magma as a singlephase fluid also leads to an unavoidable separation of the Strombolian and Vulcanian regimes (based on the free escape of gas bubbles from essentially immobile magma) and the Plinian regime (based on the immobility of bubbles relative to the magmatic liquid), precluding the study of the frequent transitions between these regimes that are observed during a single eruptive event.9 These considerations led to the introduction of a second set of Navier–Stokes equations for the gas phase, thus constituting a two-phase system with mass exchange.10–14 In such two-phase systems, the other phase represents an ideal mixture of crystals suspended into a silicate liquid. The effect of mass exchange is far from trivial because it reinforces the effect of gas expansion during ascent, causing gas volume fraction to span the entire range between 0 and 1. Note that the nucleation of gas bubbles, which is not directly addressed in this work, owes a lot of its complexity to the presence of crystals that act as nucleation sites (e.g., Ref. 15). Since the pioneering work of Sparks16, the growth dynamics of bubbles nucleated in viscous liquids has been explored in depth experimentally (e.g., among our works: Refs. 17, 18) and with numerical models (e.g., Refs. 19–23). These studies demonstrated that growth occurs under the combined effect of volatile diffusion from the liquid into the bubble and of the viscous relaxation of the spherical bubble due to gas compressibility. Integrating these two microscopic growth processes into a macroscopic description of the flowing magma has been surprising difficult. Leaving viscous relaxation aside, Mason et al.24 incorporated diffusion degassing with an analytical diffusion equation valid at small Peclet number into a single-phase flow and Lensky et al.25 included a quasi-static diffusion equation into a stick-slip model of shallow magma plug. Instead, most models assume thermodynamical equilibrium bubble growth, whether with a single phase magma (e.g., Refs. 26–28), or a two-phase magma (e.g., Refs. 11, 12). Such growth neglects viscous relaxation by using a single pressure that represents both gas and liquid pressures (the algebraic closure of the review paper (Section 4.2.1)29 and in Part I (Section 2)30)and assumes that the dissolved volatiles (usually water) are always at saturation. As a result, the dissolved volatile concentration is algebraically linked to pressure by, for instance, Henry’s law of solubility. Although in two-phase models, the gas, liquid, and total masses are conserved, this last assumption precludes mass conservation of the volatile species. An example of a two-phase flow model with a single pressure and assuming equilibrium bubble growth is that of Degruyter et al.12 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 4of54 BURGISSER et al. We use the Degruyter model as a reference for comparison because it is regularly applied to natural cases31–33 and it has also been used as a benchmark to establish a new two-phase, isothermal model.34 Some warnings35–37 have been issued about the 1.5D assumption, which considers the conduit as a vertical cylinder. The transition to 2D, however, has proven to be difficult. Only a few models include gas mass exchange to date: Dufek and Bergantz38, Dartevelle and Valentine,39 Collier and Neuberg,40 and Collombet.41 Such 2D models have the capability of studying lateral distribution of parameters such as pressure gradients, gas and heat loss at conduit walls (such as done by Wong and Segall42 in 1D for gas losses), and rheology variations caused by shear zones and shear heating. Our aim was to establish in Part I30 a two-phase, three-component transient conduit flow model that ensures (1) mass and volatile species conservation, (2) disequilibrium degassing taking into account both viscous relaxation and volatile diffusion, and (3) dissipation of total energy. We did not enforce total energy conservation as some terms are negligible and thus simplify the system to be solved. The requirement of energy dissipation corresponds to the fact that we consider the magma without porous surrounding. In this work, we propose a 1.5D simplification of this model. The dimension reduction is done by assuming a parabolic Newtonian velocity distribution horizontally and it is thus sometimes referred to as a 1.5D model. We then propose a numerical implementation and a comparison with the steady-state Degruyter model that highlights to which degree our transient model recovers equilibrium degassing. We finally present a new transient behavior that emerges as a result of relaxing the dual assumption of disequilibrium degassing. 2GENERAL MODEL Here, we present the 3D formulation of the conduit flow model that is based on that presented in Part I, Section 2 of Narbona-Reina et al,30 which is a two-phase, three-component model for a gas and liquid mixture based on mass, momentum, and energy conservation for each phase. The subscripts 𝑙and 𝑔denote the property for each phase, liquid and gas, respectively. The gas phase has only one component, water vapor. The liquid phase has two components, (dry) silicate melt and dissolved water. The unknowns of the system are the gas volume fraction 𝜑, the gas density 𝜌𝑔, the water concentration in the liquid 𝐶𝑙, the velocities 𝑢𝑙,𝑢 𝑔, the pressures 𝑝𝑙,𝑝 𝑔, and the temperatures 𝑇𝑙,𝑇 𝑔(see Table 1). As the scaling of the model against parameter ranges expected in nature has shown in Part I, Section 4 of Narbona-Reina et al,30 some simplifications are acceptable. We thus neglected the effects of surface tension and heat diffusion within each phase. The gas phase is considered inviscid in its own momentum equation but exerts a viscous drag on the liquid. We also assume a fixed bubble number density, 𝑁, and a constant liquid density, 𝜌𝑙.Thislast assumption is common (e.g., Refs. 12, 14) and can be justified by the fact that the liquid density is a weak function of 𝐶𝑙. For instance, the density of a dacitic liquid increases by only 1.3–1.6% for each weight percent of water added.43 Magma fragmentation is not taken into account as we focus on the dynamics of effusive gas escape. This limits model application to cases where 𝜑is smaller than ∼0.8.27 Under these assumptions, the mass, momentum, and heat balances of our system are: 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 5of54 TABLE 1 Main variables list. The constant parameters are listed in Table 2. (a) Latin symbols Variable Definition (unit) 𝑇𝑔Gas temperature (K) 𝑇𝑙Liquid temperature (K) 𝑅𝐻2𝑂Mass exchange term (kg/m3s) 𝐶𝑙Melt water concentration 𝑐𝑝𝑚,𝑐 𝑝𝑛 Bulk heat capacities (J/kg K) 𝑐𝑣𝑔 Gas heat capacity at constant volume (J/kg K) 𝑝𝑔Gas pressure (Pa) 𝑝𝑙Liquid pressure (Pa) 𝑢𝑔Gas velocity (m/s) 𝑢𝑙Liquid velocity (m/s) 𝑢𝑇Bulk heat transport velocity based on 𝑐𝑝𝑚 (m/s) 𝑢𝑛Bulk heat transport velocity based on 𝑐𝑝𝑛 (m/s) 𝑔Gas viscous stress tensor (Pa) 𝑙Liquid viscous stress tensor (Pa) 𝐷𝐻Shear heating (Pa/s) 𝐾𝑑Drag factor (Pa s/m2) 𝑘Magma permeability (m2) 𝑁𝑢 Nusselt number 𝑅𝑏Bubble radius (m) 𝑄Mass discharge rate (kg/s) (b) Greek symbols Variable Definition (unit) 𝜌𝑔Gas density (kg/m3) 𝜑Gas volume fraction or porosity 𝜂𝑙Liquid viscosity (Pa s) 𝛾Heat transfer coefficient (J/ K m3s) 𝜒Viscous factor (1/ Pa s) 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 6of54 BURGISSER et al. The viscous stress tensors for the liquid and gas phases are written as 𝑙=𝜂 𝑙(1 3div𝑢𝑙Id + 2𝐷(𝑢𝑙))𝐷(𝑢𝑙)=1 2(∇𝑢𝑙+(∇𝑢 𝑙)𝑇), (2) 𝑔=0. (3) The term 𝐷𝐻is the shear heating. Its expression corresponds to 𝑢𝑙times the viscous stress term div((1 − 𝜑)𝑙)(see Part I, Section 2.5 of Narbona-Reina et al30): 𝐷𝐻=(1−𝜑)𝜂 𝑙(1 3(div𝑢𝑙)2+2𝐷(𝑢 𝑙)∶∇𝑢 𝑙). For this system, the calculated variables are the liquid pressure and the gas density given by, respectively, a simplified version of the Rayleigh–Plesset equation and the ideal gas relationship: 𝑝𝑙=𝑝 𝑔−1 𝜑𝜒(div𝑢𝑙+𝑅𝐻2𝑂 𝜌𝑙(1 − 𝜑)),(4) 𝜌𝑔=𝑝𝑔 𝑐0𝑇𝑔,(5) where 𝑐0is the ideal gas constant for water. The gravity vector is 𝑔=(0,0, 𝑔) (  𝑔 = −9.81 m∕s2) and 𝜒= 3 4𝜂𝑙,𝐾 𝑑=𝜂𝑔 𝑘,𝜌=𝜑𝜌 𝑔+(1−𝜑)𝜌 𝑙,(6a) where 𝑘is the permeability of the magma. Many complexities of how gas bubbles coalesce to yield an interconnected permeable network remain unclear, but a number of constitutive equations relate gas permeability to other geometrical variables of the bubble network. In the original Degruyter model, the permeability, 𝑘, entering the drag coefficient, 𝐾𝑑, depends on both 𝜑and the bubble radius, 𝑅𝑏. As a result, 𝐾𝑑depends on the bubble number density, 𝑁, because 𝑅3 𝑏=3𝜑 4𝜋𝑁(1 − 𝜑). We used instead a permeability relationship valid for effusive eruptions that depends only on 𝜑44: 𝑘=10 −17(100𝜑)3.4.(6b) This choice removes one variable, 𝑁, from the intermodel comparison. Mueller et al.44 determined that the exponent in (6b) varies between 3 and 3.8; we chose the average value of 3.4 for most runs. The intermodel comparison focuses on a generic case of andesitic volcano. The liquid viscosity, 𝜂𝑙, is thus that of a rhyolitic melt45 (see Box 1for a general explanation on empirical relationships): log10(𝜂𝑙) = −3.545 + 0.833 log10(100𝐶𝑙)+ 9601 − 2368 log10(100𝐶𝑙) 𝑇𝑙− 195.7 − 32.25 log10(100𝐶𝑙),(6c) and the exchange term 𝑅𝐻2𝑂is defined following Part I, Section 2.2.2 of Narbona-Reina et al.30: 𝑅𝐻2𝑂=𝜌 𝑙𝛼𝜑2∕3(1 − 𝜑)2∕3 1−𝜑 1∕3 with 𝛼=3 1∕3(4𝜋𝑁)2∕3𝐷(𝐶𝑙−𝑘 ℎ√𝑝𝑔), (6d) 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 7of54 where 𝐷is the diffusion coefficient of water in the silicate liquid and 𝑘ℎis Henry’s constant. In the temperature equations, 𝑐𝑝𝑔 and 𝑐𝑝𝑙 are the constant heat capacity of gas and liquid, respectively, and 𝛾is the coefficient of heat transfer between gas and liquid. We chose for 𝛾a parameterization often used in gas–particle systems that depends on the Nusselt number Nu, on the local bubble radius 𝑅𝑏and on the constant heat conductivity of the liquid46: 𝛾=24𝜅𝑙𝜑Nu 𝑅2 𝑏 , Nu=7−10(1−𝜑)+5(1−𝜑) 2,𝑅 3 𝑏=3𝜑 4𝜋𝑁(1 − 𝜑).(6e) BOX 1 On empirical and semiempirical relationships Empirical and semiempirical relationships are frequently used in fluid mechanics and thermodynamics when either the processes are not well understood or when the effects of many parameters are lumped together into fitting coefficients. While empirical relationships have a random accuracy beyond their calibrated range and may not ensure a correct behavior when derived, semiempirical relationships are often build to recover some physically sound limit behavior. In other words, extrapolating an empirical relationship is hazardous whereas extrapolating a semiempirical relationships might yield meaningful trends and/or results. Both of these types of relationship excel at performing interpolation between (generally) experimental data points and they are accompanied by an estimate of their accuracy. The fitting coefficients are customarily given with a precision that is commensurate with the level of accuracy of the relationship. In other words, truncating the precision of the coefficients would most commonly cause the relationship to fail in delivering the promised accuracy. Leaving aside the semiempirical ideal gas law (4), we use herein two empirical relationships that govern, respectively, liquid viscosity and water diffusivity. Liquid viscosity. As reviewed by Russell et al.48, the simplest temperature dependence shown by silicate liquid viscosity is Arrhenian: 𝜂𝑙=𝐴 𝐴+𝑒 (𝐸𝐴 𝑅𝐴𝑇𝑙), where 𝑅𝐴is the universal gas constant, 𝐸𝐴is the activation energy for viscous flow, and 𝐴𝐴represents the viscosity at infinite temperature. The parameters 𝐴𝐴and 𝐸𝐴are commonly adjusted to fit the behavior of individual liquid compositions, although the theory suggests that 𝐴𝐴is a constant independent of liquid composition. The Vogel–Fulcher– Tammann (VFT) equation is modification of the Arrhenian relationship that in practice fits better viscosity data: log10(𝜂𝑙)=𝐴 𝑉+𝐵𝑉 𝑇𝑙−𝐶 𝑉, where 𝐴𝑉,𝐵 𝑉,𝐶 𝑉are adjustable parameters specific to individual liquids. As before, 𝐴𝑉determines the viscosity at infinite temperature and 𝐵𝑉is a pseudoactivation energy associated with the structural rearrangement of the liquid. The 𝐶𝑉parameter is the temperature at which viscosity becomes infinite. The composition of magmatic liquids is a function of 12 major components (SiO2, TiO2,Al 2O3,FeO,MnO,MgO, 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 8of54 BURGISSER et al. CaO, Na2O, K2O, P2O5,F,andH2O). In our work, 𝐶𝑙is H2Oand 1−𝐶 𝑙represents all the other components. The Hess and Dingwell45 relationship (6c) singles out the H2Ocompo- nent and thus expands the 𝐴𝑉,𝐵 𝑉,𝐶 𝑉parameters to make the dependence on 𝐶𝑙explicit. One of the most general 𝜂𝑙parameterization is that of Giordano et al.49, which is based on the VFT equation but has lengthy expressions for the three fitting parameters that each depends on the 12 major components listed above. Its average misfit is only ±0.4 log units for a range of liquid composition covering most terrestrial magmas. Water diffusivity. In Section 3, we use the semiempirical relationship of Ni and Zhang50 to describe how 𝐷changes as a function of 𝑝𝑔,𝑇𝑙,and𝐶𝑙. Once dissolved into a silicate liquid, total water (H2Ot) dissociates into two species, molecular water (H2Om) and hydroxyl (OH). H2Omis the dominating diffusing species, and OH is almost immobile. Neglecting OH yields a relationship between the total and molecular diffusivities: 𝐷H2Ot=𝐷 H2Om 𝑑𝑋𝑚 𝑑𝑋𝑙 , where 𝑋𝑚and 𝑋𝑙are the mole fractions of H2Otand H2Om, respectively. The H2Om diffusivity depends exponentially on the H2Otcontent: 𝐷H2Om=𝐷 0𝑒𝐴𝐷𝑋𝑙, where 𝐷0and 𝐴𝐷are parameters depending on temperature and pressure. As 𝑑𝑋𝑚 𝑑𝑋𝑙 has a lengthy expression, Ni and Zhang50 recast these two equations into a single, simplified fit of 𝐷H2Ot, which we call 𝐷Ni in Section 3. The Nusselt number correlation is that of Gunn47 and it is applicable to gas volume fractions 0–0.65 at small Reynolds number based on the relative velocity 𝑢𝑔−𝑢 𝑙. 2.1 Equivalent equations for numerical resolution There are three types of terms that would cause difficulties in the numerical implementation of the system (1) (see Part I, Section 2 of Narbona-Reina et al30): The first type of term is the presence of 𝑝𝑙in (4) which is given through a dilatancy-like constraint, the second is the presence of both 𝜕𝑡𝑝𝑔and 𝜕𝑡𝑇𝑔in (1f), and the third is the temperature coupling through the heat exchange factor 𝛾in (1f)–(1g). Both the liquid momentum equation (1d) and the dilatancy relationship (4) relate the velocity 𝑢𝑙to the two pressures in the system, 𝑝𝑙and 𝑝𝑔. These three variables are thus strongly coupled. By replacing the expression for 𝑝𝑙from (4)into(1d), we simplify the system, resulting in only 𝑢𝑙 and 𝑝𝑔appearing in the momentum equation. As a result of this partial decoupling, 𝑝𝑙becomes an explicitly given variable that is determined by the dilatancy relationship in terms of the other unknowns. We use the expression of the liquid pressure in (4) and the definition above 𝜒= 3 4𝜂𝑙 to write (1 − 𝜑)𝑝𝑙Id = ((1 − 𝜑)𝑝𝑔−𝑅𝐻2𝑂 𝜌𝑙𝜑𝜒)Id − 4(1 − 𝜑)𝜂𝑙 3𝜑 div𝑢𝑙Id. (7) 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 9of54 We inject this expression in the momentum equation (1d) related to 𝑢𝑙, which yields 𝜕𝑡((1 − 𝜑)𝜌𝑙𝑢𝑙)+div((1−𝜑)𝜌 𝑙𝑢𝑙⊗𝑢 𝑙)−div((𝑅𝐻2𝑂 𝜌𝑙𝜑𝜒 +4(1 − 𝜑)𝜂𝑙 3𝜑 div𝑢𝑙)Id + (1 − 𝜑)𝑙) +(1−𝜑)∇𝑝 𝑔−𝐾 𝑑𝜑(1 − 𝜑)(𝑢𝑔−𝑢 𝑙) = (1 − 𝜑)𝜌𝑙𝑔−𝑢𝑔+𝑢 𝑙 2𝑅𝐻2𝑂. (8) This equation has the advantage that 𝑝𝑙is no longer present in it but that it satisfies the condition on 𝑝𝑙set by (4). This substitution has introduced a dissipative term for the velocity in (8), which contributes to the numerical stability of the equation. We addressed the last two types of potential sources of numerical instability with the same goal of writing a simpler system to solve. The presence of two time derivatives in the same equation makes the stability conditions related to the time step size more difficult to manage (we used a variable-order, variable-step-size backward differentiation formula, BDF, Appendix A.2). This simplification is quite straightforward because both quantities, 𝑇𝑔and 𝑝𝑔, are related by the ideal gas law. Finally, the term 𝛾(𝑇𝑔−𝑇 𝑙)may create numerical issues if the magnitudes of the terms in the equation differ significantly, requiring greater precision in the calculation. The coefficient 𝛾takes values on the order of 1010, so when the other terms in the equation are smaller (𝐷𝐻and 𝑐𝑝𝑙(1 − 𝜑)𝜌𝑙𝜕𝑡𝑇𝑙typically span 101–105and 102–103, respectively), this relationship becomes a relaxation equation for 𝛾(𝑇𝑔−𝑇 𝑙), resulting in 𝑇𝑙∼𝑇 𝑔. Eliminating the exchange term from one of the temperature equations has thus the dual advantage of minimizing rounding errors on that equation and of easing the derivation of a single temperature system (see Section 2.2). Another set of equivalent equations to replace (1f)–(1g) can be obtained in terms of 𝜌𝑔,𝑝 𝑔,𝑇 𝑔 from (1f) when using (1a) and the gas law (5). First, we write 𝜕𝑡(𝜑𝑐0𝜌𝑔𝑇𝑔)+div(𝜑𝑐 0𝜌𝑔𝑇𝑔𝑢𝑔)−𝜑(𝜕 𝑡𝑝𝑔+𝑢 𝑔⋅∇𝑝 𝑔)=𝜕 𝑡(𝜑𝑝𝑔)+div(𝜑𝑝 𝑔𝑢𝑔)−𝜑(𝜕 𝑡𝑝𝑔+𝑢 𝑔⋅∇𝑝 𝑔) =𝑝 𝑔(𝜕𝑡𝜑+div(𝜑𝑢 𝑔)) =−𝑝 𝑔𝜕𝑡(1 − 𝜑) + 𝑝𝑔div(𝜑𝑢𝑔) =𝑝𝑔 𝜌𝑙 𝑅𝐻2𝑂+𝑝 𝑔div(𝜑𝑢𝑔+(1−𝜑)𝑢 𝑙). (9) Thus, 𝑐0(𝜕𝑡(𝜑𝜌𝑔𝑇𝑔)+div(𝜑𝜌 𝑔𝑇𝑔𝑢𝑔))−𝑝𝑔 𝜌𝑙𝑅𝐻2𝑂−𝑝 𝑔div(𝜑𝑢𝑔+(1−𝜑)𝑢 𝑙)=𝜑(𝜕 𝑡𝑝𝑔+𝑢 𝑔⋅∇𝑝 𝑔). (10) Equation (10) provides a decomposition of the term 𝜑(𝜕𝑡𝑝𝑔+𝑢 𝑔⋅∇𝑝 𝑔), which causes temperature changes due to gas expansion/contraction. This decomposition will yield an alternative equation for the gas temperature that is more amenable to numerical resolution because it will contain only one time derivative on 𝑇𝑔. The inconvenience is that the expansion/contraction has a clear physical meaning, which is not the case of individual terms of the decomposition (e.g., 𝑝𝑔div(𝜑𝑢𝑔+(1−𝜑)𝑢 𝑙)). Embedding the left-hand-side terms of (10)in(1f)yieldsanalternate 𝑇𝑔equation with the time derivative of a single variable (using the definition of heat capacity at constant volume 𝑐𝑣𝑔 =𝑐 𝑝𝑔 −𝑐 0for conciseness): 𝑐𝑣𝑔𝜑𝜌𝑔(𝜕𝑡𝑇𝑔+𝑢 𝑔⋅∇𝑇 𝑔)+𝑝 𝑔div(𝜑𝑢𝑔+(1−𝜑)𝑢 𝑙)=𝛾(𝑇 𝑔−𝑇 𝑙)−𝑝𝑔 𝜌𝑙𝑅𝐻2𝑂−𝑐 𝑣𝑔𝑇𝑔𝑅𝐻2𝑂.(11) 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 16 of 54 BURGISSER et al. 5.1 Initial conditions Initial values must be prescribed for the system (16); we will use the superscript 𝑖to denote them. The initial guess on the distribution 𝑝𝑖 𝑔is linear and corresponds to the pressure gradient of a bubble-free magma: 𝑝𝑖 𝑔=𝑝 top −𝜌 𝑙 𝑔(𝐻 − 𝑧), (34) where 𝑝top is the fixed pressure at the top boundary of the conduit. There are different ways to choose the initial value for 𝜑𝑖. Ideally, 𝜑can be initialized with a steady-state profile, such as the output of the Degruyter model. Here, to avoid circular reasoning, we specify the initial 𝜑distribution as a linear function of depth: 𝜑𝑖=𝜑 0𝑧 𝐻+𝜑 1,(35) where 𝜑0and 𝜑1are constants. These initial guesses of 𝑝𝑖 𝑔and 𝜑𝑖must be compatible with each other and with the algebraic equations in (16). To minimize incompatibility, we use the same approach as in the analysis of the relaxed system in Part I, Section 5 of Narbona-Reina et al,30 by canceling mass exchange (𝑅𝐻2𝑂=0) and pressure difference (16i). Together with (16h), these assumptions yield: 𝐶𝑖 𝑙=𝑘 ℎ√𝑝𝑖 𝑔,𝑝 𝑖 𝑙=𝑝 𝑖 𝑔,𝜌 𝑖 𝑔=𝑝𝑖 𝑔 𝑐0𝑇𝑖,(36) with 𝑇𝑖a fixed constant temperature 𝑇in, 𝑇𝑖=𝑇 in.(37) When two temperatures are considered, 𝑇𝑖 𝑔=𝑇 𝑖 𝑙=𝑇 𝑖. Finally, 𝑢𝑖 𝑙=𝑢 𝑖 𝑔=𝑢 𝑖to minimize the coupling of the momentum equations by drag. Hence, (𝐶𝑙)|𝑡=0 =𝐶 𝑖 𝑙;(𝜑) |𝑡=0 =𝜑 𝑖;(𝑢 𝑙)|𝑡=0 =(𝑢 𝑔)|𝑡=0 =𝑢 𝑖;(𝑇 𝑙)|𝑡=0 =(𝑇 𝑔)|𝑡=0 =𝑇 in; (𝜌𝑔)|𝑡=0 =𝜌 𝑖 𝑔;(𝑝 𝑙)|𝑡=0 =(𝑝 𝑔)|𝑡=0 =𝑝 𝑖 𝑔.(38) All computational values will be specified later in Table 2. 5.2 Boundary conditions The 1.5D model presented here does not include thermal conductivity. Also, the viscous part of the gas stress tensor, 𝑔, is neglected (see Equation (3)) and the liquid viscosity tensor, 𝑙, is replaced by the damping term  𝛽𝑢𝑙(see Equation (19)). So even if there is a viscous component in Equation (16d) coming from the liquid pressure term through (16i) (see also Equation A.3d): 𝜕𝑧((1 − 𝜑)𝑝𝑙)=𝜕 𝑧((1 − 𝜑)𝑝𝑔−4𝜂𝑙 3𝜑 𝑅𝐻2𝑂 𝜌𝑙−4𝜂𝑙(1 − 𝜑) 3𝜑 𝜕𝑧𝑢𝑙),(39) following the analysis in Part I, Section 3 of Narbona-Reina et al,30 boundary conditions must be imposed in the same fashion as if the system were inviscid. 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 17 of 54 For the 1D case, at the inflow boundary there must be five conditions (two for gas, three for liquid) with one additional condition on 𝜑if it is not already considered for the liquid. At the outflow boundary, two conditions (one for gas, one for liquid) must be set (see table 5 in Part I, Section 3 of Narbona-Reina et al30). We consider the options specified in eq. (88) in Part I of Narbona-Reina et al.30 for this type of systems that we reproduce here: ∙For the inlet boundary, two options are considered by imposing the velocity (BCU) or the pressure (BCP) value. For both of them, fixed values of the concentration 𝐶𝑙and the volume fraction 𝜑must be imposed. For BCU type, velocities must be given such that 𝑢𝑙>0,𝑢 𝑔>0andwehave to prescribe the value of the liquid temperature and to choose between the density 𝜌𝑔or the temperature 𝑇𝑔for the gas phase. For BCP type, pressures and temperatures are imposed for the two phases. ∙For the outlet boundary, only the pressures have to be defined. In the case of BCP, notice that the pressure values must be set to ensure that the flux enters the domain at the inlet boundary and there is no reflection at the outlet boundary. The specific values of the parameters introduced in this section are in Table 2. 5.2.1 Boundary conditions at the inlet (Γin ∶𝑧=0) We use BCP type to impose an inflow condition at the inlet. We must prescribe values for the pressures 𝑝𝑔and 𝑝𝑙and the temperatures 𝑇𝑙and 𝑇𝑔together with the concentration 𝐶𝑙and volume fraction 𝜑. We implemented the numerical solution by keeping conservative variables for Equations (16a), (16b), and (16c) (Appendix A). So, instead of 𝐶𝑙, we must prescribe 𝜌𝑙(1 − 𝜑)(1 − 𝐶𝑙). Similarly, 𝜑𝜌𝑔and 𝜌𝑙(1 − 𝜑) will be prescribed instead of 𝜑and 𝜌𝑔individually. Thus, compatibility with the system and the initial conditions must be carefully looked over. The boundary conditions read: (𝑝𝑙)|Γin =(𝑝 𝑔)|Γin =𝑝 in;(𝑇 𝑙)|Γin =(𝑇 𝑔)|Γin =𝑇 in;(𝜑𝜌 𝑔)|Γin =𝜑 in𝑝in∕(𝑐0𝑇in); (𝜌𝑙(1 − 𝜑)(1 − 𝐶𝑙))|Γin =𝜌 𝑙(1 − 𝜑in)(1 − (𝐶𝑙)in); (𝜌𝑙(1 − 𝜑))|Γin =𝜌 𝑙(1 − 𝜑in), (40) where the following variables have the same values as for the initial conditions: 𝑝in =𝑝 𝑖 𝑔,𝜑in =𝜑 𝑖, and (𝐶𝑙)in =𝐶 𝑖 𝑙. The boundary condition of 𝜑𝜌𝑔has been set to ensure compatibility with (16h)as in (36). These BCP conditions are more favorable than the BCU type for the inlet because setting (𝐶𝑙)in to 𝐶𝑒𝑞 𝑙is the only way to ensure that 𝑅𝐻2𝑂=0if 𝜑>0. While this is straightforward for the BCP type, it creates difficulties for BCU because the corresponding equilibrium porosity (32) is the porosity such that 𝐶𝑙=𝐶 𝑒𝑞 𝑙. So both 𝐶𝑒𝑞 𝑙and 𝜑𝑒𝑞 depends on 𝑝𝑔. To illustrate why it is not possible to leave 𝑝𝑔unspecified and ensure 𝑅𝐻2𝑂=0in BCU, let us consider an inflow condition for the two velocities and two specified temperatures (see Part I, Section 3 of Narbona-Reina et al30 for a justification of this starting choice and a wider discussion). As the drag term is a function of 𝜑, setting equal velocities is a natural choice to maximize decoupling between the two momentum equations. As a result, the total water content 𝐶𝑇will be a transported quantity at the boundary (see Equation 30). Thus, setting 𝐶𝑇at the boundary would be a compatible, stable choice. The corresponding, alternate equation solved, however, is the 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 18 of 54 BURGISSER et al. conservation of 𝜌𝑙(1 − 𝜑)(1 − 𝐶𝑙)(16a), which is equivalent to 𝜌(1 − 𝐶𝑇). While this equation does not depend on any source term, it depends on 𝑝𝑔through the product (1 − 𝜑)(1 − 𝐶𝑙). In addition, the most common choice for gas systems is to impose the gas density 𝜑𝜌𝑔, which is a product that also depends on 𝑝𝑔. Eliminating 𝑝𝑔from these products is not feasible algebraically because the three conditions (1 − 𝜑)(1 − 𝐶𝑙)=𝐶 1,𝜑𝜌𝑔=𝐶 2,and𝐶𝑙=𝑘 ℎ√𝑝𝑔(where 𝐶1,𝐶 2are set constants) yield a third-degree equation on 𝑝𝑔, precluding the option of letting it free. Switching (16a) for (30) would not alleviate that issue. In other words, BCU implies that 𝑅𝐻2𝑂≠0, which in our experience can only be stabilized by an additional numerical boundary condition (see below). An additional limitation of letting 𝑝𝑔free is that 𝑝𝑔cannot decrease below the critical value (𝐶𝑒𝑞 𝑙∕𝑘ℎ)2because then 𝜑𝑒𝑞 <0. Another advantage of BCP conditions is that if the inlet pressure is a prescribed function of time, such as 𝑝in(𝑡), an algorithm preserving 𝑅𝐻2𝑂=0at the inlet at all times can be implemented by calculating 𝐶𝑇(𝑡0)at the initial time 𝑡0using (29)with𝜑in(𝑡0). Then, as 𝑡>𝑡 0,𝜌𝑔is updated with 𝑝in(𝑡) and (16h), (𝐶𝑙)in is updated with (28), and 𝜑in is updated using (32). Finally, the conservative variables are updated according to (40). This implementation can be used to couple the volcanic conduit flow to a deeper reservoir source (e.g., Refs. 27, 52). 5.2.2 Boundary conditions at the outlet boundary (Γout ∶𝑧=𝐻) To avoid nonphysical results when 𝑝top is atmospheric (see Section 6.1.3) and thus ensure a realistic model intercomparison, we imposed pressure equality in most runs by using: (𝑝𝑔)|Γout =(𝑝 𝑙)|Γout =𝑝 top.(41a) A more accurate outflow boundary condition is to consider that the gas phase equilibrates with the atmosphere while 𝑝𝑙follows the dilation relationship (16i): (𝑝𝑔)|Γout =𝑝 top,(𝑝 𝑙)|Γout =𝑝 top −(1 𝜒𝜑𝜌𝑙)|Γout ((𝜌𝑙𝜕𝑧𝑢𝑙)|Γout +(𝑅 𝐻2𝑂)|Γout )(41b) with (𝑅𝐻2𝑂)|Γout =3 1∕3(4𝜋𝑁𝜑(1 − 𝜑))2∕3 |Γout 𝐷𝜌𝑙 (𝐶𝑙)|Γout −𝑘 ℎ√𝑝top 1−(𝜑) 1∕3 |Γout . Finally, a third possible outflow boundary condition is to consider the mixture momentum equation (see Part I, Section 4 of Narbona-Reina et al30) and impose an average pressure: (𝜑𝑝𝑔+(1−𝜑)𝑝 𝑙)|Γout =𝑝 top.(41c) Thanks to (16i), this condition reads (𝑝𝑔)|Γout =𝑝 top +(1 𝜒𝜑𝜌𝑙)|Γout (((1−𝜑 )𝜌𝑙𝜕𝑧𝑢𝑙)|Γout +(𝑅𝐻2𝑂)|Γout ) and (𝑝𝑙)|Γout =𝑝 top −(1 𝜒(1−𝜑 )𝜌𝑙)|Γout (((1−𝜑 )𝜌𝑙𝜕𝑧𝑢𝑙)|Γout +(𝑅𝐻2𝑂)|Γout ). 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 19 of 54 In case, the velocity is incoming at this boundary, that is, 𝑢𝑙⋅ 𝑛<0or 𝑢𝑔⋅ 𝑛<0, we should impose the values for the other variables. Similarly than before, they are given by (𝑇𝑙)|Γout =𝑇 𝑖,(𝜌 𝑙(1 − 𝜑)(1 − 𝐶𝑙))|Γout =𝜌 𝑙(1 − 𝜑𝑖)(1 − 𝑘ℎ√𝑝top), (𝜌𝑙(1 − 𝜑))|Γout =𝜌 𝑙(1 − 𝜑𝑖)IFF 𝑢𝑙⋅ 𝑛<0, (𝑇𝑔)|Γout =𝑇 𝑖,(𝜑𝜌 𝑔)|Γout =𝜑𝑖𝑝top 𝑐0𝑇𝑖IFF 𝑢𝑔⋅ 𝑛<0. (41d) This case never occurred in the performed numerical simulations even though it was implemented. 5.2.3 Additional numerical boundary conditions As discussed in the literature,53–55 the physical and mathematical boundary conditions are not always enough to solve the problem numerically and additional numerical (also called “soft”) boundary conditions may be considered. For example, Poinsot and Lele54 considered the case of imposing both velocity and temperature at the inlet, and found that an additional numerical condition on the density must be added to keep compatibility with its corresponding equation. In most cases, these numerical conditions consist of an extrapolation of the concerned variables at the boundary, that it is interpreted as a zero normal derivative through this boundary (Neumann condition). The restrictive outlet boundary condition of equal pressures (41a) involves numerical boundary conditions because of the peculiar closure expression for the liquid pressure (16i). As pointed out at the beginning of this section, the liquid pressure involve an equivalent diffusion term for the liquid velocity (Equation 39, see also Equation A.3d). For such equation type, we usually have to prescribe a Neumann condition to solve the second derivative of the velocity. In our case, this condition comes from the liquid pressure since considering 𝑝𝑙=𝑝 𝑔at the boundaries means that the diffusive term (39) is zero at the boundary: (4𝜂𝑙 3𝜑 𝑅𝐻2𝑂 𝜌𝑙+4𝜂𝑙 3𝜑 (1 − 𝜑)𝜕𝑧𝑢𝑙)|Γ𝑗 =0, for any boundary 𝑗 = in, out. (42) This is the additional numerical condition considered at both boundaries to correctly solve the velocity 𝑢𝑙. Notice that simply imposing 𝜕𝑧𝑢𝑙=0would be incorrect because it yields a wrong liquid pressure value. This condition is not an additional boundary condition needed to close the system but it is needed for algorithmic reasons in the numerical implementation. This condition just duplicates information without overspecifying the problem. 5.3 The Degruyter model The numerical resolution of the Degruyter model was performed using the Octave v6.2 software. We refer the reader to Degruyter et al.12 for details on implementation. It is a two-point boundary value system with 𝑝|Γout =𝑝 top and 𝑝|Γin =𝑝 in =𝜌 𝑙 𝑔𝐻. The algorithm does not need initial values for 𝐶𝑙, as the equivalent variable, 𝐶eq 𝑙,isgivenby(28). Similarly, only an inlet value of 𝜑|Γin =𝜑 1is needed. The inlet velocities are (𝑢𝑙)|Γin =(𝑢 𝑔)|Γin =𝑢 in. The solution is found thanks to a bisection 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 20 of 54 BURGISSER et al. algorithm, so two initial bounding values for 𝑢in are chosen so that they yield one positive and one negative residual, respectively. The bounding values are selected manually through the constant mass flow rate 𝑞: 𝑢in =𝑞 (𝑛in (𝜌𝑔)in +1−𝑛 in 𝜌𝑙)with (𝜌𝑔)in =𝑐0𝑇 𝑝in ,𝑛 in =𝐶𝑐−𝐶 eq 𝑙 1−𝐶 eq 𝑙 , 𝐶𝑐=𝜑1(𝜌𝑔)in +(1−𝜑 1)𝜌𝑙𝐶eq 𝑙 𝜑1(𝜌𝑔)in +(1−𝜑 1)𝜌𝑙 , where 𝑇is a constant temperature. 6RESULTS Table 2lists the constant parameters used in the model intercomparison. These values are typical of an effusive eruption at an andesitic volcano. In addition to these parameters, several constant values for 𝐷(0, 10−13,and10−8 m2∕s) plus the variable diffusion, 𝐷Ni,givenby(21)–(22)were chosen to cover the range typical of rhyolitic liquids and the end-member of no mass exchange. TABLE 2 Constant parameters used in the intermodel comparison. In addition, 𝐷varies according to (21)in some transient runs. (a) Parameters common to both models Parameter Value Definition (unit) 𝐻5000 Conduit height (m)  𝑔−9.81 Gravity acceleration (m/s2) 𝜌𝑙2450 Liquid density (kg/m3) 𝜂𝑔2×10 −5 Gas viscosity (Pa s) 𝑐0461.4 Gas constant for water (J/kg K) 𝑘ℎ4.11 × 10−6 Solubility coefficient (Pa−1/2) 𝑟𝑐15 Conduit radius (m) 𝑢𝑖0Initial speed (m/s) 𝑇in,𝑇 1138 (Inlet) temperature (K) 𝑝top 0.1 Outlet pressure (MPa) 𝑝in 120 Inlet pressure (MPa) 𝜑15×10 −3 Inlet porosity (b) Parameters only present in the transient model Parameter Value Definition (unit) 𝐷0–10−8 Melt water diffusion (m2/s) 𝑁1012 Bubble number density (m−3) 𝑐𝑝𝑙 1000 Liquid heat capacity (J/kg K) 𝑐𝑝𝑔 2030 Gas heat capacity (J/kg K) 𝜅𝑙2.3 Liquid heat conductivity (W/m K) 𝜑00.1–0.7 Outlet porosity minus 𝜑1 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 21 of 54 FIGURE 1 Gas volume fraction (𝜑) as a function of height in the conduit for the Degruyter run and for the transient model at steady state with several values of diffusion (isothermal) and 𝐷Ni (variable, single temperature). (A) Full conduit length (curves for 𝐷=10 −8,𝐷 Ni, and Degruyter almost overlap). (B) Closeup of the last 200 m of the conduit before the surface. Note that, for water vapor, we used the relationship between the heat capacity 𝑐𝑝𝑔 and the constant 𝑐0given by 𝑐𝑝𝑔∕𝑐0=4.4. All transient runs except one were done by assuming 𝑇𝑔=𝑇 𝑙. When initialized with the values of Table 2and an initial guest for 𝑞, the Degruyter model converges toward a constant mass flow rate of 5.3 × 105kg/s, which corresponds to an inlet velocity of 0.3 m/s. Such an ascent rate corresponds to effusive conditions near the transition between effusive and explosive behavior32. These conditions are consistent with the fact that, near the top boundary (surface), 𝜑nearly reaches 0.8 before decreasing of a few percent when reaching the surface. These conditions, close to magma fragmentation, allowed us to maximize the range of 𝜑values explored in this comparison. We used the steady-state solution of the Degruyter model to obtain the spatial distribution of 𝜑. We then approximated that distribution by adjusting 𝜑0in (35) to obtain the reference initial condition for the transient model. The steady states reached by the transient model will be shown to be independent of 𝜑0. 6.1 Steady state 6.1.1 Constant temperature For clarity, system (1) is referred to as the transient model and system (24) is referred to as the Degruyter model. All steady-state solutions of the transient model were obtained with the COMSOL Multiphysics stationary solver that ignored all the time derivatives of the equation system. The solver was initialized by transient solutions that will be described in Section 6.2. The value of 𝐷significantly changes the steady state reached by the transient model. Figure 1A shows steady-state profiles of gas volume fraction (𝜑, here referred to as porosity) for progressively vanishing values of 𝐷plus the variable case with 𝐷Ni.At𝐷=0m 2∕s, which suppresses mass exchange by setting 𝑅𝐻2𝑂=0, the steady-state porosity evolution is only due to gas expansion and permeable flow. As a result, porosity remains low from depth to 4.5 km, after which it increases 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 22 of 54 BURGISSER et al. sharply to reach ∼90 vol% at the surface. The water dissolved in the melt, 𝐶𝑙, is constant because it is imposed at the base of the conduit and that boundary value is advected to the top without mass exchange. The range of values covered by 𝐷Ni from the base to the top of the conduit is 10−11–10−12 m2∕s, respectively. So when 𝐷≥10−12 m2∕s, the steady-state porosity profiles from the transient model are close to that of the Degruyter model. When focusing on the shallowest part of the conduit (Figure 1B), only the case with 𝐷=10 −8 m2∕s reproduces the porosity reduction near the surface of the Degruyter model. Such a diffusion value is unrealistically high for natural melts50 suggesting that perfect equilibrium degassing during the whole ascent is unattainable. The shallow porosity reduction in the Degruyter model is due to permeable flow that induces 𝑢𝑔>𝑢 𝑙and gas evacuation12. This gas/liquid separation is affected by mass transfer, which is driven by the assumption of equilibrium degassing. In the Degruyter model, equilibrium degassing is enforced by the dual condition of pressure equilibrium, 𝑝𝑔=𝑝 𝑙, and no water supersaturation, 𝐶𝑙=𝐶 eq 𝑙. In the transient model, the first condition is nearly fulfilled regardless of the diffusion value because the gas overpressure with respect to the liquid pressure (𝑝𝑔−𝑝 𝑙)is small everywhere except near the surface, where it remains smaller than 0.003 MPa (Figure 2A). There is a thin boundary layer in the distribution of the gas overpressure near the outlet. It caused by the assumption that both outlet pressures are equal to atmospheric, and Section 6.1.3 presents the consequences of relaxing this assumption. How well the transient model solution fulfills the second condition, water supersaturation (i.e., 𝐶𝑙≠𝐶eq 𝑙), can be measured two ways. One measure is the difference 𝐶𝑙−𝐶 eq 𝑙(Figure 2B). The other is the difference between the saturation pressure (𝐶𝑙∕𝑘ℎ)2and 𝑝𝑔because this measure of supersaturation can be directly related to the likelihood of a (new) nucleation event (Figure 2C). Bubble nucleation is a dynamic process but it is often simplified in the case of heterogeneous nucleation (i.e., bubble nucleation on preexisting crystals) by setting a threshold in pressure drop (>5 MPa for oxides, >30 MPa for biotite56,57) above which the rate of nucleation is large enough to modify 𝑁significantly over the time period of interest. Both measures yield similar trends, so we only focus on the supersaturation pressure for the different values of diffusion (Figure 2C).Thecasewith𝐷Ni has a supersaturation pressure ∼1 MPa over most of the conduit and between 1 and 5 MPa at both conduit ends, possibly triggering a second nucleation event if oxide crystals were present. There is thus more dissolved water in the shallow part of the conduit of the transient model than in that of the Degruyter model. The effect of the separation of gas and liquid driving the shallow porosity reduction is dampened in the transient model by a slower rate of mass transfer compared to the forced equilibrium degassing of the Degruyter model, thus delaying porosity reduction. Provided that 𝐷is high enough (>10 −11 m2∕s), equilibrium degassing can thus be achieved by the transient model over most of the conduit to within 0.012 of 𝜑eq. Conduit ends, however, remain supersaturated with water, although not to a degree able to induce heterogeneous nucleation on any minerals commonly present in intermediate or silicic magmas. The Degruyter model is a two-point boundary value system where the single pressure, 𝑝,isset at both conduit ends. The transient model also sets pressure at both boundaries but it considers two pressures. As they have close values (Figure 2A), it suffices to compare 𝑝𝑔and 𝑝. The interplay between the dynamic pressure controlled by 𝑢𝑙and the magmatic pressure controlled by 𝜑causes 𝑝𝑔(and 𝑝) to depart from a linear pressure profile for both the Degruyter model and the transient model with 𝐷>10 −11 m2∕s (Figure 3). Cases with 𝐷=10 −13 and 𝐷=0m 2∕s have nearly linear pressure profiles because 𝜑is small over most of the conduit length (Figure 1). Both models yield very similar distributions of total water content (Figure 4), velocities (Figure 5), and thus mass discharge rate (5.7 vs. 5.3 ×10 5kg/s). For simplicity, we only use the 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 23 of 54 FIGURE 2 Measures of disequilibrium degassing as a function of height in the conduit for the transient model at steady state with several values of diffusion (isothermal) and 𝐷Ni (variable, single temperature). (A) Gas overpressure (𝑝𝑔−𝑝 𝑙). (B) Oversaturation (𝐶𝑙−𝐶 eq 𝑙). (C) Supersaturation pressure |(𝐶𝑙 𝑘ℎ )2−𝑝 𝑔|. The two vertical lines indicate approximate nucleation thresholds for oxides (red) and biotite (cyan). transient model solution with 𝐷Ni as a reference to characterize how water weight fractions change as a function of height in the conduit. The total water content at the inlet, 𝐶𝑇, is the same in both models (i.e., 𝐶𝑇=𝐶 𝑐, Section 5.3). In the Degruyter model, 𝐶𝑐is assumed constant along the conduit length. Figure 4shows that 𝐶𝑇in both models have the same decreasing trend toward the surface, which indicates that total water content is in fact not constant along the conduit. That water mass conservation is not maintained is also apparent in the behavior of the ratios between exsolved water and total mass, 𝑛and 𝑛eq, because they strongly diverge above 4.5 km. This divergence happens where 𝑢𝑔becomes much larger than 𝑢𝑙, decoupling outgassing from gas exsolution (Figure 5). At the local scale of a single gas bubble, disequilibrium degassing has been characterized by different regimes of bubble growth58. In a Lagrangian framework attached to one bubble, 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 24 of 54 BURGISSER et al. FIGURE 3 Gas pressure (𝑝𝑔)asa function of height in the conduit for the Degruyter model and the transient model at steady state with several values of diffusion. FIGURE 4 Various measures of water weight fraction as a function of height in the conduit for the Degruyter model (𝐶𝑐) and the transient model with 𝐷Ni at steady state (𝐶𝑇,𝑛,𝑛eq). Note that the three corresponding curves of the Degruyter model for 𝐶𝑇,𝑛,and 𝑛eq would overlap exactly those of the transient model at the scale of the figure. these regimes can be identified by two dimensionless parameters, Θ𝐷and Θ𝑉, that are linked, respectively, to the diffusion time scale 𝜏dif =𝑅2 𝑏 𝐷and the viscous time scale 𝜏vis =4𝜂𝑙0 𝑝𝑙0 20,58: Θ𝐷=𝜏dif 𝜏=𝑅2 𝑏 𝜏𝐷 and Θ𝑉=𝜏vis 𝜏=4𝜂𝑙0 𝜏𝑝𝑙0 ,(43) 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 25 of 54 FIGURE 5 Gas (𝑢𝑔) and liquid (𝑢𝑙) velocities as a function of height in the conduit for the Degruyter model and the transient model with 𝐷Ni at steady state with constant and variable, single temperature. where 𝜏=𝑝 𝑙∕ 𝑃is a characteristic time scale of pressure change and  𝑃is the local decompression rate given by (in an Eulerian framework):  𝑃=𝜕 𝑡𝑝𝑙+𝑢 𝑙𝜕𝑧𝑝𝑙.(44) These characteristic times come from the microscopic framework (see also Part I, Section 2.2 of Narbona-Reina et al30); 𝜏dif is the time for the Lagrangian of the concentration (see eq. 26 in Narbona-Reina et al30)and𝜏vis comes from the Rayleigh–Plesset equation for the bubble radius (see eq. 17 in Narbona-Reina et al30). In Part I, Section 4 of Narbona-Reina et al.30, we have proposed a scaling for the mass conservation (16b)–(16c) and the 𝜑equation (50) below in drift-flux formulation. The Lagrangian scaling (43) can be related to our Eulerian scaling by assuming that the pressure changes are dominated by the vertical motion (𝑢𝑙term in Equation 44) so that 𝜏=𝐿 0∕𝑢0, where 𝐿0and 𝑢0are the characteristic length and speed, respectively. In two-phase formulation (i.e., with 𝑢𝑙and 𝑢𝑔), the dimensionless mass conservation (16b)and𝜑transport (50) can be expressed in terms of Θ𝐷 and Θ𝑉: St 𝜕 𝑡(1 − 𝜑) + 𝜕 𝑧((1 − 𝜑)  𝑢𝑙)=−St 3𝐶𝑙0 Θ𝐷(𝜑2∕3  𝑁)2∕3  𝐶𝑙−√ 𝑝𝑔 1−𝜑 1∕3 , St 𝜕 𝑡𝜑+ 𝑢𝑙𝜕 𝑧𝜑=St 3 Θ𝑉 𝜑(1 − 𝜑)  𝜂𝑙(𝑝𝑔0 𝑝𝑙0  𝑝𝑔− 𝑝𝑙), where the diacritic notation with a tilde is used for dimensionless variables (that is, 𝜙=𝜙 0 𝜙where 𝜙0is some reference value) and St =𝐿 0∕(𝑢0𝑡0)is the Strouhal number. The Θ𝐷–Θ𝑉framework of Forestier-Coste et al.20 can thus be used to characterize the bubble growth regimes occurring during ascent in our simulations. 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 32 of 54 BURGISSER et al. FIGURE 10 Gas overpressure with respect to the liquid pressure (𝑝𝑔−𝑝 𝑙) in the upper part of the conduit for the transient model at steady state with 𝐷Ni, a single temperature, and an outlet pressure of 0.1 MPa. The outlet has either an imposed pressure equality (black, 𝑝top =𝑝 𝑔=𝑝 𝑙), or a liquid pressure controlled by dilatancy (cyan, 𝑝top =𝑝 𝑔), or an imposed average pressure (green, 𝑝top =𝜑𝑝 𝑔+(1−𝜑)𝑝 𝑙). (Figure 2A). As near the outlet the liquid is accelerating and gas exsolution occurs, one can expect that (𝑝𝑙)|Γout <𝑝 top. Now, 𝑝top being set to the (rather small) atmospheric value of 0.1 MPa, both acceleration and exsolution cause (𝑝𝑙)|Γout to reach the negative value of −5 MPa (Figure 10). Using the averaged pressure condition with (𝜑𝑝𝑔+(1−𝜑)𝑝 𝑙)|Γout =𝑝 top given by (41c) reduces the pressure difference at the outlet but does not alleviate negative 𝑝𝑙values: (𝑝𝑙)|Γout = −1.26 MPa and (𝑝𝑔)|Γout =0.6MPa. These nonphysical results are linked to the breakdown of the 1.5D conduit geometry at the outlet; a realistic simulation of an effusive eruption with atmospheric gas pressure at the outlet needs to take the free surface of the liquid into account. One way to set realistic outlet conditions in a 1.5D geometry is to consider that the cylindrical conduit is topped by a magma accumulation (such as a dome) that increases the upper boundary pressure. A simple example without magma accumulation above the upper boundary would be a 100-m thick overburden that causes 𝑝top =2.5MPa. Runs done with (𝑝𝑔)|Γout =(𝑝 𝑙)|Γout =2.5 MPa and runs done with the alternate outlet boundary conditions (41b)and(41c) yield again almost identical steady-state results on all the variables but 𝑝𝑙, which has the realistic outlet value of 1.91 MPa with (41b)and2.36MPawith(41b)((𝑝𝑔)|Γout =2.57MPa in that last case). Compared to the variable temperature solutions with no overburden, the addition of an overburden increases the mass discharge rate from 2.9 to 3.4 × 105kg/s regardless of the outlet pressure condition but yields similar trends as those shown in Figures 1–5. 6.1.4 Simplified system In steady state, the largest change induced by our expression of the energy balance is a reduction of the mass discharge rate by a factor 2 compared to the isothermal case. This reduction is caused by the increased viscosity due to liquid cooling. Although temperature variations cannot be ignored, some simplifications do not affect much the steady-state mass discharge rate. We assess the impact of neglecting various physical processes to find the minimum model that still captures the essential physics of steady-state magma ascent. We use the mass discharge rate: 𝑄=𝜋𝑟 2 𝑐(𝜑𝜌𝑔𝑢𝑔+(1−𝜑)𝜌 𝑙𝑢𝑙) 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 33 of 54 as a diagnostic quantity because it can be shown by adding (1a)and(1b) that 𝑄is constant in steady state (Remark 1 in Narbona-Reina et al30). The physical processes considered are shear heating, liquid inertia, mass exchange, and permeable flow. Taking shear heating into account in the temperature equation increases the mass discharge rate by 13 %. A run with the transient model using (46) instead of (16e) increases the mass discharge rate from 2.9 to 3.2 × 105kg/s (8 %). Omitting the 𝑅𝐻2𝑂and inertial terms in the momentum equations (16d)–(16e) and the 𝑅𝐻2𝑂term in the 𝑝𝑙equation (16i) decreases mass discharge rate by only 0.1%. These four physical processes have thus relatively low influence on 𝑄 and thus could potentially be neglected to obtain a simplified system. As 𝜑is generally small at the conduit inlet, we use the scaling variable 𝑄0=𝜌 𝑙0𝜋𝑟2 𝑐𝑢0∼𝑄 |Γ𝑖𝑛. We relax the constraint of energy dissipation and we arbitrarily set that the 𝑄0of the simplified system should remain within 10% of that of the full system. As we have explored only a small subset of conduit conditions (Table 2), we use the dimensionless analysis of Section 4.3 of Narbona-Reina et al30) to design regime diagrams. We start each analysis by assuming that the dimensionless number of interest should be much lower that the equation leading-order term(s). We then use the steady-state solution of the transient model with 𝐷Ni and a single temperature (e.g., Figure 9) to determine the actual threshold value(s) of that dimensionless number beyond which physical processes could be neglected. Finally, the thresholds are used to map the dynamical regimes over a wide range of natural parameters. The dimensionless factor of the shear heating term is 𝐹sh =𝜂 𝑙0𝑢0∕(𝜌𝑙0𝐿0𝑇0𝑐pl). Figure 11A shows that it is well below 1 over the whole ascent despite that shear heating affects mass discharge rate by >10%. To determine the actual threshold below which shear heating can be neglected, we use the condition that the mass discharge rate without shear heating, 𝑄1=𝜌 𝑙0𝜋𝑟2 𝑐𝑢1, should be within 10% of the reference rate 𝑄0: 𝑓=0.1=𝑄1−𝑄 0 𝑄0=𝑢1 𝑢0−1. To the main order (Section 4.3 of Narbona-Reina et al30), the liquid speed is proportional to 𝜂𝑙 because the flow can be approximated by the Poiseuille equation (i.e., an incompressible Newtonian fluid flowing in a laminar fashion through a cylindrical pipe of constant radius). So 𝑢1∕𝑢0≈𝜂 𝑙1∕𝜂𝑙0 and 𝜕𝑇(log10 𝜂𝑙)∼𝜂𝑙1 −𝜂 𝑙0 Δ𝑇 log10 𝑒 𝜂𝑙0 ∼log10 𝑒 Δ𝑇 (𝑢1 𝑢0−1 )=log10 𝑒 Δ𝑇 𝑓 where Δ𝑇 is the temperature difference between the two cases. In other words Δ𝑇 = 𝑓log 10 𝑒 𝜕𝑇(log10 𝜂𝑙) is the maximum temperature difference beyond which liquid heating (or cooling) causes mass discharge rate fluctuation in excess of 100𝑓 %. In the comparison shown on Figure 8,evaluation of Δ𝑇 with 𝑓=0.1and  𝛽yields 3–5 K. This range of Δ𝑇 occurs around 20–30 MPa, which corresponds to an 𝐹𝑠ℎ range of 1.4–2.4 ×10−6 (Figure 11A). The number of degrees of freedom controlling 𝐹𝑠ℎ can be reduced by relating 𝜂𝑙0 to 𝑝0by using (6c)and(28). Using that relationship, Figure 11A shows that two extreme values of 𝑢0frame the behavior of 𝐹𝑠ℎ as a function of pressure (Table 3). As these extreme values are only an order of magnitude apart, 𝐹𝑠ℎ changes during ascent are mostly due to 𝜂𝑙0, the variation of which spans several orders of magnitude. Converting 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 34 of 54 BURGISSER et al. FIGURE 11 Scaling analysis and regime diagrams. (A, C, E, G) Dimensionless parameters as a function of pressure for the transient run with 𝐷Ni and a single temperature (black curves). Individual threshold values are shown as gray symbols or gray horizontal lines. (B, D, F) Threshold liquid viscosity (𝜂𝑙0) as a function of mass discharge rate (𝑄0). (H) Threshold drift speed (𝑤0) as a function of 𝑄0. (A) Shear heating scaling (𝐹𝑠ℎ). Red and blue curves frame the transient run for extreme velocities (𝑢0). (B) The gray area indicates where shear heating can be neglected. (C) Liquid inertia scaling (𝑀𝑎,𝑅𝑒). Red and blue curves frame the transient run for extreme 𝑢0 values. (D) The gray area shows the parameter space where liquid inertia can be neglected. The vertical colored lines represent various pressures (𝑝0) if the pressure gradient dominates the momentum equation (𝑀𝑎). The gray 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 35 of 54 TABLE 3 Scaling variables and their ranges in the transient run. Parameter Minimum Maximum 𝑢00.1 1 𝑤010−9 10 𝜂𝑙0 104108 𝜌𝑙0 2450 2450 𝐿030 30 𝑇01138 1138 𝑝00.1 × 106120 × 106 𝑢0to 𝑄0and using the threshold values of 𝐹𝑠ℎ, Figure 11B shows the maximum 𝜂𝑙0 values below which shear heating can be neglected for a range of natural 𝑄0values. By convention, the scaling analysis yields dimensionless numbers such that the inertial term of the momentum equation is 𝑂(1) (i.e., there is a factor 1 in front of 𝜕 𝑧((1 − 𝜑)  𝜌𝑙 𝑢2 𝑙)in the dimensionless version of Equation 16d). To decide whether this term can be neglected, it needs to be compared to the other terms of the same equation. Focusing on the two largest terms that survive at the main order, the pressure gradient and the viscous resistance term, one can multiply the dimensionless momentum equation (see Section 4.3 and the mixture velocity equation in Appendix E of Narbona-Reina et al30) by either 𝑀𝑎2, where 𝑀𝑎2=𝜌 𝑙0𝑢2 0∕𝑝0is the Mach number, or by 𝑅𝑒, where 𝑅𝑒 = 𝜌𝑙0𝑢0𝐿0∕𝜂𝑙0 is the Reynolds number, respectively, to evaluate whether the inertial term can be neglected. Figure 11C shows that the inertial term can be neglected if 𝑀𝑎2<0.1and 𝑅𝑒 < 5. It also suggests that 𝑀𝑎 changes during ascent mostly because of 𝑝0,𝑢0playing a subordinate role. Using the same reasoning as for 𝐹𝑠ℎ,𝑅𝑒 changes during ascent are mostly due to 𝜂𝑙0 with 𝑢0also playing a minor role in controlling 𝑅𝑒 as a function of pressure. Converting 𝑢0in 𝑄0, Figure 11D suggests that inertial terms can only be neglected if 𝑄<3×10 6kg/s and at progressively lower liquid viscosity below that mass discharge rate. The same reasoning and accompanying multiplication by either 𝑀𝑎2,or𝐹𝐾=𝑘 0𝑢2 0𝜌𝑙0∕ (𝜂𝑔0𝐿0𝑤0)with the weighed drift velocity equation on 𝑢𝑔−𝑢 𝑙(velocity equation on 𝑌𝑤 in Appendix E of Narbona-Reina et al30) yields dimensionless parameters that can be grouped in three distinct factors for each scaling: 𝐹𝑚1 =𝑤 0𝑢0∕(𝑇0𝑐0),𝐹𝑚2 =𝑤 2 0∕(𝑇0𝑐0),and𝐹𝑚3 = 𝜂𝑙0𝑢0∕(𝜌𝑙0𝐿0𝑐0𝑇0)for 𝑀𝑎2and 𝐹𝑘1 =𝑝 0𝑢0𝑘0∕(𝑇0𝑐0𝜂𝑔0𝐿0),𝐹𝑘2 =𝑝 0𝑤0𝑘0∕(𝑇0𝑐0𝜂𝑔0𝐿0),and𝐹𝑘3 = 𝑝0𝜂𝑙0𝑢0𝑘0∕(𝜌𝑙0𝐿2 0𝑐0𝑇0𝜂𝑔0𝑤0)for 𝐹𝐾. Figure 11E shows how these six factors change during ascent. The largest factor is either 𝐹𝑚3, or 𝐹𝑘3, with a common threshold of 0.1. In the case of 𝐹𝑚3, and considering the parameters kept constant in this analysis (Table 3), there is a direct relationship between 𝜂𝑙0 and 𝑄0that is shown in Figure 11F. In the case of 𝐹𝑘3, the maximum viscosity below which Darcy’s law (46) is valid is 𝜂𝑙0 =𝐹 𝑘3 4𝜌2 𝑙0𝑇0𝑐0𝐿4 0𝜂𝑔0 𝜋𝑄0 𝑤0 𝑘0𝑝0 . line is valid if viscous stresses dominate the momentum equation (𝑅𝑒). (E) Darcy scaling. Curves represent the scaling parameters (𝐹𝑖𝑗). (F) Gray areas indicate where Darcy’s equation is valid when viscous resistance (𝐹𝑘3, dark gray) or inertia (𝐹𝑚3, dark and light gray) dominates. (G) Mass exchange scaling (𝐹𝑟). Red and blue curves frame the transient run for extreme 𝑢0and 𝑤0values. (H) The gray area indicates where liquid inertia can be neglected. 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 36 of 54 BURGISSER et al. Considering the parameters kept constant in Table 3and setting 𝐹𝑘3 =0.1,𝜂𝑙0 depends on 𝑄0 and on the dimensional quantity 𝑤0∕(𝑘0𝑝0). Figure 11F shows this maximum viscosity for three representative values of 𝑤0∕(𝑘0𝑝0): 1 (e.g., 𝑤0=10 −6 m∕s,𝑘0=10 −14 m2,𝑝0= 100 MPa), 102 (e.g., 0.01 m∕s,10−12 m2,100 MPa), and 104(e.g., 10 m∕s,10−10 m2,10 MPa). Overall, Figure 11F suggests that the drift velocity equation (46) is valid over the whole range of mass discharge rates when 𝜂𝑙0 <10 4and when 𝜂𝑙0 <10 8and 𝑤0∕(𝑘0𝑝0)≥104. Remark 3. The maximum values of 𝑤0considered in Table 3might imply that the Darcian permeability relationship (6), which only takes viscous effects into account, should be extended to a generalized Forchheimer relationship that includes the inertial effects caused by gas recirculation within the microscopic cavities of the bubble network (e.g., Refs. 12, 62): 𝐾𝑑(𝐹𝑜) =𝜂𝑔 𝑘+𝜌𝑔 𝑘2|𝑢𝑔−𝑢 𝑙|, where 𝑘2is an inertial permeability coefficient that quantifies how the geometry of the permeable network alters the gas velocity and pressure distributions. As a result, the linear Darcy interaction (46) between gas and liquid becomes quadratic at large drift velocities. This generalized Forchheimer relationship is dominated by the inertial effects when 𝑝0𝑤0𝑘0∕(𝜂𝑔0𝑐0𝑇0𝑘20)≫1, where 𝑘20 is a reference inertial permeability coefficient. The dimensionless factor of the 𝑅𝐻2𝑂terms in the momentum equations (16d)–(16e)and𝑝𝑙 equation (16i)is𝐹𝑟=𝑝 0𝑤0∕(𝜌𝑙0𝑢0𝑐0𝑇0). Figure 11G shows that 𝐹𝑟is well below 1 over the whole ascent, which is consistent with the fact that neglecting those terms does not affect much the mass discharge rate. It also suggests that 𝑤0and 𝑝0are the strongest controls of 𝐹𝑟. Assuming a𝐹𝑟threshold of 0.1, Figure 11H shows the maximum values of 𝑤0as a function of 𝑄0for characteristic values of 𝑝0. Neglecting the contribution of mass transfer in the momentum and 𝑝𝑙 equations is thus possible over the whole range of mass discharge rates for 𝑤0<0.01m/s at low pressure. At high pressure, where porosity is low and 𝑘∼10 −13–10−15 m2, Darcy’s law (46) suggests that 𝑤0∼10 −6–10−4 m/s. As a result, mass transfer terms can be neglected for 𝑄0>10 3kg/s at high pressure. Although a more complete parametric study is needed to explore the consequences of varying the parameters kept constant in this analysis (Table 3), we can propose a minimal set of equations by cumulating all the above simplifications. Omitting shear heating and the inertial and mass exchange terms in the liquid momentum equation, using Darcy’s law to calculate 𝑢𝑔,and ignoring mass exchange in the dilatancy equation yields a simplified system that can approximate the steady-state mass discharge rate of the full system within 10% if the conditions shown in Figure 11B,D,F,H are simultaneously fulfilled: 𝜕𝑡((1 − 𝜑)(1 − 𝐶𝑙))+𝜕 𝑧((1 − 𝜑)(1 − 𝐶𝑙)𝑢𝑙)=0, (51a) 𝜕𝑡((1 − 𝜑)𝜌𝑙)+𝜕 𝑧((1 − 𝜑)𝜌𝑙𝑢𝑙)=−𝑅 𝐻2𝑂,(51b) 𝜕𝑡(𝜑𝜌𝑔)+𝜕 𝑧(𝜑𝜌𝑔𝑢𝑔)=𝑅 𝐻2𝑂,(51c) 𝜕𝑡((1 − 𝜑)𝜌𝑙𝑢𝑙)+𝜕 𝑧((1 − 𝜑)𝑝𝑙)+ 𝛽𝑢𝑙−𝑝 𝑔𝜕𝑧(1−𝜑)−𝐾 𝑑𝜑(1 − 𝜑)(𝑢𝑔−𝑢 𝑙) = (1 − 𝜑)𝜌𝑙𝑔, (51d) 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 37 of 54 FIGURE 12 Temporal outputs from the transient model with 𝐷=10 −13 m2∕s and a single temperature. Times are also indicated in fraction of total ascent time of the liquid at steady state. The difference between the solution at 72 min and the stationary solution is not apparent at the figure scale. (A) Dissolved water content (𝐶𝑙). (B) Gas volume fraction (𝜑). 𝑢𝑔=𝑢 𝑙−1 (1 − 𝜑)𝐾𝑑 𝜕𝑧𝑝𝑔,(51e) 𝑐𝑝𝑚(𝜕𝑡𝑇+𝑢 𝑇𝜕𝑧𝑇)+𝑝 𝑔𝜕𝑧(𝜑𝑢𝑔+(1−𝜑)𝑢 𝑙)+(𝑐 𝑣𝑔 −𝑐 𝑝𝑙)𝑇𝑅𝐻2𝑂=− 4𝜂𝑙 3𝜑𝜌𝑙𝑅𝐻2𝑂𝜕𝑧𝑢𝑙,(51f) 𝑝𝑔=𝑐 0𝜌𝑔𝑇, (51g) 𝑝𝑙=𝑝 𝑔−4𝜂𝑙 3𝜑 𝜕𝑧𝑢𝑙.(51h) Importantly, this simplified system (51) no longer guarantees energy dissipation (Appendix B). The only term that could be neglected while ensuring energy dissipation is the shear heating term. 6.2 Transient behavior Reaching steady state from an initial condition with a porosity increasing linearly along the conduit length takes time because of disequilibrium degassing. Steady-state solutions feature small differences between 𝑝𝑔and 𝑝𝑙so that gas overpressure is not the limiting factor. Dissolved water, however, follows mass conservation and its distribution is controlled by a combination of advection by 𝑢𝑙and exsolution through 𝑅𝐻2𝑂. The initialization of 𝐶𝑙at 𝐶eq 𝑙can thus be far from the steady-state values, which can only be reached when 𝐶𝑙has been advected over the whole length of the conduit. This advection can easily be seen when diffusion is small, such as 𝐷=10 −13 m2∕s (Figure 12A). 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 38 of 54 BURGISSER et al. FIGURE 13 Gas volume fraction (𝜑) as a function of height in the conduit for different times of the transient model with a single temperature. Times are also indicated in fraction of total ascent time of the liquid at steady state. (A) 𝐷Ni.(B)𝐷=10 −8 m2∕s. The inlet speed in steady state for the 𝐷=10 −13 m2∕s run is ∼1 m/s, and thus 76 min are needed for a magma batch to cover the 5 km of the conduit length assuming constant speed. This calculation gives an estimate of how long it takes for steady state to be established. A more precise ascent time can be obtained by integrating the liquid travel time across each computational node, which yields an integrated ascent time of 72 min (absolute times are also reported as fraction of this integrated time). The stationary solution yields differences with the transient solution at 76 min of <10% on 𝜑and <5% on the other variables (𝐶𝑙,𝑝 𝑔,𝑝 𝑙,𝑢 𝑙,𝑢 𝑔,and𝑇). These differences between transient and steady states decrease to <2% when the transient solution reaches 150 min. The transient solutions at the transfer time of 𝐻∕(𝑢𝑙)|Γin thus offer a good representation of steady state. This is likely to be true under most conditions because 𝑢𝑙tends to increase toward the surface, making the transfer time based on the inlet velocity a minimum estimate. The effect of 𝐶𝑙advection on the porosity evolution is shown in Figure 12B. The transition between the constant value of 𝐶𝑙advected from the conduit base and the initial 𝐶𝑙distribution causes a kink in porosity. This kink migrates from 1 km at 33 min to 2 km at 50 min and to 3.5 km at 67 min. It is no longer apparent in the quasi-steady-state solution at 72 min because it has been advected past the shallowest part of the conduit. An interesting and persistent behavior occurs at higher 𝐷values. Figure 13 shows porosity losses at shallow depth and early times. This loss is due to a complex interplay between gas/liquid separation and exsolution. This low porosity region is transient with 𝐷Ni, disappearing after 17 min to yield a monotonous porosity increase over the whole conduit length that eventually reaches steady state after ∼6 h because the average ascent rate is ∼0.2 m/s (Figure 13A). The low porosity region is more marked at high 𝐷values and causes the system to be sensitive to initial conditions, sometimes hindering the progression from a given set of initial conditions to steady state 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 39 of 54 (Figure 13B). This temporary porosity drop is not due to numerical instabilities (it persists with a much finer grid, see Appendix A) and its origin lies in disequilibrium degassing. The behavior of the initial porosity loss is best understood with the transport equation on porosity (50): 𝜕𝑡𝜑+𝑢 𝑙𝜕𝑧𝜑= 3 4𝜂𝑙𝜑(1 − 𝜑)(𝑝𝑔−𝑝 𝑙) = (1 − 𝜑)𝜕𝑧𝑢𝑙+𝑅𝐻2𝑂 𝜌𝑙. The low porosity region corresponds to vanishing 𝑢𝑙, which suggests that the porosity loss can be understood by only looking at the right-hand-side terms. Unlike in the case of steady-state solutions (Figure 2), the term 𝑅𝐻2𝑂∕𝜌𝑙is the smallest of the three terms on the right-hand sides of (50). In other words, the behavior of porosity at shallow level shown in Figure 13 is due to gas underpressure instead of supersaturation like in the steady-state case. More precisely, what controls porosity reduction is the product 𝜑(1 − 𝜑)(𝑝𝑔−𝑝 𝑙), which is minimized when 𝜑∼0.5 and 𝑝𝑔<𝑝 𝑙. The gas underpressure is a consequence of gas escape toward the outlet because of the large permeability at moderate and high 𝜑, which causes 𝑢𝑔>𝑢 𝑙and, in turn, 𝑝𝑔<𝑝 𝑙.The porosity loss is thus a consequence of efficient liquid–gas separation near the vent. All steady-state solutions presented in Section 6.1 except that with 𝐷=10 −8 m2∕s were obtained by starting from an initial linear porosity profile reaching 50 vol% at the surface (i.e., 𝜑0=0.5, Figure 13A). The case with 𝐷=10 −8 m2∕s yields a porosity loss such that 𝜑reaches 0 at a depth of a few tens of meters (Figure 13B). This creates an impermeable plug where 𝑢𝑙→0. Exploring this transient regime is beyond the scope of this intermodel comparison focused on steady states, but, once 𝜑∼0at one computational cell around 360 s, further calculations to ∼366 s show that this region of vanishing porosity extends vertically, increasing plug thickness. This thickening occurs because the gas above the plug is being expelled toward the surface (i.e., 𝑢𝑔≥𝑢𝑙), which maintains gas underpressure almost all the way to the surface (it vanishes at the outlet according to the boundary condition (41a)). Changing the initial surface porosity coefficient, 𝜑0, within the range 0.1–0.7 did not preclude plug formation at 𝐷=10 −8 m2∕s, hereby excluding any steady state from these initial conditions. Instead, to reach steady state with that 𝐷value, the model was initialized with 𝜑0=0.25and 𝐷=10 −13 m2∕s for the first 5000 s. This allowed the porosity loss to be resorbed in a fashion similar to that depicted in Figure 13A. Then, 𝐷was changed to the target value of 10−8 m2∕s and the simulation was run until steady state was achieved. This procedure does not represent a natural scenario; its only purpose is to reach steady state at that 𝐷value. We applied the same initialization procedure to all the other cases with smaller 𝐷values (𝜑0=0.25and 𝐷=10 −13 m2∕s until 5000 s before setting 𝐷to the target value) and we obtained the same steady states as when setting 𝐷at the target value at 𝑡=0s. Also, varying 𝜑0did not change the respective steady states reached. The appearance of a low porosity region is not restricted to arbitrary initial conditions that are not compatible with the two-phase system (as in Figure 13). Magma permeability as deduced from eruptive products covers a broad range of values. In the permeability relationship (6b), we used the average exponent value of 3.4, but the full range of permeability from effusive products is covered by varying the exponent between 3 and 3.844. Figure 14 shows a transient run with 𝐷Ni that started from the steady-state solution of Figure 13A at 𝑡=0s. The exponent of the permeability equation (6b) was changed from the original, average value of 3.4 to the largest value of 3.8 in 1 s. This abrupt change can represent a sudden increase of permeability near the surface, such as one due to local shear zone development or fracturing. At 𝜑=0.6, this change causes permeability to 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 40 of 54 BURGISSER et al. FIGURE 14 Gas volume fraction (𝜑) as a function of height in the conduit for the transient model with 𝐷Ni and a single temperature. The condition at 𝑡=0s is the steady state and the later times shows the model response to a change of the permeability exponent in Equation (6b)from3.4to3.8. vary from 6×10 −11 to 10−11 m2. A shallow plug similar to those depicted in Figure 13 is created in ∼1 min. Thus, an abrupt but moderate increase in the shallow permeability of the magma can cause plug formation. 6.2.1 Alternate outlet boundary condition Similarly to the stationary case (Section 6.1.3), relaxing equal pressures at the outlet and setting instead (𝑝𝑔)|Γout =𝑝 top and (𝑝𝑙)|Γout according to (41b) yields almost identical transient, variable temperature results on all variables (e.g., Figures 12 and 13)but𝑝𝑙(cf. Figure 10). When this alternate boundary conditions is applied to the case of an overburden (𝑝top =2.5MPa), the change of the permeability exponent from 3.4 to 3.8 creates a plug directly centered on the outlet (i.e., 𝜑has only one local maximum at shallow depth instead of two like in Figure 14). This behavior, however, is not due to the boundary condition; the same 𝜑behavior occurs with equal outlet pressure and an overburden. This change in plug formation is instead caused by the fact that the overburden displaces the outlet position 100 m deeper than the case with 𝑝top =0.1MPa, effectively localizing the outlet inside the plug of Figure 14 where 𝜑→0. The average pressure condition (41c) is very sensitive to velocity changes near the outlet because both 𝑝𝑔and 𝑝𝑙depend on the liquid velocity divergence (gradient in 1.5D). As 𝑝𝑙is given explicitly in (16i) and does not appear in the other equations of the system (3) when using (7) (see the implemented system (A.3)), negative 𝑝𝑙values do not impede the numerical resolution. This is not true for 𝑝𝑔, which is present in most equations, including in the 𝑅𝐻2𝑂expression (6d) under 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 41 of 54 a square root. As a result, this condition is numerically unstable during the transient plug formation when 𝑝top =0.1MPa because the plug causes a cycle of acceleration–deceleration near the outlet that yields (𝑝𝑔)|Γout <0during the iterative process. This precluded us to obtain transient solutions similar to those of Figures 12 and 13. To obtain steady-state solutions with 𝑝top =0.1 MPa, we ran the model with 𝑝top =2.5MPa during the first 5000 s and then decreased 𝑝top down to the target value of 0.1 MPa following a smooth step function over 2000 s to finally reach steady state at ∼9000 s. The average pressure condition is thus unwieldy to study transient behavior with low outlet pressure and does not ensure strictly positive pressures. 7DISCUSSION We presented a 1.5D model of volcanic conduit flow with two phases, gas and liquid, that conserves water mass and ensures total energy dissipation. We compared the outputs of this transient model to those of the well-established steady-state model of Degruyter et al.12 Weassumedpuremelt, which is a restrictive assumption as most magmas contain crystals. If a two-phase framework is kept (i.e., solids are assumed to have the same velocity as the liquid), such solids mostly affect constitutive equations and have a limited effect on the structure of the transport equations presented herein. Crystals mostly change the rheology by inducing non-Newtonian behavior, and change magma permeability 𝑘. The reader is referred to Melnik and Sparks27 and Costa et al.63 for examples of how to take into account the effects of an additional solid phase on both the magma viscosity and density. Adopting a three-phase framework would be essential to predict solidinduced jamming and shear localization. A more immediate goal for future work than extending our model to three phases is to implement the 3D system (1) because it has the potential to address lateral distribution of pressure gradients, gas and heat loss at conduit walls, and to set a robust framework to study rheology changes due to shear localization. The isothermal comparison between our model and that of Degruyter et al.12 under conditions typical of an effusive eruption at an andesitic volcano suggests that strict equilibrium degassing is unreachable with a realistic water diffusion coefficient. Assuming equilibrium degassing recovers a realistic porosity evolution over most of the ascent, except near the surface where porosity is underestimated (Figure 1B). This is also the region where water balance is not respected (Figure 4). Overall, assuming equilibrium has a small influence on conduit flow pressure, ascent velocities, and thus on mass discharge rate. We have shown in Part I, Section 5 in Narbona-Reina et al.30 that the relaxed system with 𝐶𝑙=𝐶 eq 𝑙,𝑇𝑔=𝑇 𝑙,and𝑢𝑔=𝑢 𝑙can be expressed as three balances: conservation of the bulk mass, transport of 𝐶𝑇, and conservation of the bulk momentum with a dual pressure gradient. This system has a Rayleigh–Plesset dilatancy that is a function of ratios of energy densities accounting for the work done by changing phases to maintain chemical equilibrium. Our results show a match between the relaxed and nonrelaxed mass exchange rates of our transient model when 𝐶𝑙=𝐶 eq 𝑙, but no match with the isothermal rate of the Degruyter model. This implies that the energy balance must be taken into account to accurately reproduce equilibrium degassing. In other words, the assumption 𝐶𝑙=𝐶 eq 𝑙cannot be made in an isothermal system. Considering two pressures has two advantages. The first is that supersaturation pressures are readily available in the model outputs, which eases the assessment of whether a second nucleation event could be triggered (Figure 2C). Such events have been shown to occur during Vulcanian explosions,64 possibly triggering them65. The second advantage is that the transient behavior of 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 48 of 54 BURGISSER et al. 39. Dartevelle S, Valentine GA. Transient multiphase processes during the explosive eruption of basalt through a geothermal borehole (Námafjall, Iceland, 1977) and implications for natural volcanic flows. Earth Planet Sci Lett. 2007;262(3):363-384. 40. Collier L, Neuberg J. Incorporating seismic observations into 2D conduit flow modeling. J Volcanol Geotherm Res. 2006;152(3):331-346. 41. Collombet M. Two-dimensional gas loss for silicic magma flows: toward more realistic numerical models. Geophys J Int. 2009;177(1):309-318. doi:10.1111/j.1365-246X.2008.04086.x 42. Wong YQ, Segall P. Numerical analysis of time-dependent conduit magma flow in dome-forming eruptions with application to Mount St. Helens 2004-2008. J Geophys Res Solid Earth. 2019;124(11):11251-11273. doi:10. 1029/2019JB017585 43. González-Mellado AO, De la Cruz-Reyna S. A simplified equation of state for the density of silicate hydrous magmas: an application to the Popocatépetl buoyancy-driven dome growth process. J Volcanol Geotherm Res. 2008;171(3):287-300. 44. Mueller S, Melnik O, Spieler O, Scheu B, Dingwell DB. Permeability and degassing of dome lavas undergoing rapid decompression: an experimental determination. Bull Volcanol. 2005;67(6):526-538. doi:10.1007/s00445- 004-0392-4 45. Hess KU, Dingwell DB. Viscosities of hydrous leucogranitic melts: a non-Arrhenian model. Am Mineral. 1996;81:1297-1300. 46. Musser JM, Carney JE. Theoretical Review of the MFIX Fluid and Two-Fluid Models. DOE/NETL-2020/2100. NETL; 2020. 47. Gunn DJ. Transfer of heat or mass to particles in fixed and fluidised beds. Int J Heat Mass Transf. 1978;21(4):467- 476. 48. Russell JK, Hess KU, Dingwell DB. Models for viscosity of geological melts. Rev Mineral Geochem. 2022;87(1):841-885. 49. Giordano D, Russell JK, Dingwell DB. Viscosity of magmatic liquids: a model. Earth Planet Sci Lett. 2008;271:123-334. 50. Ni H, Zhang Y. H2O diffusion models in rhyolitic melt with new high pressure data. Chem Geol. 2008;250(1):68- 78. 51. Schlichting H, Gersten K. Boundary-Layer Theory.Springer;2017. 52. Kozono T, Koyaguchi T. Effects of gas escape and crystallization on the complexity of conduit flow dynamics during lava dome eruptions. J Geophys Res. 2012;117. doi:10.1029/2012JB009343 53. Rudy DH, Strikwerda JC. Boundary conditions for subsonic compressible Navier-Stokes calculations. Comput Fluids. 1981;9(3):327-338. 54. Poinsot TJ, Lele SK. Boundary conditions for direct simulations of compressible viscous flows. JComputPhys. 1992;101(1):104-129. 55. Choudhary A, Roy CJ, Luke EA, Veluri SP. Code verification of boundary conditions for compressible and incompressible computational fluid dynamics codes. Comput Fluids. 2016;126:153-169. 56. Hurwitz S, Navon O. Bubble nucleation in rhyolitic melts: experiments at high pressure, temperature, and water content. Earth Planet Sci Lett. 1994;122:267-280. 57. Gardner JE, Denis MH. Heterogeneous bubble nucleation on Fe-Ti oxide crystals in high-silica rhyolitic melts. Geochim Cosmochim Acta. 2004;68(17):3587-3597. 58. Lensky NG, Navon O, Lyakhovsky V. Bubble growth during decompression of magma: experimental and theoretical investigation. J Volcanol Geotherm Res. 2004;129(1):7-22. 59. Llewellin EW, Mader HM, Wilson SDR. The constitutive equation and flow dynamics of bubbly magmas. Geophys Res Lett. 2002;29(24):23-1-23-4. 60. Lyakhovsky V, Hurwitz S, Navon O. Bubble growth in rhyolitic melts: experimental and numerical investigation. Bull Volcanol. 1996;58(1):19-32. doi:10.1007/s004450050122 61. La Spina G, de’ Michieli Vitturi M, Clarke AB. Transient numerical model of magma ascent dynamics: application to the explosive eruptions at the Soufrière Hills Volcano. J Volcanol Geotherm Res. 2017;336:118-139. 62. Ruth D, Ma H. On the derivation of the Forchheimer equation by means of the averaging theorem. Transp Porous Media. 1992;7:255-264. 63. Costa A, Melnik O, Sparks RSJ. Controls of conduit geometry and wallrock elasticity on lava dome eruptions. Earth Planet Sci Lett. 2007;260:131-151. 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 49 of 54 64. Formenti Y, Druitt TH. Vesicle connectivity in pyroclasts and implications for the fluidisation of fountaincollapse pyroclastic flows, Montserrat (West Indies). Earth Planet Sci Lett. 2003;214(3-4):561-574. 65. Burgisser A, Arbaret L, Martel C, Forien M, Colombier M. The role of oxides in the shallow vesiculation of ascending magmas. J Volcanol Geotherm Res. 2020;406:107072. 66. Costa A, Melnik O, Vedeneeva E. Thermal effects during magma ascent in conduits. J Geophys Res Solid Earth. 2007;112(B12). doi:abs/10.1029/2007JB004985 How to cite this article: Burgisser A, Collombet M, Narbona-Reina G, Bresch D. Two-phase magma flow with phase exchange: Part II. 1.5D numerical simulations of a volcanic conduit. Stud Appl Math. 2024;153:e12747. https://doi.org/10.1111/sapm.12747 APPENDIX A: COMPUTATION DETAILS This appendix contains details of how the transient model is implemented in the finite-element commercial software COMSOL Multiphysics v5.6. For user-friendliness, we kept within the Appendix notations consistent with the software documentation. The spatial variable is the height along the conduit that we denoted by 𝑧in the main text. Nevertheless, for 1D problems in COMSOL Multiphysics, the independent variable is denoted by 𝑥. To avoid confusion in the computational code notation that is used in this appendix, we assume the notation 𝑥for the independent variable, keeping in mind that, when referring to the main text, 𝑧≡𝑥. Also, derivatives in this software are calculated as 𝑑(𝑓,𝑥) for any function 𝑓. We use the “Coefficient Form PDE” module of the software, which implements a general equation of the form: 𝑒𝑎𝜕2 𝑡2𝑈+𝑑 𝑎𝜕𝑡𝑈 + div(−𝑐∇𝑈 − 𝛼𝑈 + 𝛾) + 𝛽 ⋅ ∇𝑈 + 𝑎𝑈 = 𝑓, (A.1) where 𝑈is the vector of unknowns and 𝑒𝑎,𝑑 𝑎are the mass coefficients, 𝑐is the diffusion coefficient, 𝛼the conservative flux convection coefficient, 𝛾the conservative flux source, 𝛽the convection coefficient, 𝑎the absorption coefficient, and 𝑓the source term. The PDE module automatically converts the above equations written in coefficient form to the weak form, and that weak form is solved by the software using the finite-element method. The Dirichlet-type boundary conditions are introduced using the corresponding module “Dirichlet Boundary Condition” and the flux conditions through the “Flux/source” module that has the general expression: −𝑛 ⋅ (−𝑐∇𝑈 − 𝛼𝑈 + 𝛾) = 𝑔 − 𝑞𝑈, (A.2) where 𝑛is the outward vector normal to the boundary. Both boundary condition types are also converted into their weak form counterparts by COMSOL. A.1 Implementation of the transient model To write the computationalsystem, the nine unknowns of the differential system (16a), (16b), (16c), (16d), (16e), (16h), (16i), and the 1.5D version of (12)and(11) are denoted as follows: 𝑈=(Clphi,rholphi,rhogphi,𝑢 𝑙,𝑇 𝑙,𝑝 𝑔,𝑢 𝑔,𝑇 𝑔,𝑝 𝑙)𝑡, where Clphi =𝜌 𝑙(1 − 𝜑)(1 − 𝐶𝑙), rholphi =𝜌 𝑙(1 − 𝜑), rhogphi =𝜌 𝑔𝜑, 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 50 of 54 BURGISSER et al. are the conservative variables. The primitive variables are obtained as 𝜑=1−rholphi 𝜌𝑙;𝐶 𝑙=1− Clphi rholphi;𝜌 𝑔=rhogphi 𝜑. The momentum equations were codedin the default nonconservative form proposed in COMSOL, because it is more amenable toimplement boundary conditions. Most of the runs presented herein are performed with considering identical gas and liquid temperatures. We implemented an unified code whether 𝑇𝑔=𝑇 𝑙or not. For clarity, we write the system with the primitive variables and in the same order as in the vector of unknown 𝑈, 𝜕𝑡((1 − 𝜑)𝜌𝑙(1 − 𝐶𝑙))+𝜕 𝑥((1 − 𝜑)𝜌𝑙(1 − 𝐶𝑙)𝑢𝑙)=0, (A.3a) 𝜕𝑡((1 − 𝜑)𝜌𝑙)+𝜕 𝑥((1 − 𝜑)𝜌𝑙𝑢𝑙)=−𝑅 𝐻2𝑂,(A.3b) 𝜕𝑡(𝜑𝜌𝑔)+𝜕 𝑥(𝜑𝜌𝑔𝑢𝑔)=𝑅 𝐻2𝑂,(A.3c) (1 − 𝜑)𝜌𝑙(𝜕𝑡𝑢𝑙+𝑢 𝑙𝜕𝑥𝑢𝑙)+𝜕 𝑥((1 − 𝜑)𝑝𝑔−𝑅𝐻2𝑂 𝜌𝑙𝜑𝜒 −4(1 − 𝜑)𝜂𝑙 3𝜑 𝜕𝑥𝑢𝑙)+ 𝛽𝑢𝑙 −𝑝 𝑔𝜕𝑥(1−𝜑)−𝐾 𝑑𝜑(1 − 𝜑)(𝑢𝑔−𝑢 𝑙) = (1 − 𝜑)𝜌𝑙𝑔−𝑢𝑔−𝑢 𝑙 2𝑅𝐻2𝑂,(A.3d) 𝑐𝑝𝑙(1 − 𝜑)𝜌𝑙(𝜕𝑡𝑇𝑙+𝑢 𝑙𝜕𝑥𝑇𝑙)+(𝑐 𝑝𝑔 −𝑐 0)𝜑𝜌𝑔(𝜕𝑡𝑇𝑔+𝑢 𝑔𝜕𝑥𝑇𝑔)+𝑝 𝑔𝜕𝑥(𝜑𝑢𝑔+(1−𝜑)𝑢 𝑙) + ((𝑐𝑝𝑔 −𝑐 0)𝑇𝑔−𝑐 𝑝𝑙𝑇𝑙)𝑅𝐻2𝑂=−𝑅𝐻2𝑂 𝜌𝑙𝜑𝜒(𝜕𝑥𝑢𝑙+𝑅𝐻2𝑂 𝜌𝑙(1 − 𝜑))+ 𝛽𝑢2 𝑙,(A.3e) 𝑝𝑔=𝑐 0𝜌𝑔𝑇𝑔,(A.3f) 𝜑𝜌𝑔(𝜕𝑡𝑢𝑔+𝑢 𝑔𝜕𝑥𝑢𝑔)+𝜕 𝑥(𝜑𝑝𝑔)−𝑝 𝑔𝜕𝑥𝜑+𝐾 𝑑𝜑(1 − 𝜑)(𝑢𝑔−𝑢 𝑙)=𝜑𝜌 𝑔𝑔−𝑢𝑔−𝑢 𝑙 2𝑅𝐻2𝑂, (A.3g) 𝐵𝑇((𝑐𝑝𝑔 −𝑐 0)𝜑𝜌𝑔(𝜕𝑡𝑇𝑔+𝑢 𝑔𝜕𝑥𝑇𝑔)+𝑝 𝑔𝜕𝑥(𝜑𝑢𝑔+(1−𝜑)𝑢 𝑙))=𝛾(𝑇 𝑔−𝑇 𝑙) −𝐵 𝑇(𝑝𝑔 𝜌𝑙+(𝑐 𝑝𝑔 −𝑐 0)𝑇𝑔)𝑅𝐻2𝑂,(A.3h) 𝑝𝑙=𝑝 𝑔−1 𝜑𝜒(𝜕𝑥𝑢𝑙+𝑅𝐻2𝑂 𝜌𝑙(1 − 𝜑)).(A.3i) As (A.3e) is the 1.5D equivalent of (13) when 𝑇𝑔=𝑇 𝑙, we introduced the binary variable 𝐵𝑇to easily switch to a single temperature: 𝐵𝑇={1when 𝑇𝑔≠𝑇𝑙, 0when 𝑇𝑔=𝑇 𝑙. We also implemented a version of the code with 𝑇𝑔=𝑇 𝑙=𝑇and without (A.3h) to reduce the size of 𝑈and the associated resolution matrices in (A.1) for computational efficiency. 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 51 of 54 Coefficients for (A.1) are as follow: ∙𝑒𝑎=0, ∙𝑑𝑎(mass coefficient) is a 9×9 matrix with nonzero values (𝑑𝑎)11 =1,(𝑑 𝑎)22 = 1, (𝑑𝑎)33 =1,(𝑑 𝑎)44 =rholphi,(𝑑 𝑎)55 =𝑐 𝑝𝑙rholphi,(𝑑 𝑎)58 =(𝑐 𝑝𝑔 −𝑐 0)rhogphi,(𝑑 𝑎)77 = rhogphi,(𝑑 𝑎)88 =𝐵 𝑇(𝑐𝑝𝑔 −𝑐 0)rhogphi, ∙𝑐(diffusion coefficient) is a 9×9diagonal matrix with just one nonzero value 𝑐44 =4(1−𝜑)𝜂𝑙 3𝜑 ∙𝛼(conservative flux convection coefficient) is a 9×9 diagonal matrix with elements in the diagonal (−𝑢𝑙,−𝑢 𝑙,−𝑢 𝑔,0,0,0,0,0,0), ∙𝛾(conservative flux source) is a vector with nonzero values 𝛾4=(1−𝜑)𝑝 𝑔−𝑅𝐻2𝑂 𝜌𝑙𝜒𝜑 ,𝛾7=𝜑𝑝 𝑔 ∙𝛽(convection coefficient) is a 9×9 matrix with nonzero values 𝛽44 =rholphi 𝑢𝑙,𝛽55 = 𝑐𝑝𝑙rholphi 𝑢𝑙,𝛽58 =(𝑐 𝑝𝑔 −𝑐 0)rhogphi 𝑢𝑔,𝛽77 =rhogphi 𝑢𝑔,𝛽88 =𝐵 𝑇(𝑐𝑝𝑔 −𝑐 0)rhogphi 𝑢𝑔, ∙𝑎(absorption coefficient) is a 9×9matrix with the following elements: 𝑎44 =𝐾 𝑑𝜑(1 − 𝜑) +  𝛽, 𝑎46 = 𝑑(𝜑, 𝑥), 𝑎47 =−𝐾 𝑑𝜑(1 − 𝜑), 𝑎55 =−𝑐 𝑝𝑙𝑅𝐻2𝑂,𝑎 56 =𝑑(𝜑𝑢 𝑔+(1−𝜑)𝑢 𝑙,𝑥), 𝑎 58 =(𝑐 𝑝𝑔 −𝑐 0)𝑅𝐻2𝑂 𝑎66 =1, 𝑎74 =𝑎 47,𝑎 76 =−𝑑(𝜑,𝑥),𝑎 77 =−𝑎 47, 𝑎86 =𝐵 𝑇(𝑑(𝜑𝑢𝑔+(1−𝜑)𝑢 𝑙,𝑥)+ 𝑅𝐻2𝑂 𝜌𝑙),𝑎 88 =−𝛾𝐵 𝑇(𝑐𝑝𝑔 −𝑐 0)𝑅𝐻2𝑂, 𝑎96 =−1,𝑎 99 =1, ∙𝑓(source term) is a vector defined by 𝑓1=0, 𝑓2=−𝑅 𝐻2𝑂, 𝑓3=𝑅 𝐻2𝑂, 𝑓4=rholphi 𝑔−𝑢𝑔−𝑢𝑙 2𝑅𝐻2𝑂, 𝑓5=−𝑅𝐻2𝑂 𝜌𝑙𝜒𝜑 (𝑑(𝑢𝑙,𝑥)+ 𝑅𝐻2𝑂 rholphi )+  𝛽𝑢2 𝑙, 𝑓6=𝑐 0𝜌𝑔𝑇𝑔, 𝑓7=rhogphi 𝑔−𝑢𝑔−𝑢𝑙 2𝑅𝐻2𝑂, 𝑓8=𝛾𝑇 𝑙, 𝑓9=−1 𝜒𝜑 (𝑑(𝑢𝑙,𝑥)+ 𝑅𝐻2𝑂 rholphi ). Initial conditions are coded following (38) for the conservative variables, so (Clphi)|𝑡=0 =𝜌 𝑙(1 − 𝜑𝑖)(1 − 𝐶𝑖 𝑙); (rholphi)|𝑡=0 =𝜌 𝑙(1 − 𝜑𝑖); (rhogphi)|𝑡=0 =𝜑 𝑖𝑝𝑖 𝑔∕(𝑐0𝑇in); (𝑢𝑙)|𝑡=0 =(𝑢 𝑔)|𝑡=0 =𝑢 𝑖;(𝑝 𝑔)|𝑡=0 =(𝑝 𝑙)|𝑡=0 =𝑝 𝑖 𝑔; (𝑇𝑔)|𝑡=0 =(𝑇 𝑙)|𝑡=0 =𝑇 in.(A.4) The Dirichlet boundary conditions at the inlet follow the values defined in (40): (𝑝𝑙)|Γin =(𝑝 𝑔)|Γin =𝑝 in;(𝑇 𝑙)|Γin =(𝑇 𝑔)|Γin =𝑇 in;(rhogphi)|Γin =𝜑 in𝑝in∕(𝑐0𝑇in); (Clphi)|Γin =𝜌 𝑙(1 − 𝜑in)(1 − (𝐶𝑙)in); (rholphi)|Γin =𝜌 𝑙(1 − 𝜑in). (A.5) The Neumann condition (42) is coded using the “Flux/source” module through the general form equation (A.2). Notice that for equations with nonzero values of any of the coefficients 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 52 of 54 BURGISSER et al. appearing in this equation, 𝑐,𝛼,𝛾, it forces no-flux contribution at the inlet. In our system, this affects variables in positions 1,2,3,4,7. Removing this forcing is done by specifying the following coefficients: ∙𝑔, vector with components 𝑔1=−𝑢 𝑙Clphi 𝑛𝑥, 𝑔2=−𝑢 𝑙rholphi 𝑛𝑥, 𝑔3=−𝑢 𝑔rhogphi 𝑛𝑥, 𝑔4= −(1 − 𝜑)𝑝𝑔𝑛𝑥, 𝑔5=0, 𝑔6=0, 𝑔7=−𝜑𝑝 𝑔𝑛𝑥, 𝑔8=0, 𝑔9=0; ∙𝑞=0. In other words, all the nonzero terms in 𝑔1,2,3,7 simply make free conditions for the corresponding variables (0=0) due to the preprogrammed presence of conservative flux 𝛼,𝛾, and diffusion 𝑐. The boundary conditions at the outlet for pressures, (41a), are imposed as Dirichlet conditions: (𝑝𝑙)|Γout =(𝑝 𝑔)|Γout =𝑝 top. The additional numerical condition (42) is coded in the same way as the inlet boundary and it is combined with the alternative conditions for an hypothetical entering flow situation in (41d). It is imposed through the coefficients: ∙𝑔, vector with components 𝑔1=−𝑢 𝑙Clphi 𝑛𝑥−(Clphi −Clphi|Γout )(𝑢𝑙𝑛𝑥<0), 𝑔2=−𝑢 𝑙rholphi 𝑛𝑥−(rholphi −rholphi|Γout )(𝑢𝑙𝑛𝑥<0), 𝑔3=−𝑢 𝑔rhogphi 𝑛𝑥−(rhogphi −rhogphi|Γout )(𝑢𝑔𝑛𝑥<0), 𝑔4= −(1 − 𝜑)𝑝𝑔𝑛𝑥, 𝑔5=−(𝑇 𝑙−(𝑇 𝑙)|Γout )(𝑢𝑙𝑛𝑥<0), 𝑔6=0, 𝑔7=−𝜑𝑝 𝑔𝑛𝑥, 𝑔8=−(𝑇 𝑔−(𝑇 𝑔)|Γout )(𝑢𝑔𝑛𝑥<0), 𝑔9=0; ∙𝑞=0, where the variables with indices |Γout are set to their initial values evaluated at 𝑥=𝐻. 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License BURGISSER et al. 53 of 54 When the (𝑝𝑙)|Γout ≠(𝑝𝑔)|Γout condition is coded, the Dirichlet condition becomes: (𝑝𝑔)|Γout =𝑝 top (𝑝𝑙)|Γout =𝑝 top −1 𝜒𝜑(dtang(𝑢𝑙,𝑥)+(𝑅 𝐻2𝑂)|Γout ∕rholphi), where (𝑅𝐻2𝑂)|Γout =3 1∕3(4𝜋𝑁𝜑(1 − 𝜑))2∕3𝐷𝜌𝑙 𝐶𝑙−𝑘 ℎ√𝑝top (1 − 𝜑1∕3) and dtang(𝑢𝑙,𝑥)is an operator that calculates the derivative of 𝑢𝑙in the 𝑥-direction (for the 1D case) at the boundary where the evaluation takes place. This formulation takes advantage of the COMSOL setup that Dirichlet conditions are set in an unidirectional fashion, applying right-hand- side terms on the left-hand-side term (here (𝑝𝑙)|Γout ) but not on any variables appearing in the right-hand side of the condition. 𝐷is evaluated at 𝑝top if 𝐷Ni is used. The component 𝑔4above is replaced by 𝑔4=((1 − 𝜑)𝑝𝑔+𝑅 𝐻2𝑂∕(𝜌𝑙𝜒𝜑) + (1 − 𝜑)𝑑(𝑢𝑙, 𝑥)∕(𝜑𝜒))𝑛𝑥. A.2 Numerical method The sensitivity of solutions to domain discretization was tested at uniform intervals of 8, 4, and 1 m over the 5000-m long domain. Solutions were independent of interval sizes for the two smaller values, so 4-m intervals were chosen for all runs except the one featured in Figure 14, where 1-m intervals were used to better resolve the low-porosity region. We used cubic Hermite elements for the discretization of variables in the system, thus ensuring the continuity of the first derivatives at the nodes. The PARDISO solver is used with the default options, a relative tolerance on the residuals of 10−3 to 2×10 −3, and an absolute tolerance of 10−3. Variables were scaled according to their initial values. In transient runs, time stepping was done using the backward differentiation formula (BDF) algorithm with a maximum time step of 0.1 s. Stored time steps were those obtained from the solver to avoid interpolation inaccuracies. APPENDIX B: ENERGY OF THE SIMPLIFIED SYSTEM Following the same steps than in Section 2.5 of Narbona-Reina et al.30, we find formally (when no boundary conditions are applied) the following energy balance for the simplified system (51): 𝜕𝑡(1 2(1 − 𝜑)𝜌𝑙𝑢2 𝑙+ 𝑔𝑧(1 − 𝜑)𝜌𝑙+𝑐 𝑝𝑛𝑇)+𝜕 𝑧((1 − 𝜑)𝑝𝑙𝑢𝑙+ 𝑔𝑧(1 − 𝜑)𝜌𝑙𝑢𝑙+𝑐 𝑝𝑚𝑇𝑢𝑇) = −𝜑(1 − 𝜑)𝜒(𝑝𝑔−𝑝 𝑙)2−𝐾 𝑑𝜑(1 − 𝜑)(𝑢𝑔−𝑢 𝑙)2− 𝛽𝑢2 𝑙 −𝑔𝑧𝑅𝐻2𝑂+1 2𝑢2 𝑙(𝜕𝑧((1 − 𝜑)𝜌𝑙𝑢𝑙)+𝑅 𝐻2𝑂)−(𝑝 𝑔−𝑝 𝑙)𝑅𝐻2𝑂 𝜌𝑙.(B.1) The first two terms on the right-hand side are the dissipative contributions also present in the energy of the complete model. The term  𝛽𝑢2 𝑙stems from neglecting the shear heating in the temperature equation. The gravity term 𝑔𝑧𝑅𝐻2𝑂is a consequence of replacing the gas momentum equation by Darcy’s law. The term 1 2𝑢2 𝑙(𝜕𝑧((1 − 𝜑)𝜌𝑙𝑢𝑙)+𝑅 𝐻2𝑂)is related to both Darcy’s law and the neglected inertial term in the liquid momentum equation. Finally, the term 𝑅𝐻2𝑂(𝑝𝑔−𝑝 𝑙)∕𝜌𝑙 stems from neglecting mass exchange in the dilatancy equation (48). 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 54 of 54 BURGISSER et al. Is the simplified system dissipative? The three first terms on the right-hand side of (B.1)yield a dissipative total contribution. Note that the term |𝑢𝑙|2(𝜕𝑧((1 − 𝜑)𝜌𝑙𝑢𝑙)∕2 + 𝑅𝐻2𝑂)on the righthand side of (B.1) comes from the fact that in (51) we have neglected the inertial quantity 𝜕𝑧((1 − 𝜑)𝑢𝑙𝜌𝑙). If we also suppress the term 𝜕𝑡((1 − 𝜑)𝜌𝑙𝑢𝑙), then |𝑢𝑙|2(𝜕𝑧((1 − 𝜑)𝜌𝑙𝑢𝑙)∕2 + 𝑅𝐻2𝑂)does not appear in (B.1) but we lose also the quantity 𝜕𝑡((1 − 𝜑)𝜌𝑙𝑢2 𝑙)∕2 on the left-hand side of (B.1). If the time derivative on (1 − 𝜑)𝜌𝑙𝑢𝑙is also suppressed, then the time dependency of 𝑢𝑙is given through the nonstationary equations on 𝜑and 𝑇. This is what is called a semistationary system. Assuming such a semistationary system, the last quantity with an undefined sign is (𝑝𝑔−(𝑝 𝑙− 𝜌𝑙𝑔𝑧))𝑅𝐻2𝑂. Recall that 𝑅𝐻2𝑂has the same sign that 𝐶𝑙−𝑘 ℎ√𝑝𝑔. Therefore, we have the following cases: ∙(𝑝𝑔−(𝑝 𝑙−𝜌 𝑙𝑔𝑧))𝑅𝐻2𝑂≥0(dissipative case) if ◦𝐶𝑙≥𝑘ℎ√𝑝𝑔and 𝑝𝑔≥𝑝𝑙−𝜌 𝑙𝑔𝑧 ◦or 𝐶𝑙≤𝑘ℎ√𝑝𝑔and 𝑝𝑔≤𝑝𝑙−𝜌 𝑙𝑔𝑧 ∙(𝑝𝑔−(𝑝 𝑙−𝜌 𝑙𝑔𝑧))𝑅𝐻2𝑂≤0(nondissipative, unstable case) if ◦𝐶𝑙≥𝑘ℎ√𝑝𝑔and 𝑝𝑔≤𝑝𝑙−𝜌 𝑙𝑔𝑧 ◦or 𝐶𝑙≤𝑘ℎ√𝑝𝑔and 𝑝𝑔≥𝑝𝑙−𝜌 𝑙𝑔𝑧 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12747 by Readcube (Labtiva Inc.), Wiley Online Library on [14/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License