Bistability of buoyancy-driven exchange flows in vertical tubes

Buoyancy-driven exchange flows are common to a variety of natural and engineering systems, ranging from persistently active volcanoes to counterflows in oceanic straits. Laboratory experiments of exchange flows have been used as surrogates to elucidate the basic features of such flows. The resulting data have been analysed and interpreted mostly through core–annular flow solutions, the most common flow configuration at finite viscosity contrasts. These models have been successful in fitting experimental data, but less effective at explaining the variability observed in natural systems. In this paper, we demonstrate that some of the variability observed in laboratory experiments and natural systems is a consequence of the inherent bistability of core–annular flow. Using a core–annular solution to the classical problem of buoyancy-driven exchange flows in vertical tubes, we identify two mathematically valid solutions at steady state: a solution with fast flow in a thin core and a solution with relatively slow flow in a thick core. The theoretical existence of two solutions, however, does not necessarily imply that the system is bistable in the sense that flow switching may occur. Through direct numerical simulations, we confirm the hypothesis that core–annular flow in vertical tubes is inherently bistable. Our simulations suggest that the bistability of core–annular flow is linked to the boundary conditions of the domain, which implies that is not possible to predict the realized flow field from the material parameters of the fluids and the tube geometry alone. Our finding that buoyancy-driven exchange flows are inherently bistable systems is consistent with previous experimental data, but is in contrast to the underlying hypothesis of previous analytical models that the solution is unique and can be identified by maximizing the flux or extremizing the dissipation in the system. Our results have important implications for data interpretation by analytical models and may also have interesting ramifications for understanding volcanic degassing.


Introduction
Buoyancy-driven, bidirectional flow in channels or tubes is relevant to many natural and industrial processes. Examples include the joint pumping of viscous oil and water through pipelines (e.g., Joseph et al. 1997), the flow of cement into drilling mud in wellbores (e.g., Frigaard & Scherzer 1998), the counterflow of currents with different densities in oceanic straits (e.g., Dalziel 1992), and magma circulation in the conduits of persistently degassing volcanoes (e.g., Stevenson & Blake 1998). If the fluids are immiscible, the emerging flow is linearly unstable (e.g., Joseph et al. 1997) and significant † Email address for correspondence: jsuckale@stanford.edu arXiv:1802.01664v1 [physics.flu-dyn] 5 Feb 2018 effort has been devoted to identifying the pertinent flow regimes and their relative stabilities (e.g., Barnea 1987;Brauner 1998).
Here, we are interested in understanding the exchange flow of two Newtonian fluids in a vertical tube at low Reynolds number. Our work is motivated primarily by a need to better understand degassing processes in persistently active volcanoes, the most common form of volcanism on Earth. Many persistently degassing volcanoes have been active for periods exceeding the historical record. While erupting comparatively small amounts of lava, they continuously emit copious amounts of volatiles and thermal energy, suggesting that the majority of magma must be recycled back into the plumbing system after decoupling from the gas phase near the surface (e.g., Francis et al. 1993). At steady state, the ascent of gas-rich, buoyant magma in the conduit would therefore approximately balance the simultaneous descent of gas-poor, dense magma, which would result in exchange flow in the volcanic conduit (Kazahaya et al. 1994;Stevenson & Blake 1998;Burton et al. 2007;Witham 2011).
Our work builds on existing laboratory experiments of exchange flow in closed vertical tubes by Stevenson & Blake (1998). We focus on this specific set of analogue experiments, as their relatively simple geometry provides a suitable starting point towards understanding more complex geometries (Huppert & Hallworth 2007;Beckett et al. 2011;Palma et al. 2011). We analyze the flow behavior observed in the laboratory by using direct numerical simulations to virtually reproduce the original experiments. Our numerical model does not require closures such as drag forces, interface stresses or rise speeds (Suckale et al. 2010b,a;Qin & Suckale 2017). Instead, these flow properties emerge self-consistently from the computation, similar to analogue experiments where the flow dynamics emerge directly from the experimental setup and the materials used.
The most commonly observed flow regime of exchange flows in closed vertical tubes has a core-annular geometry, where one fluid flowing in the center of the tube (core) is surrounded by a film of the other fluid wetting the tube walls (annulus), as schematically illustrated in figure 1 (Bai et al. 1992;Petitjeans & Maxworthy 1996;Stevenson & Blake 1998;Scoffoni et al. 2001;Kuang et al. 2003). The result that core-annular flow appears to be the dominant flow regime for bidirectional flow in the presence of substantial viscosity contrasts may seem counter-intuitive. Core-annular flow is prone to instability, as has been shown through linearized stability analysis by numerous studies reviewed in Joseph et al. (1997). However, high viscosity contrasts and surface tensions may suppress the growth rate of interface instabilities sufficiently to stabilize core-annular flow (e.g., Hickox 1971). If the two fluids are miscible, even a small degree of mixing at the interface can further stabilize the flow configuration over a wide range of conditions (Meiburg et al. 2004).
The prominence of the core-annular regime in exchange flow is fortuitous, because coreannular flow is amenable to analytical analysis (e.g., Russell & Charles 1959;Kazahaya et al. 1994;Ullmann & Brauner 2004;Huppert & Hallworth 2007). Here, we build on Ullmann & Brauner (2004) to derive a simple analytical model of steady core-annular flow in a vertical tube. Our analysis generalizes the work of Kazahaya et al. (1994) by allowing for a mobile interface between the ascending and descending fluids. This generalization resolves the mismatch between theoretical predictions and laboratory observations reported by Stevenson & Blake (1998). Our analytical model improves upon the formulation of Huppert & Hallworth (2007) by eliminating the approximate representation of pressure and employing an alternative non-dimensionalization. The model predicts the existence of two distinct solutions at steady state and our numerical simulations confirm that both solutions can be realized in practice. This finding is consistent with laboratory experiments in more complex geometries (e.g., Beckett et al. 2011), which allow for flow in and out of an open vertical tube. We conclude that exchange flow in vertical tubes is a bistable system-an aspect of core-annular flow that could have important implications for understanding magma circulation in volcanic conduits.

Model Description
We combine two distinct and complementary model components: an analytical model of core-annular flow at steady state, and a direct numerical model of exchange flow in two dimensions (2D). An advantage of the analytical model is its applicability in both 2D and 3D, which is important for a detailed comparison of 3D laboratory experiments with numerical simulations in 2D. An advantage of the numerical model is the ability to capture the various flow regimes that arise for both the miscible fluids employed in the laboratory experiments and the immiscible fluids assumed in the analytical model. By combining the two approaches, we are able to leverage the flexibility of the numerical model in concert with the simplicity and explanatory power of the analytical model.

Derivation
At low Reynolds number, the fully-developed, steady-state core-annular flow of two immiscible and incompressible fluids in a pipe of inclination α from the vertical direction can be described by the 1D Navier-Stokes equations along the radial coordinate r ∈ [0, R] (see figure1), where the subscripts (·) d and (·) a denote the descending and the ascending fluids, respectively, R is the pipe radius and δ represents the unknown location of the interface between the ascending and descending fluids. Unlike in the core-annular model proposed by Huppert & Hallworth (2007) where the form of the pressure drop is postulated (Huppert & Hallworth 2007, Eq. (4.4)), dp/dz is an unknown constant to be determined here. The boundary conditions are no-slip at the tube wall, and continuity of velocity and shear stress across the fluid-fluid interface, Instead of a phase flow-rate scaling (Ullmann & Brauner 2004), we define the nondimensional quantitiesr where P = (dp/dz + ρ d g cos α)/(g∆ρ) is the non-dimensional pressure drop driving the descending phase, M= µ d /µ a is the viscosity ratio, and ∆ρ = ρ d − ρ a is the density difference. We integrate equations (2.7), while accounting for appropriately nondimensionalized boundary conditions, to obtain is the vertical flow speed at the interface given by The ascending flux in a closed tube with incompressible fluids must exactly balance the descending flux, where Te is the dimensionless flux or Transport number (e.g., Ullmann & Brauner 2004;Huppert & Hallworth 2007). The Transport number is related, but not identical, to the Poiseuille number, Ps, defined in Stevenson & Blake (1998), where u is the (dimensional) terminal rise speed. The Ps number therefore represents the non-dimensional rise speed, whereas the Te number captures the non-dimensional flux and hence Te = Ps πδ 2 . (2.12) We focus on Te in our analysis, becuase the flux is balanced between ascending and descending fluids and does not hinge on selecting either a characteristic rise speed from the spatially variable function, u a (r), or a core radius, δ. To fulfill the constraint of no-net flux in the tube, we substitute (2.8) into (2.10) and solve for the driving force P that depends only on the dimensionless parameters of the problem (i.e., δ, M, α), (2.13) We can now express Te as a function of P and δ Te = 2π P 16 (δ 2 − 1) 2 + δ 2 8 cos α(δ 2 − 1 − 2δ 2 log δ) . (2.14) We note that once P and δ are defined, Te is uniquely specified. The opposite, however, is not true since Te is a quadratic function of δ. For any Te below the maximum possible flux or flooding point, there are mathematically valid solutions. Despite differences in some of its technical ingredients, the model by Huppert & Hallworth (2007) also admits multiple solutions, but the authors implicitly impose uniqueness by maximizing the flux or, alternatively, extremizing the dissipation in the system. In §3, we numerically explore the validity of this uniqueness hypothesis and investigate whether both solutions may be realized.

Special Case: 2D Geometry
The analytical model we have now laid out for a general tube geometry in 3D can also be formulated in 2D. The latter is useful to compare our 2D numerical simulations to the 3D laboratory experiments. In the 2D case, the non-dimensional vertical speed in the two fluids is where 0 y 1 is the Cartesian direction perpendicular to the side boundaries and δ is the half-thickness of the ascending fluid. We will later use these expressions as boundary conditions to force bidirectional flow in our 2D numerical model (see §3.3).

Governing Equations
Our numerical model solves for conservation of mass and momentum. We assume that both fluids are incompressible. The governing equations are hence the incompressibility condition, (2.17) and the Navier-Stokes equation, where v is the velocity, p the pressure, g the gravitational acceleration, ρ the density, and µ the dynamic viscosity. We assume that both fluids abide by a Newtonian rheology. Density and viscosity vary in space and time to reflect the different properties of the two fluids. We consider the two contrasting scenarios of immiscible and miscible fluids. The immiscible case is useful for a direct comparison with the analytical solution, which is strictly valid only for two immiscible fluids. The miscible case represents the experimental setup by Stevenson & Blake (1998), which involves two miscible fluids with low diffusivity.

Immiscible fluids
In the immiscible case, a sharp interface separates the two fluids. We advect the curve, Γ , which represents this interface, with the flow field according to The interface deforms in response to the hydrodynamic forces acting on it. The material properties change discontinuously at the interface µ(x) = µ a in the ascending fluid µ d in the descending fluid , (2.20) core annulus  and may entail a jump in the pressure and normal stresses at the interface: Here, square brackets [·] denote a jump at the fluid-fluid interface, I is the identity tensor, T the viscous stress tensor, σ the surface tension, κ the curvature of Γ , n the unit normal vector on Γ pointing from the ascending towards the descending fluid, and t 1 and t 2 the two unit tangential vectors on Γ .

Miscible fluids
To allow for mixing, we introduce a continuous variable, c, which represents the concentration of the buoyant fluid. Initially, the concentration is set to unity in the buoyant fluid and zero in the heavy fluid. The concentration evolves over time due to advection and diffusion, where D is the diffusion coefficient. In the experiments, Stevenson & Blake (1998) used miscible fluids, such as syrup, dilute syrup, glycerol and water, but they did not provide the fluid diffusivities. In the absence of detailed measurements, we use a constant diffusivity of D = 10 −10 m/s 2 in both fluids for all simulations, a value motivated by the diffusivity measured for corn syrup in distilled water (Ray et al. 2007). We assume that density and viscosity depend linearly on the concentration c, such that (2.25) We note that additional complexity may arise in the vicinity of the interface if a nonlinear dependence of viscosity on concentration is assumed as is the case in some prior studies (e.g., Tan & Homsy 1986;Goyal & Meiburg 2006) (see online supplement for details).
To initialize the miscible simulations, we assume that the initial concentration field has an interface with a finite thickness of 1.5∆x, where ∆x is the coarse grid resolution (see §2.2.4 and figure 2). Hence, no discontinuous material contrasts arise in our miscible computations.

Numerical Methods
We discretize the governing equations, (2.17) and (2.18), using the numerical method derived, verified and validated in Qin & Suckale (2017). Building on Suckale et al. (2010b), and Suckale et al. (2012), we developed this computational technique specifically for multiphase flow problems with large viscosity contrasts at low Reynolds number. It consists of three main components. First, a multiphase Navier-Stokes solver that handles substantial and discontinuous differences in the properties of material phases by adopting an implicit implementation of the viscous term, time-step splitting, and an approximate factorization of the sparse coefficient matrices for computational efficiency. Second, a level-set-based interface solver that tracks the motion of a fluid-fluid interface through an iterative topology-preserving projection (Qin et al. 2015). The sharp interface solver is only pertinent if the two fluids are immiscible. In the computationally simpler case The aspect ratio of all laboratory tubes was 1:100 despite different physical lengths and widths. Here and after, only central part of numerical domain shown for better visibility.
of two miscible fluids, we solve an advection-diffusion equation for concentration. The third component is an adaptive grid refinement algorithm that increases resolution in the vicinity of the moving interface where most of the salient physical processes originate. Accurate interface advection is highly dependent on grid resolution (e.g., Sethian 1996), particularly in flows that hinge sensitively on the growth of interface instabilities. To maximize numerical resolution at the interface, we adopt an adaptive mesh refinement strategy that tracks the interface position over time. Figure 2 shows a close-up of the computational domain around the interface to schematically illustrate the grid refinement at the initial condition (figure 2A) and at a later time (figure 2), both shown at intentionally coarse resolution for easier visualization. The refined zone extends as the interface is stretched by flow. While the computational challenges associated with tracking diffusive interface are less pronounced, we adopt the same grid refinement strategy for the miscible case (figures 2C and D). For more details regarding the numerical technique, including the various benchmarks performed to verify and validate the numerical method, please refer to Qin & Suckale (2017) and the online supplement to this study. Table 1: Summary of analogue experimental data, and numerical and analytical results. Viscosity ratio (M ), density contrast (∆ρ), tube radius (R), tube length (L), front rise speed (u), Poiseuille number (Ps), dimensionless core radius (δ) and flow regime (Reg) observed in experiments (EXP) and numerical models (NUM). Thin-core, δ thin , and thick-core, δ thick , radii, and flow reversal flux ratio Terev/Te in analytical models (ANA).
Bold numbers indicate the solution observed in corresponding experiments and simulations.

Virtual Reproductions of Analogue Experiments
Stevenson & Blake (1998) initiate flow by inverting closed tubes filled with the denser, more viscous fluid (dark blue in our figures) in the lower half, and the buoyant, less viscous fluid (light blue in our figures) in the upper half. We initiate our models with the two fluids in an unstable density stratification and a slight cosine perturbation along their interface (see figure 2). We have verified that an initial condition of different symmetry yields qualitatively equivalent model outcomes, as shown in the online supplement. To represent the rigid glass walls of the experimental test tubes, we set all four side boundaries to noslip (v = 0). All simulations are computed on a 40×4000 grid with a factor 4 refinement applied to resolve the interface.
We numerically reproduce all eleven analogue experiments by Stevenson & Blake (1998), based on the detailed material properties they reported (see table 1). Figure  3 shows the resulting flows at times, where the front rise speed of the overturn flow has reached an approximate steady state; these are t = 200 × t 0 for experiments #1,2 and #5-9, and t = 1000 × t 0 for #3,4 and #10,11, where t 0 = R/U is the dimensionless, viscous time scale.
We generally observe the same behavior as reported by Stevenson & Blake (1998), which they classified into three different overturn styles. In figure 4, we illustrate these three overturn styles in the reproduced experiments #8, #9, and #10. At high viscosity ratios (M > 300) , the flow configuration is characterized by stable core-annular flow (flow regime I; figures 4A1-3). With decreasing viscosity contrast, interface waves become more pronounced. However, at intermediate viscosity ratios (10 < M < 300) the amplitude of interface waves remains small enough to avoid wave bridging and disintegration of the core (flow regime II; figures 4B1-3). Rather, the descending fluid intermittently rips off the tube wall, while the ascending fluid forms a coherent core in the center of the tube. Once the viscosity contrast becomes small (M < 10), wave bridging occurs and, as a consequence, both fluids sink or rise as separate batches in the center of the tube (flow regime III; figures 4C1-3). We note that all three flow regimes exhibit interface waves, but their amplitude and dynamic significance decrease as the viscosity contrast increases in agreement with theory (Hickox 1971). The transition between these three flow regimes is thus more gradual than this categorization into three distinct regimes may suggest. Figure 4 shows the results of miscible simulations (figures 4A2, B2, and C2) in comparison to three equivalent but immiscible cases (figures 4A3, B3, and C3). The three flow regimes appear qualitatively similar for both miscible and immiscible fluids. Since the thickness of the mixing boundary layer is small in comparison to the flow features, this finding is not particularly surprising. The miscible flows tends to have less variability along the interface, which confirms that even a small degree of miscibility notably dampens the growth rate of interface instabilities (e.g., Meiburg et al. 2004).
The availability of an analytical solution, at least in the limit of steady core-annular flow, provides a further opportunity to validate our numerical method beyond the benchmarks and tests presented in Qin & Suckale (2017). In figure 5, we compare the virtual reproduction of experiment #8 (figure 5A) for miscible fluids with its equivalent but immiscible case (figure 5B), and with the analytical solution calculated for the same flow parameters. To compare the numerical and analytical results, we take horizontal profiles of the vertical flow speed across the numerical domain in an area of well-developed core-annular flow. The flow profile of the immiscible simulation agrees remarkably well with the analytical result (figure 5C). The near fit between the immiscible numerical and the analytical solution is particularly encouraging considering that the latter was derived for steady state, while the former clearly arises in a transient flow.
The vertical speed profile, however, is significantly modified by fluid miscibility (figure 5C). The two reasons are a gradual change in the density and viscosity across the mixing layer and a more distributed shear stress in the interfacial zone. Despite the small for numerical simulations (green stars), analytical model predictions (black circles for 3D, blue squares for 2D), and the range of analogue experiments (grey shading) by Stevenson & Blake (1998). Analytical solution by Kazahaya et al. (1994) shown for comparison (black line).
amount of mixing here, the dynamic consequences are notable both in the increased maximum rise speed in the core, as well as in the narrower upwelling portion of the flow. We find that an exponential variation of viscosity with composition further magnifies this effect (see online supplement).

Analysis of Virtual Experiments
One of the main findings in Stevenson & Blake (1998) is that the Poiseuille number in their experiments is essentially constant (i.e., Ps ≈ 0.065) at finite viscosity contrast (approximately M > 10). This behavior is surprising given the five orders of magnitude variation in viscosity contrast tested in the experiments. It also stands in stark contrast to the theory of Kazahaya et al. (1994), which predicts a monotonic increase of Ps with M. Figure 6 compares the range of experimentally determined Ps values with the numerical values resulting from our simulations, and the theoretical values calculated from our analytical model in 2D and 3D. We also show the theoretical values of Kazahaya et al. (1994) for comparison. Since our analytical model does not allow us to quantify the rise speed of a transient interface as done in Stevenson & Blake (1998), we instead compute Ps from the average vertical flow speed in the core phase. For consistency, we use the same measure to calculate Ps from the numerical results as well. The online supplement provides a comparison of these two metrics for rise speed and shows that they produce comparable values in our simulations. Figure 6 shows that our 3D analytical model (black dots) agrees well with the range of observed Ps numbers (grey shaded) reported by Stevenson & Blake (1998) for M > 10. It also quantifies the difference between Ps in 2D (blue squares) versus 3D (black circles). In 3D, the rise speed is about twice as fast as in 2D. The 2D analytical estimates (blue squares) agree well with the 2D numerical estimates (green stars). All three pieces, the measured data, the analytical values and the numerical outcomes, indicate that Ps levels off with increasing M, in contrast to Kazahaya et al. (1994). Stevenson & Blake (1998) hypothesized that the disconnect between observed Ps values To quantify flow reversal more systematically, we compute the ratio between the flux in a given phase, Te, and the flow-reversal flux in the same phase Te rev , which we define as Te rev = δ δrev 2πru a (r) dr . (3.1) where δ rev is the point at which the vertical flow speed in the ascending fluid crosses zero, u a (δ rev ) = 0. With this definition, Te rev quantifies the flux of ascending fluid that is trapped in the flow reversal zone. Figure 8 shows the analytically predicted flow reversal flux as a function of dimensionless model parameters. Figure 8A shows its dependence on the core radius, δ, for different viscosity contrasts. For viscosity contrasts M > 10, the experimentally observed core radii cluster just below the values resulting in maximum flow reversal. In figure 8B, we plot the ratio of flow reversal to total flux, again as a   function of δ and M. We find that Te rev /Te is mostly insensitive to the viscosity ratio M for δ > 0.4. Thus, we find that for all experiments in flow regimes I or II, the flow reversal flux should constitute about 10% of the net flux. The insight that the interface between the two fluids is inherently dynamic does not yet explain why Ps remains approximately constant with increasing M. To find an explanation, we return to (2.14) in §2.1. It states that Te is a quadratic function of δ, which implies that two valid solutions for δ may exist for a given flux. The existence of multiple solutions is a common feature of multiphase flows, yet it does not necessarily imply that both solutions are realized in the laboratory (e.g., Landman 1991;Brauner 1998;Ullmann et al. 2003;Picchi & Poesio 2016), or indeed in natural systems.
In figure 9, we illustrate the two valid core-annular flow solutions for experiment #5 (see table 1). At the experimentally observed dimensionless flux, Te = 0.075, our model predicts two solutions for the core radius, the thick-core (δ thick , blue diamond) and the thin-core solution (δ thin , red triangle). The two corresponding flow profiles, shown from the center of the tube to the wall in figures 9B and C, highlight that the same overall flux can be accommodated either by a thin, rapidly ascending core with a thick annulus, or through a thick, slowly ascending core with a thin annulus. Interestingly, both solutions are far removed from the point of maximum flux or flooding point (yellow ×). Table 1 lists the two analytically computed core radii at the Te values inferred from the reported terminal rise speeds in all eleven experiments. The solutions compatible with their experimental and our numerical outcomes are printed in bold. We find that all experiments with viscosity contrasts of M > 10 select the thick-core solution, δ thick . The thin-core solution, δ thin , may be pertinent for the experiments with comparable viscosities (experiments #3, #4, #10 and #11), but these cases do not exhibit stable core-annular flow and the applicability of the analytical model is thus questionable. Figure 10 shows analytical flux solutions for all experiments as a function of core radius and viscosity contrast. For each Te-δ curve in figure 10A, we mark the solution that is realized in experiments (δ thick , red triangle; δ thin , blue diamond). We show curves for experiments #3, #4, #10 and #11 as light grey lines to convey that our analytical model is questionable for these cases. The experiments with M > 3 (solid lines) all cluster around δ ≈ 0.61 and Te ≈ 0.075. The scatter in Ps number in Stevenson & Blake (1998) (see figure 6) might hence not be entirely related to uncertainty of measurement, but reflect the fact that Te values at M > 3 are predicted to be very similar but not identical (see figure 10B). For experiments with a viscosity ratio close to unity (light grey lines), Te assumes much lower values (figure 10B), and the two possible solutions are no longer clearly distinct.
This analysis suggests that the dichotomy in Ps number observed by Stevenson & Blake (1998) and reproduced in figure 6 reflects a shift from thick-core core-annular flow (flow regimes I and II) to unstable overturn flow (flow regime III). Contrary to the thincore solution, the thick-core solution exhibits approximately constant core thickness over a large range of M (see figure 10A), which means that a constant Te translates to a constant Ps (2.12). If any of the experiments were exhibiting the thin-core solution, Ps would increase with M.
The result that only the thick-core solution, δ thick , is realized at finite viscosity contrast in Stevenson & Blake (1998) raises the question why this flow configuration appears to be preferable in this setting. This question is easiest to answer at high viscosity contrasts (i.e., M > 10, 000). As detailed in table 1, the predicted core radius for the thin-core solution, δ thin , becomes as small as a few percent of the tube radius at very high viscosity contrasts (e.g., experiments #2, #7 and #8). It is not surprising that a flow configuration with such a thin core is unstable. We recall that immiscible core-annular flow is always linearly unstable and only becomes metastable if the viscosity contrast is large enough to suppress interface waves to amplitudes that remain small compared to the core radius (e.g., Hickox 1971). At very high M, the thin-core solutions are hence highly prone to wave bridging and flow collapse (e.g., Barnea 1987).

Simulations of Persistent Exchange Flow
Determining the physical relevance of the two valid solutions based on core radius alone is unsatisfactory at intermediate viscosity contrast, when the two solutions, δ thin and δ thick , become increasingly comparable. In this section, we explore the possibility that the physically pertinent solution may be controlled not only by the material parameters but also by the domain boundaries, an idea raised but not further explored by Barnea & Taitel (1985). The experiments by Stevenson & Blake (1998) were performed in closed test tubes, which is a significant difference to natural systems, where exchange flow is typically the consequence of continued flux.
To generalize insights obtained from closed to open systems, we perform simulations with open boundaries at the top and the base of the model domain. All of the simulations discussed in this section are forced by a time-independent inflow condition imposed along the base of the domain and set according to the analytical speed profile in 2D (see (2.15)-2.16). We test both thin-and thick-core solutions by applying u| base (r) = u thin (r) or u| base (r) = u thick (r), respectively. We continue to enforce a no-slip condition (v| side = 0) along the side walls. For the outflow condition and the initial interface position we consider four different cases (see figure 11).
For the first case (figures 11A1-C), we impose a fixed outflow condition along the top boundary. We set the same analytical flow profile for the top and the base of the domain, i.e., u| base (r) = u| top (r) = u thick (r), or u| base (r) = u| top (r) = u thin (r), respectively. We also initiate the concentration field to the corresponding geometry extended through the whole domain (see figures 11A1 and B1). This choice implies that the interface is pinned both at the top and the bottom of the domain. While this setup is clearly contrived, it is an interesting end-member case for understanding the respective stability of the two solutions. Intuitively, one might expect that, if forced by the analytical solution on both ends, the flow field in the domain will remain close to that solution. Figure 11C shows that this is clearly not the case. While both interfaces are slightly wavy initially (figures 11A1 and B1), the thick-core solution stabilizes ( figure 11A3). However, the thin-core case begins transitioning to the thick-core flow field from the top boundary downwards almost immediately after the onset of flow (figures 11B2 and B3). The speed profiles taken in the middle of the domain after simulations have reached an approximately steady state confirm that the flow field for both simulations-the one driven by δ thick and the one driven by δ thin -closely approach the analytical thick-core solution ( figure 11C).
The second case (figures 11D1-F) represents a scenario with open outflow across a stress-free top boundary, which we enforce by setting p| top = const., and ∂v/∂z| top = 0. The flow field across the upper boundary is thus allowed to evolve over time, and the interface will move freely with the outflow. The initial position of the interface is identical to the first case (figures 11D1 and E1). The emerging flow fields correspond closely to the analytical solution imposed at the base of the domain. A thick-core inflow leads to stable thick-core configuration ( figure 11D3) and a thin-core inflow leads to stable thin-core configuration ( figure 11E3). This case demonstrates that the thin-core solution is a physically relevant flow configuration, at least at a viscosity contrast of M ≈ 30. Together with case one above, these simulations confirm that boundary conditions indeed have a profound influence on which mathematically valid flow configuration is realized in practice, as previously suggested in Barnea & Taitel (1985).
One advantage of the variable over the fixed outflow case is that we can now consider a different initial interface position. In our third case (figures 11G1-J), we initiate an interface confined to the vicinity of the lower boundary (figures 11G1 and H1), ensuring that the core radius corresponds to the analytical solution enforced at inflow. In these simulations, the interface moves through the domain in a similar way as in the experiments (figures 11G2 and H2), but eventually leaves the computational domain through the upper boundary. Immediately after onset, both the thick-and thin-core flows approach the thick-core solution ( figure 11G2 and H2). However, as soon as the interface intersects with the upper boundary, the flow in both simulations collapses back to the thin-core solution (figures 10G3 and H3). This case demonstrates that the preferable flow configuration does not only depend on the boundary conditions, but that a transient effect on one boundary, such as the passage of the interface through the top, can trigger a switch in the flow configuration propagating across the entire domain.
Since the analytical solution strictly only applies to the immiscible limit, we have so far only considered immiscible fluids in this section. However, we have shown above (figure 5) that miscible flows are comparable to immiscible ones if the thickness of the mixing layer remains small compared to the flow features. In the fourth case (figures 11K1-M), we therefore repeat case three with miscible fluids. We observe a similar, but less pronounced behavior. Initially, the thin core thickens immediately (figure 11L2), but both flows collapse back to a thin-core solution once the interface passes the upper boundary (figures 11K3 and L3). The steady-state vertical speed profiles are reminiscent of the thin-core solution but more spread out ( figure 11M). The interface also remains less wavy throughout, which supports the theoretical expectation that mixing along the interface stabilizes core-annular flow against interface waves (e.g., Meiburg et al. 2004).

Theoretical ramifications
Laboratory experiments have provided invaluable insights into the dynamics of buoyancy-driven exchange flows, but they inevitably simplify the more complex flow problem they intend to represent. One aspect in which many analogue models of exchange flow differ fundamentally from natural or industrial systems is that bidirectional flow only occurs as transient behavior until a steady state of complete inversion of the two fluids is reached (Stevenson & Blake 1998;Huppert & Hallworth 2007;Beckett et al. 2011;Palma et al. 2011). The steady state in the laboratory is hence very different from  Figure 11: Snapshots of numerical simulations forced with thick-and thin-core analytical solutions at the inlet. Material properties are identical to experiment #9 by Stevenson & Blake (1998). Results are shown for: fixed outlet boundary, enforcing the analytical model at the top boundary for thick-(A1-A3) and thin-core (B1-B3) solution; free outlet, stress-free top boundary allowing free flow through the top, starting from fully developed bidirectional configuration (D1-E3), and with a transient front moving through the domain (G1-H3); the same as previous case but with miscible fluids (K1-L3). Horizontal profiles of vertical speed (C, F, J, M) show numerical solutions approaching either thin-or thick-core analytical solutions depending on boundary conditions and fluid miscibility. the steady state in many natural or industrial problems of interest. Our analysis suggests that this difference may be consequential, because closed systems only realize a subset of the possible flow solutions that occur in open systems. Experiments in closed domains may hence systematically underestimate the dynamic variability of open-system flow.
We derive a simple analytical model that allows us to characterize the two steady-state solutions for core-annular flow. These two solutions refer to the same magnitude of flux, or Te, but differ in their core radii, δ thick and δ thin , flow speed profile, and degree of flow reversal. We also predict the flux curve, Te(δ), based only on the fluid properties and tube radius. The model does not require a fitting parameter to match experimental or numerical data, but does not predict the total flux nor the solution that is realized to achieve that flux. In fact, our analysis implies that it is not possible to predict these two outcome variables based on the fluid properties and the tube geometry alone, because the boundary conditions, above all pressure, play an important role in determining these. Huppert & Hallworth (2007) suggested that buoyancy-driven exchange flow in vertical pipe tends to maximize flux, which would be equivalent to maximizing Te. This argument implies that the flux should at the flooding point for most or all of the experiments (figure 10A), which is at odds with both the experimental and the numerical results discussed here. By inspection of table 1 and figure 10, we find that the flow system in this setup preferentially adjusts to the thick-core solution, despite the flux clearly not being maximized. Confirming Beckett et al. (2011), we conclude that viscous dissipation is not useful to distinguish the two possible flow configurations. In fact, we find that the non-dimensional viscous dissipation is directly proportional to the flux magnitude Te, and thus depends on δ in the same fashion (see online supplement).
Our finding, that the boundaries and flow history select the realized solution, is consistent with laboratory experiments that mimic open-system, buoyancy-driven exchange flow by connecting a vertical tube to fluid reservoirs at both ends (e.g., Huppert & Hallworth 2007;Beckett et al. 2011;Palma et al. 2011). The flow patterns that arise in this more complex geometry are more varied than in Stevenson & Blake (1998). We have reanalyzed the measured velocity profiles shown in figure 9 from Beckett et al. (2011) and find that their experiments #20, #15 and #11 (their figure 9 a,b, and c) correspond to the thin-core solution, experiment #9 (their figure 9 d) is close to maximum flux where only one solution exists, and experiments #17 and #16 (their figure 9 e and f) realize the thick-core solution. Our analysis of their data is shown in figure 12, where we give the Te(δ) curves and mark realized solutions for all experiments in their figure 9 (figure 12A). Our analytical model matches the measured speed profiles well for both the thin-and the thick-core example (figure 12B-C). The experimental data by Beckett et al. (2011) hence support our conclusion that both solutions are pertinent for understanding exchange flow in open systems.

Ramifications for Volcanic Systems
In its current form, our model is not yet suitable for a direct quantitative comparison with field data from a specific volcano. Some of our insights, however, may inform the interpretation of field data on a qualitative level. One pertinent insight is that a change at either the base of a conduit or its opening in a volcanic crater could potentially trigger a switch in the flow regime throughout the conduit. The effect that different flow regimes in the volcanic conduit have on eruptive surface processes is well explored and reviewed by Vergniolle & Mangan (2000). Our results suggest that the reverse is possible as well, i.e., eruptive surface processes can affect the flow regime in the volcanic conduit. For example, a disruption at the free surface, as might arise during an eruption or other events such as mass movements in the crater, could potentially trigger a switch in the flow configuration that is realized in the conduit. Of course, our simplistic boundary conditions (figures 11D1-E3 and G1-H3) do not adequately represent eruptive processes at a free surface. Nonetheless, our results demonstrate the potentially significant role of surface conditions in selecting a flow regime in the conduit.
We argue that the existence of two different, stable flow configurations could be reflected in erupted field samples from persistently degassing volcanoes, potentially from different stages of the same eruption. A switch from thick-to thin-core flow could increase the magma ascent rate by more than an order of magnitude, which may be detectable in microanalytical data. If a change in the estimated or measured ascent speed of the magma is detected, our results show that this observation does not necessarily require additional magma supply or a change in the volatile influx at depth. It could also be related to flow switching in the volcanic conduit.
Another potentially relevant insight of our study is the significant flow reversal present in all of the three flow regimes identified by Stevenson & Blake (1998). This result suggests that flow reversal might be the norm in volcanic conduits, rather than the exception. Our analytical model predicts that flow reversal occurs in the less viscous phase, a finding that is corroborated by our simulations (see figure 13 and additional results in the online supplement). Typically, the buoyant, volatile-rich magma has the lower viscosity, as was assumed in the experiment by Stevenson & Blake (1998). In that case, flow reversal flux is oriented downward, which raises the question whether the magma trapped in flow reversal is simply cycled back into a crustal reservoir never to be sampled by eruption, or whether some magma may circulate in the conduit for some time before being erupted. In the latter scenario, continued cycling along with mixing between different magma batches may lead to fundamentally different compositional evolution of both the magma and its volatile load than would be expected from a straight decompression path.
Flow reversal may also increase mixing between volatile-poor and volatile-rich magmatic melt, because the two fluids move in the same direction in some portion of the conduit. As demonstrated in figure 5B and C, even a small amount of mixing could have important dynamic ramifications for conduit flow. As pointed out by Witham (2011), magma mixing during conduit convection could be relevant for understanding why observed melt-inclusion trends from persistently degassing volcanoes rarely coincide with modeled degassing trends, as recently reviewed in Métrich & Wallace (2008).

Conclusions
In this study, we reproduce, explain and generalize laboratory experiments of buoyancy-driven exchange flow in vertical tubes. We derive an analytical model for coreannular flow-the most commonly observed configuration of bidirectional flow-that is consistent with laboratory observations and direct numerical simulations. Our primary purpose in this paper is to understand the flow behavior observed in the laboratory, but our model may also provide a suitable starting point for integrating additional complexity needed for quantifying conduit flow in persistently degassing volcanoes. The key finding from our analysis is that core-annular flow is bistable at finite viscosity contrast. This result implies that buoyancy-driven exchange flow is not uniquely determined by the material properties of the fluids and geometric parameters of the tube. The pressure and fluxes at the boundaries of the domain, along with the transient history of the flow, play an important role in selecting the realized flow solution.

S. Additional discussion of numerical model
This appendix provides: (1) An investigation on the impact of the initial condition on the bidirectional flow regimes; (2) A convergence test to assess the spatial and temporal resolution of the numerical model and a brief review of our adaptive grid refinement strategy; (3) A comparison of the average rise velocity used in the analytical model to the frontal rise speed used in analogue experiments; (4) A comparison of the linear dependence of viscosity on concentration to the nonlinear one; (5) The expression for calculating the viscous dissipation for concentric core-annular flows.
S.1. Impact of the numerical initial condition on the bidirectional flow We reproduce the analogue experiments with two kinds of perturbations applied to the interface between the two liquids, a sine and a cosine function. The flow regimes with three different viscosity contrast are shown in figure S1. At high viscosity contrast (experiment #8), steady core-annular flow forms from both initial perturbations. The descending annulus phase, however, is asymmetric when using a sine function initially (A2) in the sense that one side stretches further down than the other, whereas a cosine initial condition leads to a symmetric flow (A3). The degree of asymmetry depends on the magnitude of the initial perturbation. At intermediate viscosity contrast (experiment #9), the simulations with both sine (plot B2) and cosine (plot B3) initial perturbations produce a similar flow regime where the ascending liquid rises in the center of the tube and the descending liquid breaks into small blobs. At low viscosity contrast (experiment #10), the ascending and descending liquid is approximately symmetric for both initial conditions (figures S1C2-3). According to these results, the shape of the initial perturbation does not determine the flow regime, and will only have a minor effect on the flow symmetry. We also find that the simulations result in similar Ps number independent of their initial perturbation, even though the sine perturbation somewhat shortens the incubation time (see figure S1). Due to the negligible effect of the initial perturbation on the model outcomes, we only present the numerical simulations with a cosine initial condition in the manuscript.

S.2. Numerical convergence test and adaptive grid refinement
In figure S2, we reproduce experiment #5 to evaluate that the simulation converges as the spatial (figure S2A) and temporal (figure S2B) resolution are increased. We test convergence by comparing the numerical rise speed to the corresponding analytical solution. The numerical rise speed is taken along a horizontal cross section in an area of well developed bidirectional flow. Both convergence rates are 1st order. To increase the accuracy of the numerical model efficiently, we apply an adaptive grid refinement strategy. In the simulation of both immiscible and miscible flow, our grid hierarchy is composed of two levels. Level 1 is a fixed, coarse grid that defines the whole computational domain. Level 2 is an adaptively generated, fine grid that overlays the coarse level. The fine grid covers the whole liquid-liquid interface to accurately simulate the evolution of flow regimes. As the bidirectional flow develops and the interface extends through the domain, the area covered by the fine grid is increased as well. Compared to a simulation using the finer grid resolution over the entire computational domain, this adaptive grid refinement strategy obtains the same high accuracy, but saves more than 50% of the computational expense. More details of our adaptive grid refinement strategy are discussed in Qin & Suckale (2017).

I II III
sin cos sin cos sin cos Figure S1: Direct numerical simulations of the three primary flow regimes observed in bidirectional pipe flow for different viscosity contrast with sine initial perturbation (A2, B2, C2) and cosine initial perturbation (A3, B3, C3) in comparison to sketches (A1, B1, C1) reproduced from the analogue experiments by Stevenson & Blake (1998). The simulation snapshots are taken at 100s (flow regime I: A2 and A3), 300s (flow regime II: B2 and B3) and 200s (flow regime III: C2 and C3), respectively. [htb] (A) (B) Figure S2: Convergence tests for the relative LI error of numerical rise speed. Convergence of the relative error is 1st order for both spatial (A) and temporal (B) resolution.

S.3. Comparison of the average rise speed to the frontal rise speed
In their analogue experiments, Stevenson & Blake (1998) calculated the Ps number using the frontal rise speed of the ascending liquid. Our analytical model assumes that core-annular flow is at steady state, which is not compatible with a propagating front. To compute Ps with our analytical model we use the average rise speed, which we define as the average of the vertical velocity component in horizontal direction across the ascending liquid in the core, and use the experimental Ps numbers to estimate the core radius. Figure S3 compares the two velocity measures for experiment #8. Figure S3(B) shows that both frontal rise speed and average rise velocity reach a similar steady state after an initial spin-up. We also compare the two rise-speed estimates for other experiments and found that the difference is always small (<5%) and will decrease with spatial and temporal resolution.

S.4. Comparison of linear to nonlinear dependence of viscosity on concentration
Several previous numerical studies of the miscible fluids applied nonlinear dependence of viscosity on concentration instead of a linear one used in the present work. For example, Meiburg et al. (2004) defined an exponential relationship between viscosity and concentration as µ = µ a e ln µ d µa c . (S.1) In figure S4, we compare the simulation with a linear viscosity profile to the one with an exponential viscosity profile for experiment #5. In both simulations, we use a linear density profile. We find that the exponential viscosity profile is associated with a smoother flow regime (S4B) and a notably different speed profile (S4C) as compared to the linear case. The exponential viscosity profile enforces an increased downwelling velocity of the buoyant fluid. Since the density profile is linear, the buoyancy force changes linearly across the interface, whereas the viscous resistance changes exponentially ( Figure S4D). That means that there is a thin boundary layer of relatively heavy fluid that has relatively low viscosity and is therefore actively sinking. This actively sinking boundary layer creates the downward velocity peaks shown in figure S4C.
S.5. Viscous dissipation for fully developed laminar concentric core-annular flow We compute the dimensionless viscous dissipation for fully-developed, laminar, and concentric core-annular flow by integrating the local viscous dissipation rate over the where u d and u a are the velocity profiles of the descending and ascending phase respectively, see eqs. 2.18-19; the expression of Transport number Te is given in eq. 2.25. The viscous dissipation is proportional to the dimensionless flux and, in case of vertical core annular flow (cos α = 1), we get Φ = Te.
S.6. Interfacial speed for different viscosity contrasts The direction of motion of the interface determines the phase in which flow reversal occurs. As evident from out analytical model, the sign of u i depends sensitively on the driving force, which in turn depends on the viscosity contrast. Figure S5 shows the dependence of the dimensionless interfacial speed u i on the core radius, δ and the viscosity contrast, M. Negative interface speeds are associated with flow reversal in the ascending fluid and positive interface speeds lead to flow reversal in the descending fluid. When the viscosity of the descending phase is significantly higher than that of the ascending core (M > 10), flow reversal occurs in the ascending phase almost across the whole range of δ, except for extremely low core radii that are likely dynamically unstable. Conversely, as M decreases, flow reversal begins to shift into the descending phase over a growing range of core radii. If, however, the viscosity ratio is reversed from the experimental setup in Stevenson & Blake (1998), i.e., if the buoyant fluid is now more viscous than the descending fluid, flow reversal occurs in the descending phase across the entire range of δ.
[t]  Figure S5: The dimensionless interfacial velocity for the core-annular flow configuration as a function of the core-thickness δ and the viscosity ratio M: negative ui implies flow reversal in the ascending (core) phase, positive ui implies flow reversal in the descending phase.