# Characterization of Transport-Enhanced Phase Separation in Porous Media Using a Lattice-Boltzmann Method

Authors: Andrea Parmigiani; Paolo Roberto Di Palma; Sébastien Leclaire; Faraz Habib; Xiang-Zhao Kong
Journal: Geofluids
 (2019)

---

## Abstract

Phase separation of formation fluids in the subsurface introduces hydrodynamic perturbations which are critical for mass and energy transport of geofluids. Here, we present pore-scale lattice-Boltzmann simulations to investigate the hydrodynamical response of a porous system to the emergence of non-wetting droplets under background hydraulic gradients. A wide parameter space of capillary number and fluid saturation is explored to characterize the droplet evolution, the droplet size and shape distribution, and the capillary-clogging patterns. We find that clogging is favored by high capillary stress; nonetheless, clogging occurs at high non-wetting saturation (larger than 0.3), denoting the importance of convective transport on droplet growth and permeability. Moreover, droplets are more sheared at low capillary number; however, solid matrix plays a key role on droplet’s volume-to-surface ratio.


---

## Body

## 1. Introduction

Fluids flowing through natural reservoirs are often solutions of several chemical species/volatiles, some of which are less abundant and, depending on thermodynamic conditions, dissolve in the host fluids. Changes of reservoir thermodynamic conditions, though, can drive the solution out of equilibrium such that the host fluids might become oversaturated with one or more chemical species [1]. Under some circumstances, this oversaturation eventually leads to volatile exsolution and/or phase separation with the emergence and growth of a second phase. These circumstances have been reported in (1) gas bubble migration in sediments and subsurfaces [2, 3], (2) CO2 sequestration [4–10] where the dissolved CO2 exsolves from solution due to depressurization and/or increase in temperature, (3) hydrothermal phase separation [11–13] where a denser brine and a lighter low-salinity vapor coexist at elevated temperatures and pressures, and (4) gas exsolutions in magma reservoirs which can be induced by depressurization episodes and/or crystallization of anhydrous mineral phases driven by cooling [14–16]. This conversion from a single-phase to a multiphase system, particularly a mixture of immiscible fluids, is known to alter the system dynamics due the emergence of, such as buoyancy, viscosity contrast, and capillarity. In particular, the capillary resistivity to the multiphase transport influences the discharges of phases by means of the competition between viscous and capillary stresses [5, 7, 10, 17, 18]. This conversion in the mid-ocean ridge hydrothermal systems is also thought to strongly influence the chemistry of vents and the convective transfer of energy [11, 19]. In shallow ocean sediments, this conversion leads to the formation of gassy sediments, due to the arising and growing of gas bubbles. These bubbles affect the mechanical properties of the sediment [20, 21] and could dramatically reduce the sediment hydraulic conductivity [22].While the aforementioned studies have elucidated some effects of phase separation on some convective systems, these studies have not focused on dynamic response of phase separation in a heterogeneous porous medium under a background pressure gradient [5, 7–9, 16]. In this paper, we investigate, via multiphase numerical modeling, the system’s hydrodynamic response to phase exsolution/separation out of a solution, where fluid flow is driven by a constant hydraulic gradient and where the buoyancy stress potentially rising from the emergence of a dispersed phase can be neglected (i.e., low Bond number). We focus our discussions on two main aspects: (1) the characterization of the evolution of the separating/dispersed non-wetting phase (including size, shape, and mobility of the emergent droplets) and (2) phase separation effects on fluid discharge (i.e., hydraulic response to the evolving capillary stress on the emergence of two-phase fluid flows). In the interest of our analysis, density and viscosity of both phases are assumed to be identical (i.e., both the capillary and viscous stresses dominate over the buoyancy stress), and the theme of this work is mainly on the effects of variations in capillary properties (surface tension) and fluid saturations on the system evolution.


## 2. Conceptual Model and Simulation Set-Up

Our numerical model is designed to investigate the hydrodynamic response of the system during the emergence, growth, and transport of non-wetting droplets in a complex porous geometry where the flow of the wetting phase is driven by a constant hydraulic gradient. The model manipulates the emergence of droplets using a diffusion process and resolves the competition between the capillary and viscous stresses at the pore scale. For convenience, we avoid modeling the complex thermodynamics of droplet nucleation by taking the out-of-equilibrium conceptual model, i.e., the so-called spinodal decomposition [23–25]. The spinodal decomposition describes a mechanism that governs a rapid demixing of a multiphase mixture into multiple coexisting/immiscible fluids. This mechanism obeys a purely diffusion process and triggers phase separation where droplets are formed without overcoming an energetic barrier. Once the formation of droplets is initiated, the growth of droplets is driven by the droplet-droplet hydrodynamic interactions and coalescence. Eventually, the system evolves into a two-phase system with a continuous wetting phase and a dispersed non-wetting phase (i.e., droplets). This simplification allows us to introduce in situ, random droplet nucleations and to focus on the physical processes of droplet growth during the simulations of two immiscible fluid flows under isothermal and isobaric conditions.In this work, a variant of the lattice-Boltzmann two-phase color-gradient method (CGM) described by Leclaire et al. [26] is employed to perform the simulations of pore-scale isothermal immiscible two-phase fluid flow in porous media. The lattice-Boltzmann method (LBM) is a finite difference computational fluid dynamic approach that, in contrary to the traditional discrete solvers of Navier-Stokes equations, solves a simplified discrete version of the Boltzmann equation [27]. The LBM excels in the modeling of fluid-fluid and fluid-solid interactions of multiphase fluids within structures of complex geometry such as porous media [28, 29]. The CGM has been shown to be capable of correctly simulating a time-scale power law which is related to the coalescence rate of spinodal decomposition [24]. In addition, the CGM has been extensively validated against various benchmarks [26, 30–32]. The CGM has also been reported to be well suited to resolve the competition between capillary and viscous stresses as occurring in porous media [26]. In order to take advantage of high parallelization efficiency, our CGM code is implemented in PALABOS [33], an open-source C++ parallel library for LBM. All reported simulations in the current work are performed using Euler, the high-performance scientific compute clusters, at ETH Zurich.For this study, the CGM performs time integrations of two ensembles of particle populations (e.g., red and blue) that describe a continuous wetting phase and a dispersed non-wetting phase, respectively. By employing a three-dimensional D3Q15 lattice, the time integration of the particle population collision-streaming scheme is divided into seven calculation steps:(1)
First, the particle populations at the computational domain boundaries are updated according to the classical fully periodic boundary conditions [34]
(2)
The classical weakly compressible Boltzmann fluid hydrodynamic is introduced into the color-gradient method with a color-blind collision between the particles. The hydraulic gradient is modeled by including an external force on the wetting phase [32, 35]
(3)
The wetting boundary condition is applied (i.e., the dispersed phase is perfectly non-wetting with a 180-degree contact angle) using ghost nodes as described by Leclaire et al. [31]
(4)
A perturbation operator is introduced to model the interfacial tension at the interface between the two immiscible phases
(5)
The finite width of the interface and the immiscibility is preserved by using an additional recoloring step
(6)
A no-slip boundary condition between the colored particles and the solid phase is applied using a classical full-way bounce-back [27]
(7)
Finally, the usual streaming step is applied to each of the colored particle populationsCompared to the CGM described by Leclaire et al. [26], the current model applies two simplifications so that the total computational expense of this study can be reduced. The first simplification is the implementation of a single-relaxation-time model which is a special case of the multirelaxation-time model with χ=1, as presented by Leclaire et al. [26]. The second simplification is the employment of the ghost node approach, instead of the advanced wetting boundary condition described by Leclaire et al. [26]. It has been demonstrated that the ghost node approach yields satisfactory wettability behaviors (i.e., a correct contact angle) for a dispersed non-wetting phase with a contact angle larger than 120 degrees (see the supplementary material in Leclaire et al. [31]). We would like to point out, however, that if the dispersed phase is more wetting than the continuous phase, i.e., the contact angle of the dispersed phase is less than 90 degrees, a much more careful attention on the implementation of the wetting boundary condition would be necessary, because nonphysical mass transport could appear along the solid surfaces [26, 31, 36]. In order to concise the text and emphasize the focus of this study, the detailed mathematical description of the proposed CGM model (i.e. steps 2, 4, 5, and 7) is documented in the appendix. The mathematical description of the boundary conditions (i.e., steps 1, 3, and 6) refers to the aforementioned references.The current CGM belongs to a so-called diffused-interface class of immiscible multiphase solver where 3-5 lattice nodes are required to resolve the fluid-fluid interface between the immiscible fluids. However, the low-dissolution properties peculiar to the recoloring method (RM) during the recoloring step allow the formation of stable small droplets (radii smaller than 5 lattice nodes) which, although underresolved, can be advectively transported by the fluid flow [29, 37]. These small droplets might coalesce with each other or with larger droplets and, therefore, fully participate in the evolving dynamics that we want to investigate. More importantly, the RM instinctively guarantees the mass conservation for both the wetting and the non-wetting fluids. In our calculations, a flag variable, Ψ=ρw−ρnw/ρw+ρnw, is defined using the density fields of the wetting fluid, ρw, and the non-wetting fluid, ρnw. Regions/lattice nodes occupied by the wetting and the non-wetting fluids are indicated by Ψ<0 and Ψ≥0, respectively (e.g., droplets in Figure 1(c) are the volume rendering of Ψ≥0). The flag Ψ is also used to retrieve the total momentum field, j, and the individual momentum fields, jw and jnw, of the wetting and non-wetting fluids, respectively. Integration of the vertical components of the momentum fields, namely, jz, jzw, and jznw, over the whole simulation domain yields the total discharge, Jztot, and its individual contributions, Jzw and Jznw, of the two immiscible fluids, respectively (Figure 1(a)).Figure 1
(a) Three-stage (i.e., stages I, II, and III) temporal evolution of total discharge,JZtot, and non-wetting fluid discharge, JZnw, for three σlb values and Snw=0.05. (b) Temporal evolution of Snw distribution (simulation time is normalized to the maximum lattice time of 6×105) as a function of porous medium conductivity for null (middle panel) and non-null (bottom panel) hydraulic gradients with σlb=0.05 and Snw=0.05. The top panel shows the histogram of the normalized velocity. (c) Droplet evolution at two different times (early and later times at stage II) during spinodal decomposition for σlb=0.05, Snw=0.05, and Flbp=2.0×10−4.
(a)
(b)
(c)The porous medium used in the current simulations is generated using a crystal nucleation-and-growth algorithm [38] that follows Hersum and Marsh [39]. The computational domain measures in total nz×ny×nx=300×200×200 lattice nodes. The porosity of this porous medium, ϕ, is set to 0.5 to ensure a high enough numerical resolution for resolving a significant size range in both pores and droplets (the characteristic pore size is about 10-20 lattice nodes). This high porosity thus preserves a statistically meaningful number of crystals and droplets (each about thousands, establishing critically necessary large representative elementary volumes (REVs)). It therefore allows us to draw conclusions on fluid discharges for a wide range of phase saturation, S.A set of simulations with a non-wetting phase saturation of0≤Snw≤0.7 and a surface tension of 5×10−4≤σ≤10−1 (in lattice units) has been performed to focus on the effect of variation in capillary properties (surface tension) and fluid volume fractions/saturations. For an individual simulation, its fluid saturation S of individual calculation is preset and remains constant during the simulation. Except for some small random perturbations (which are needed to trigger the phase separation process), the two fluids are initially uniformly mixed and distributed throughout the void space of the porous medium. The viscosity of both phases is set to vlb=0.1667 (in lattice units), which corresponds to a relaxation time of τ=1.0 for all the simulations. Periodic boundary conditions are applied in all directions. In all simulations, a constant hydraulic gradient is implemented by applying a body force of Flbp=2×10−4 on the wetting phase only [40]. Therefore, the total driving force on the system is proportional to the wetting phase volume, Vw. The dispersed phase is assumed to be completely non-wetting (i.e., its contact angle is equal to 180°). Under such conditions, the transport of the non-wetting droplets depends solely on the competition between the shear (viscous) and capillary stresses.


## 3. Results and Discussions

### 3.1. Transport Stage and Mechanism

A typical three-stage evolution of the phase separation is clearly demonstrated in our simulations (Figure1): transport of mixture in stage I, transport of growing droplets in stage II, and steady-state transport of two phases in stage III. During stage I, the two-fluid mixture quickly establishes an early/initial steady-state discharge, jztot, which depends on Flbp and the wetting fluid volume fraction Sw (Figure 1(a)). At this stage, the two phases are still completely mixed (i.e., droplets are not yet formed). This complete mixing is indicated by the zero discharge, jznw, of the non-wetting phase. At the beginning of stage II, phase separation is triggered and droplets start to develop via the spinodal decomposition, dominated by a diffusion growth (Figure 1(a)). A few thousand iterations later, droplets hydrodynamically interact with each other, and they grow via coalescence. As soon as the droplets reach sizes similar to the ones of pores and throats, capillary stresses start affecting the fluid flow. These observations are supported by the confirmation of the increase in jznw (due to the growth and transport of droplets) and the decrease in total discharge, jZtot, induced by the capacity reduction of flow pathways which results in higher energy dissipation. At late times of stage II, droplets can reach sizes similar to or even bigger than the ones of pores and throats. A typical example of droplet evolution is reported by taking snapshots at early and late times during stage II in one of our simulations (Figure 1(c)). A more comprehensive understanding of the phase separation process at early times can be gained by watching movies included in our supplementary materials (available here). In stage III, both the total and non-wetting discharges eventually reach their steady-state values. During this stage, transport of large droplets is strongly controlled by the competition between capillary and shear stresses. At high surface tension (indicated by the red line in Figure 1(a)), large droplets are mostly trapped and they act as immobile obstacles to the flow of the wetting phase. At small surface tensions (indicated by the blue and green curves in Figure 1(a)), droplets are likely to be deformed and they can migrate under the driving of the wetting phase. The latter case results in much higher total and non-wetting discharges. Nonetheless, both the total and non-wetting discharges fluctuate around their steady-state values due to the capture and release of droplets. Apparently, this fluctuation is lessened due to snap-off and reconnection of droplets as surface tension decreases, which will be further explained in the subsequent subsections.Since phase separation is occurring while the fluid mixture is being transmitted, droplet growth is expected to be strongly influenced by the connectivity of the porous medium and the advection of fluids. This observation is confirmed in our calculations that growth of droplets is governed by a transport-enhanced mechanism, as described by Zuo et al. [10]. Such behavior is clearly visible (Figure 1(b)) when comparing the evolution of saturation distribution Snw of the non-wetting phase between two particular scenarios, one with a null hydraulic gradient (i.e., droplet growth is diffusion-dominated) and one with a non-null hydraulic gradient (i.e., droplet growth is influenced by fluid advection). To interpret such behavior, for both scenarios, the saturation of the non-wetting phase, Snw, is assembled as a function of the normalized z-component of fluid velocity. The fluid velocity is obtained in a steady-state single-phase fluid flow calculation with the same geometry and hydraulic gradient. It is rational to approximate that the strength/capacity of advection is proportional to velocity magnitude, i.e., the higher the z-velocity, the higher the advection capacity. For the scenario of null hydraulic gradient (i.e., diffusion-controlled and affected at a lesser extent by geometrical properties), Snw quickly becomes a roughly uniform distribution, where non-wetting droplets mainly reside in regions of median advection capacity (Figure 1(b)). For the scenario of non-null hydraulic gradient, a much more rich temporal evolution of Snw is registered, especially at the early times when small droplets are transported in the highly advective regions (i.e., the most well-connected pores and throats) but eventually coalesce to form large droplets which reside in regions of low (but non-null) advection capacity and likely to be capillary-trapped.


### 3.2. Connectivity of Wetting Fluid

As stated in previous sections, a steady-state velocity field is reached during stage III. This steady-state velocityuw of the wetting phase is utilized to compute the continuous streamlines of the wetting phase only. Here, the continuous streamline is defined as a streamline continuously spanning from the inlet to outlet boundaries along the main flow direction. By definition, in the lattice-Boltzmann method, the velocity uwx,y,z is calculated as the ratio of momentum over fluid density, i.e., uw=jw/ρw. Then, a seepage velocity of the wetting phase in the main flow direction (z-direction) U¯zw is calculated along the streamlines. Note that streamlines are generally tortuously passing through the porous medium and thus have a length longer than the 22domain size in the z-direction [41]. Mathematically, the z-direction seepage velocity of the wetting phase is formulated as
(1)U¯zw=1N∑s=1N∑i=1Lsuzwxi,yi,zi,where Ls is the length of the s-th continuous streamline and N is the total number of the continuous streamlines that can be identified and pass through the porous medium from the initial cross section z=0 in an individual simulation. A total velocity, U¯w, along all the continuous streamlines is calculated accordingly:
(2)U¯w=1N∑s=1N∑i=1Lsuxwxi,yi,zi2+uywxi,yi,zi2+uzwxi,yi,zi2.Starting from the initial cross sectionz=0, velocities along the streamlines are sampled at the equidistant points, using an integration step equals to 0.1 of a lattice to ensure a good estimation of the continuous streamlines passing through the pore space. The presence of a notable amount of continuous streamlines indicates that (1) continuous pathways still exist in the wetting phase (i.e., it is physically meaningful to use body force to formulate a background hydraulic gradient) and (2) the porous medium is not yet completely clogged by the non-wetting phase. A completely clogged stage is reached when the wetting phase retains very limited number of continuous streamlines (up to 100). Note that there are more than 8000 continuous wetting-fluid streamlines for a single-phase flow in the current porous medium. At the completely clogged stage, the trapped or slowly moving non-wetting fluid breaks up the wetting phase into several disconnected parts inside the porous medium and thus impedes the transmit of the wetting phase. Figure 2(b) shows the completely clogged stages for five different surface tensions σlb, indicated by arrows with different colors. Our results clearly indicate that the completely clogged stage occurs at higher wetting-fluid saturation Sw for higher σlb, due to higher resistance introduced by capillary stresses.Figure 2
(a) Magnitude of the norm of total velocity throughout the porous medium (solid in back and gray) as well as non-wetting droplets (in red) atSnw=0.05 and σlb=0.05 at a given time snapshot (3×105 iterations). (b) Mean seepage velocity, U¯zw, averaged along the continuous streamlines of the wetting phase. (c) and (d) Magnitude of the norm of seepage velocity (solids are rendered transparently in gray) and non-wetting droplets (in red) for Snw=0.05 (left) and Snw=0.20 (right) with σlb=0.05 at 3×105 iterations. For better visualization, not all streamlines are reported.
(a)
(b)
(c)
(d)Since the viscosity of both phases is fixed, the competition between capillary and viscous stresses is therefore explored by uniquely varyingσlb. This competition can be cast into a dimensionless capillary number Ca=ucwvlb/σlb, where ucw is the characteristic velocity of the wetting phase. If ucw is defined as the mean seepage velocity, U¯zw, of the wetting phase, i.e., before the system is completely clogged, Ca is yielded between 2.8×10−3 and 1.63. If ucw is defined as the mean velocity of uw from discharge jw, Ca is yielded between 7.8710 and 0.03. Although there are discussions on the definition validity of the pore-scale capillary number [42, 43], a transition of a capillary-dominated regime to a viscous-dominated one is indeed observed in our simulations.


### 3.3. Droplet Shape and Size Distribution

It is expected that the volume and shape of droplets reveal the effect of the competition between the viscous and capillary stresses on droplet mobilities. Such a competition is qualitatively depicted in Figure3(a), where droplet distributions of three scenarios at Snw=0.05 at a later stage are reported with respect to σlb=0.0005, 0.005, and 0.05. When the capillary stresses dominate over the viscous stresses (at σlb=0.05, rendered in red in Figure 3(a)), droplets are mainly rounded. Under such circumstances, droplets with radii similar to the ones of pores in which they reside are likely to be capillary-trapped. At lower surface tensions (rendered in blue and green in Figure 2(a)), the viscous stresses tend to overwhelm the capillary ones. The dominant of viscous stress facilitates the elongation of medium-size droplets along the main flow direction, while large droplets can be heavily squeezed and stretched during their creeping movements from one pore to the other. As a result, the snap-off of droplets often occurs, as indicated by Rossen [44]. In summary, the competition between viscous and capillary stresses governs the droplet mobilities, which influence the overall mass discharge, as reported in Figure 1(a). We can conclude that the smaller the surface tension, the lower the energy dissipation and therefore the higher the total discharge. Moreover, as surface tension decreases, snap-off and creeping movement of droplets ease the sawlike disturbance in fluid discharge due to the capture and release of droplets.Figure 3
(a) Snapshots of non-wetting phase (droplets) at steady state for different surface tensions (σlb=5×10−2, 5×10−3, and 5×10−4) and Snw=0.05. (b) Characterization of droplet shape at different surface tensions using volume-to-surface-area, Vnw/Sfnw, versus “spherical radius” of individual droplets, 3Vnw/4π1/3, where Vnw and Sfnw are the volume and surface area of individual droplets. For a spherical droplet, Vnw/Sfnw=3Vnw/4π1/3, indicated by the black solid line. Results are plotted for all simulations with saturations of 0.8≤Sw≤0.95. Inset: individual droplet volumes normalized by the averaged pore size of the medium are plotted against the apparent Capillary number, Ca∗=SwFlbpSfnw/σlb, for these three surface tensions. The size of the symbols is proportional to the reciprocal of droplet saturation; i.e., the higher the droplet saturation, the smaller the symbol.
(a)
(b)By further varyingSnw, two distinguishable dynamic behaviors of the non-wetting phase can be identified: droplet percolation and ganglion snap-off. Droplet percolation occurs at low non-wetting phase saturations (Snw<0.1−0.15, depending on σlb). In this regime, small droplets are passively transported by the wetting fluid and they moderately interact with the solid geometry. In contrast, large droplets strongly interact with the solid geometry, move slowly through the pores, or can be even capillary-trapped. Overall, large droplets act as collectors of smaller and faster ones. As Snw increases, large ganglia (i.e., collections of connected droplets) can be formed, some of which will eventually percolate the entire system. Such ganglia are likely to snap off and reconnect again over time. Our simulations suggest that for Snw>0.2 (depending on σlb), ganglion snap-off and redevelopment represent the most important mechanism of droplet mobilities. Moreover, snap-off is more frequent at lower surface tension (low Ca, see movies in the supplementary material). Our observation is in line with Deng et al. [45], stating that high capillary stresses inhibit snap-off.Figure3(b) reports the distributions of volume Vnw to surface Sfnw ratio of droplets with respect to the hypothetical volume-to-surface ratio for a perfect sphere (black line with a slope of 1). Such analyses are performed for σlb=0.0005, 0.005, and 0.05 and for Snw up to 0.2. This is because at higher Snw, droplets mostly develop into large ganglia (i.e., collections of coalesced droplets) that strongly skew these distributions. Droplet statistics are evaluated when the system reaches a time-averaged steady state (either at 3×105 or 6×105 iterations, depending on σlb). As expected, the lower the surface tension, the smaller the droplet volume-to-surface ratio (i.e., the higher the droplet deformation). In contrast, high surface tension results in more spherical droplets which follow the black line with a slope of 1. The deviation away from the black line is induced by the size of the droplets and the saturation of the non-wetting phase (up to Snw≥0.15), where droplets are sensibly larger than the average pore size and therefore more deformed. This deviation enhances the coalescence of droplets [46] and forces the volume-to-surface ratio to deviate from the black line in Figure 3(b).The effect of pore topology on droplet volume distribution is documented in the inset of Figure3(b), where the normalized droplet volume (to the averaged pore volume) is plotted versus an apparent Capillary number Ca∗ which is formulated as Ca∗=SwFlbpSfnw/σlb. Here, Ca∗ is an analog of the dimensionless Bond number, Bo=Δρgr2/σ (Δρ is the density difference between fluid phases, g is the gravitational acceleration, r is the droplet radius, and σ is the surface tension), where r2 is approximated by the droplet surface area Sfnw and Δρg is replaced by SwFlbp which represents the effective hydraulic gradient applied to the system and exerted as a body force on the wetting phase only. A regression analysis suggests that the normalized droplet volume is proportional to Ca∗, where γ=1.259, 1.294, and 1.375, with respect to σlb=5×10−4, 5×10−3, and 5×10−2. Note that the exponent γ is bounded by γ=1.50 which is obtained through a theoretical calculation on a spherical droplet. The decrease in exponent γ due to the decrease in surface tension further confirms the effect of lower surface tensions on droplet stretching and breaking (into smaller droplets, i.e., snap-off). However, the change in γ is rather subtle such that a much higher weight can be likely attributed to the pore topology effect on γ. Our interpretation is in agreement with Shimizu and Tanaka [25], discussing the dominant effect of the solid-matrix topology on phase separation in porous media.


### 3.4. Implications to Geothermal, Hydrothermal, and Magmatic Systems

Up till now, we have demonstrated the following observations that can be concluded from our simulations on phase separation driven by a background hydraulic gradient in a porous medium:(a)
The exsolved phase (hereafters droplet) prefers its settlements in regions of low advection capacity, suggesting a transport-enhanced mechanism for droplet growth. Because of the trapping of droplets, the overall resistance to flow increases
(b)
System clogging varies with capillary-viscous stress competition, but overall clogging prefers at high non-wetting phase saturation. In other words, the higher the phase proportion of droplets, the higher the potential of system clogging
(c)
The fluid-matrix interaction strongly influences the normalized droplet volume distribution that follows a power law. The geometry of the pore space plays an important role in the geometry of the non-wetting fluid and hence in the cloggingIt is well-known that phase separation has been identified as a critical process that strongly influences the efficiency of heat and mass transfer in geothermal (particularly high-enthalpy) [47–49] and hydrothermal systems [12, 50–53]. In geothermal and hydrothermal systems, phase separation may take place under a variety of circumstances, typically including subcritical (i.e., boiling) and supercritical (condensation) phase separations. The former one occurs below the critical point where a vapor will be separated from the brine, whereas the latter one occurs above the critical point where a dense, highly saline brine will be separated [51]. Despite the evidence and the importance of clogging due to phase separation, the impact of clogging on geofluid circulation during phase separation has not been considered. One of the key impacts due to clogging is the reduction in permeability which has been identified as a critical parameter for geofluid circulation. Our results suggest that the reduction in permeability depends on the velocity distribution (due to the transport-enhanced mechanism), fluid composition (phase saturation), and pore space geometry.In magmatic systems, a so-called gas-driven filter-pressing has been proposed to elucidate the influence of bubble exsolution and growth on the melt extraction from crystal-rich zones [54–57]. As the exsolved bubbles nucleate and grow in the pore space in magma reservoirs (largely made of crystals and silicate melt in the interstices, called “crystal mushes”), the volume taken by the exsolved bubbles is expected to push out the melt and form the so-called crystal-poor rhyolite/obsidian flow. However, the clogging effect of bubbles is certainly not taken into consideration in these gas-driven filter-pressing scenarios [54–57]. Our results show that clogging by bubbles is important and must be considered in estimating the efficiency of gas-driven filter-pressing (particularly for low Bond number situations). In particular, the more the bubbles form, the slower the phase separation will be, ultimately limiting the efficiency of phase separation.


## 3.1. Transport Stage and Mechanism

A typical three-stage evolution of the phase separation is clearly demonstrated in our simulations (Figure1): transport of mixture in stage I, transport of growing droplets in stage II, and steady-state transport of two phases in stage III. During stage I, the two-fluid mixture quickly establishes an early/initial steady-state discharge, jztot, which depends on Flbp and the wetting fluid volume fraction Sw (Figure 1(a)). At this stage, the two phases are still completely mixed (i.e., droplets are not yet formed). This complete mixing is indicated by the zero discharge, jznw, of the non-wetting phase. At the beginning of stage II, phase separation is triggered and droplets start to develop via the spinodal decomposition, dominated by a diffusion growth (Figure 1(a)). A few thousand iterations later, droplets hydrodynamically interact with each other, and they grow via coalescence. As soon as the droplets reach sizes similar to the ones of pores and throats, capillary stresses start affecting the fluid flow. These observations are supported by the confirmation of the increase in jznw (due to the growth and transport of droplets) and the decrease in total discharge, jZtot, induced by the capacity reduction of flow pathways which results in higher energy dissipation. At late times of stage II, droplets can reach sizes similar to or even bigger than the ones of pores and throats. A typical example of droplet evolution is reported by taking snapshots at early and late times during stage II in one of our simulations (Figure 1(c)). A more comprehensive understanding of the phase separation process at early times can be gained by watching movies included in our supplementary materials (available here). In stage III, both the total and non-wetting discharges eventually reach their steady-state values. During this stage, transport of large droplets is strongly controlled by the competition between capillary and shear stresses. At high surface tension (indicated by the red line in Figure 1(a)), large droplets are mostly trapped and they act as immobile obstacles to the flow of the wetting phase. At small surface tensions (indicated by the blue and green curves in Figure 1(a)), droplets are likely to be deformed and they can migrate under the driving of the wetting phase. The latter case results in much higher total and non-wetting discharges. Nonetheless, both the total and non-wetting discharges fluctuate around their steady-state values due to the capture and release of droplets. Apparently, this fluctuation is lessened due to snap-off and reconnection of droplets as surface tension decreases, which will be further explained in the subsequent subsections.Since phase separation is occurring while the fluid mixture is being transmitted, droplet growth is expected to be strongly influenced by the connectivity of the porous medium and the advection of fluids. This observation is confirmed in our calculations that growth of droplets is governed by a transport-enhanced mechanism, as described by Zuo et al. [10]. Such behavior is clearly visible (Figure 1(b)) when comparing the evolution of saturation distribution Snw of the non-wetting phase between two particular scenarios, one with a null hydraulic gradient (i.e., droplet growth is diffusion-dominated) and one with a non-null hydraulic gradient (i.e., droplet growth is influenced by fluid advection). To interpret such behavior, for both scenarios, the saturation of the non-wetting phase, Snw, is assembled as a function of the normalized z-component of fluid velocity. The fluid velocity is obtained in a steady-state single-phase fluid flow calculation with the same geometry and hydraulic gradient. It is rational to approximate that the strength/capacity of advection is proportional to velocity magnitude, i.e., the higher the z-velocity, the higher the advection capacity. For the scenario of null hydraulic gradient (i.e., diffusion-controlled and affected at a lesser extent by geometrical properties), Snw quickly becomes a roughly uniform distribution, where non-wetting droplets mainly reside in regions of median advection capacity (Figure 1(b)). For the scenario of non-null hydraulic gradient, a much more rich temporal evolution of Snw is registered, especially at the early times when small droplets are transported in the highly advective regions (i.e., the most well-connected pores and throats) but eventually coalesce to form large droplets which reside in regions of low (but non-null) advection capacity and likely to be capillary-trapped.


## 3.2. Connectivity of Wetting Fluid

As stated in previous sections, a steady-state velocity field is reached during stage III. This steady-state velocityuw of the wetting phase is utilized to compute the continuous streamlines of the wetting phase only. Here, the continuous streamline is defined as a streamline continuously spanning from the inlet to outlet boundaries along the main flow direction. By definition, in the lattice-Boltzmann method, the velocity uwx,y,z is calculated as the ratio of momentum over fluid density, i.e., uw=jw/ρw. Then, a seepage velocity of the wetting phase in the main flow direction (z-direction) U¯zw is calculated along the streamlines. Note that streamlines are generally tortuously passing through the porous medium and thus have a length longer than the 22domain size in the z-direction [41]. Mathematically, the z-direction seepage velocity of the wetting phase is formulated as
(1)U¯zw=1N∑s=1N∑i=1Lsuzwxi,yi,zi,where Ls is the length of the s-th continuous streamline and N is the total number of the continuous streamlines that can be identified and pass through the porous medium from the initial cross section z=0 in an individual simulation. A total velocity, U¯w, along all the continuous streamlines is calculated accordingly:
(2)U¯w=1N∑s=1N∑i=1Lsuxwxi,yi,zi2+uywxi,yi,zi2+uzwxi,yi,zi2.Starting from the initial cross sectionz=0, velocities along the streamlines are sampled at the equidistant points, using an integration step equals to 0.1 of a lattice to ensure a good estimation of the continuous streamlines passing through the pore space. The presence of a notable amount of continuous streamlines indicates that (1) continuous pathways still exist in the wetting phase (i.e., it is physically meaningful to use body force to formulate a background hydraulic gradient) and (2) the porous medium is not yet completely clogged by the non-wetting phase. A completely clogged stage is reached when the wetting phase retains very limited number of continuous streamlines (up to 100). Note that there are more than 8000 continuous wetting-fluid streamlines for a single-phase flow in the current porous medium. At the completely clogged stage, the trapped or slowly moving non-wetting fluid breaks up the wetting phase into several disconnected parts inside the porous medium and thus impedes the transmit of the wetting phase. Figure 2(b) shows the completely clogged stages for five different surface tensions σlb, indicated by arrows with different colors. Our results clearly indicate that the completely clogged stage occurs at higher wetting-fluid saturation Sw for higher σlb, due to higher resistance introduced by capillary stresses.Figure 2
(a) Magnitude of the norm of total velocity throughout the porous medium (solid in back and gray) as well as non-wetting droplets (in red) atSnw=0.05 and σlb=0.05 at a given time snapshot (3×105 iterations). (b) Mean seepage velocity, U¯zw, averaged along the continuous streamlines of the wetting phase. (c) and (d) Magnitude of the norm of seepage velocity (solids are rendered transparently in gray) and non-wetting droplets (in red) for Snw=0.05 (left) and Snw=0.20 (right) with σlb=0.05 at 3×105 iterations. For better visualization, not all streamlines are reported.
(a)
(b)
(c)
(d)Since the viscosity of both phases is fixed, the competition between capillary and viscous stresses is therefore explored by uniquely varyingσlb. This competition can be cast into a dimensionless capillary number Ca=ucwvlb/σlb, where ucw is the characteristic velocity of the wetting phase. If ucw is defined as the mean seepage velocity, U¯zw, of the wetting phase, i.e., before the system is completely clogged, Ca is yielded between 2.8×10−3 and 1.63. If ucw is defined as the mean velocity of uw from discharge jw, Ca is yielded between 7.8710 and 0.03. Although there are discussions on the definition validity of the pore-scale capillary number [42, 43], a transition of a capillary-dominated regime to a viscous-dominated one is indeed observed in our simulations.


## 3.3. Droplet Shape and Size Distribution

It is expected that the volume and shape of droplets reveal the effect of the competition between the viscous and capillary stresses on droplet mobilities. Such a competition is qualitatively depicted in Figure3(a), where droplet distributions of three scenarios at Snw=0.05 at a later stage are reported with respect to σlb=0.0005, 0.005, and 0.05. When the capillary stresses dominate over the viscous stresses (at σlb=0.05, rendered in red in Figure 3(a)), droplets are mainly rounded. Under such circumstances, droplets with radii similar to the ones of pores in which they reside are likely to be capillary-trapped. At lower surface tensions (rendered in blue and green in Figure 2(a)), the viscous stresses tend to overwhelm the capillary ones. The dominant of viscous stress facilitates the elongation of medium-size droplets along the main flow direction, while large droplets can be heavily squeezed and stretched during their creeping movements from one pore to the other. As a result, the snap-off of droplets often occurs, as indicated by Rossen [44]. In summary, the competition between viscous and capillary stresses governs the droplet mobilities, which influence the overall mass discharge, as reported in Figure 1(a). We can conclude that the smaller the surface tension, the lower the energy dissipation and therefore the higher the total discharge. Moreover, as surface tension decreases, snap-off and creeping movement of droplets ease the sawlike disturbance in fluid discharge due to the capture and release of droplets.Figure 3
(a) Snapshots of non-wetting phase (droplets) at steady state for different surface tensions (σlb=5×10−2, 5×10−3, and 5×10−4) and Snw=0.05. (b) Characterization of droplet shape at different surface tensions using volume-to-surface-area, Vnw/Sfnw, versus “spherical radius” of individual droplets, 3Vnw/4π1/3, where Vnw and Sfnw are the volume and surface area of individual droplets. For a spherical droplet, Vnw/Sfnw=3Vnw/4π1/3, indicated by the black solid line. Results are plotted for all simulations with saturations of 0.8≤Sw≤0.95. Inset: individual droplet volumes normalized by the averaged pore size of the medium are plotted against the apparent Capillary number, Ca∗=SwFlbpSfnw/σlb, for these three surface tensions. The size of the symbols is proportional to the reciprocal of droplet saturation; i.e., the higher the droplet saturation, the smaller the symbol.
(a)
(b)By further varyingSnw, two distinguishable dynamic behaviors of the non-wetting phase can be identified: droplet percolation and ganglion snap-off. Droplet percolation occurs at low non-wetting phase saturations (Snw<0.1−0.15, depending on σlb). In this regime, small droplets are passively transported by the wetting fluid and they moderately interact with the solid geometry. In contrast, large droplets strongly interact with the solid geometry, move slowly through the pores, or can be even capillary-trapped. Overall, large droplets act as collectors of smaller and faster ones. As Snw increases, large ganglia (i.e., collections of connected droplets) can be formed, some of which will eventually percolate the entire system. Such ganglia are likely to snap off and reconnect again over time. Our simulations suggest that for Snw>0.2 (depending on σlb), ganglion snap-off and redevelopment represent the most important mechanism of droplet mobilities. Moreover, snap-off is more frequent at lower surface tension (low Ca, see movies in the supplementary material). Our observation is in line with Deng et al. [45], stating that high capillary stresses inhibit snap-off.Figure3(b) reports the distributions of volume Vnw to surface Sfnw ratio of droplets with respect to the hypothetical volume-to-surface ratio for a perfect sphere (black line with a slope of 1). Such analyses are performed for σlb=0.0005, 0.005, and 0.05 and for Snw up to 0.2. This is because at higher Snw, droplets mostly develop into large ganglia (i.e., collections of coalesced droplets) that strongly skew these distributions. Droplet statistics are evaluated when the system reaches a time-averaged steady state (either at 3×105 or 6×105 iterations, depending on σlb). As expected, the lower the surface tension, the smaller the droplet volume-to-surface ratio (i.e., the higher the droplet deformation). In contrast, high surface tension results in more spherical droplets which follow the black line with a slope of 1. The deviation away from the black line is induced by the size of the droplets and the saturation of the non-wetting phase (up to Snw≥0.15), where droplets are sensibly larger than the average pore size and therefore more deformed. This deviation enhances the coalescence of droplets [46] and forces the volume-to-surface ratio to deviate from the black line in Figure 3(b).The effect of pore topology on droplet volume distribution is documented in the inset of Figure3(b), where the normalized droplet volume (to the averaged pore volume) is plotted versus an apparent Capillary number Ca∗ which is formulated as Ca∗=SwFlbpSfnw/σlb. Here, Ca∗ is an analog of the dimensionless Bond number, Bo=Δρgr2/σ (Δρ is the density difference between fluid phases, g is the gravitational acceleration, r is the droplet radius, and σ is the surface tension), where r2 is approximated by the droplet surface area Sfnw and Δρg is replaced by SwFlbp which represents the effective hydraulic gradient applied to the system and exerted as a body force on the wetting phase only. A regression analysis suggests that the normalized droplet volume is proportional to Ca∗, where γ=1.259, 1.294, and 1.375, with respect to σlb=5×10−4, 5×10−3, and 5×10−2. Note that the exponent γ is bounded by γ=1.50 which is obtained through a theoretical calculation on a spherical droplet. The decrease in exponent γ due to the decrease in surface tension further confirms the effect of lower surface tensions on droplet stretching and breaking (into smaller droplets, i.e., snap-off). However, the change in γ is rather subtle such that a much higher weight can be likely attributed to the pore topology effect on γ. Our interpretation is in agreement with Shimizu and Tanaka [25], discussing the dominant effect of the solid-matrix topology on phase separation in porous media.


## 3.4. Implications to Geothermal, Hydrothermal, and Magmatic Systems

Up till now, we have demonstrated the following observations that can be concluded from our simulations on phase separation driven by a background hydraulic gradient in a porous medium:(a)
The exsolved phase (hereafters droplet) prefers its settlements in regions of low advection capacity, suggesting a transport-enhanced mechanism for droplet growth. Because of the trapping of droplets, the overall resistance to flow increases
(b)
System clogging varies with capillary-viscous stress competition, but overall clogging prefers at high non-wetting phase saturation. In other words, the higher the phase proportion of droplets, the higher the potential of system clogging
(c)
The fluid-matrix interaction strongly influences the normalized droplet volume distribution that follows a power law. The geometry of the pore space plays an important role in the geometry of the non-wetting fluid and hence in the cloggingIt is well-known that phase separation has been identified as a critical process that strongly influences the efficiency of heat and mass transfer in geothermal (particularly high-enthalpy) [47–49] and hydrothermal systems [12, 50–53]. In geothermal and hydrothermal systems, phase separation may take place under a variety of circumstances, typically including subcritical (i.e., boiling) and supercritical (condensation) phase separations. The former one occurs below the critical point where a vapor will be separated from the brine, whereas the latter one occurs above the critical point where a dense, highly saline brine will be separated [51]. Despite the evidence and the importance of clogging due to phase separation, the impact of clogging on geofluid circulation during phase separation has not been considered. One of the key impacts due to clogging is the reduction in permeability which has been identified as a critical parameter for geofluid circulation. Our results suggest that the reduction in permeability depends on the velocity distribution (due to the transport-enhanced mechanism), fluid composition (phase saturation), and pore space geometry.In magmatic systems, a so-called gas-driven filter-pressing has been proposed to elucidate the influence of bubble exsolution and growth on the melt extraction from crystal-rich zones [54–57]. As the exsolved bubbles nucleate and grow in the pore space in magma reservoirs (largely made of crystals and silicate melt in the interstices, called “crystal mushes”), the volume taken by the exsolved bubbles is expected to push out the melt and form the so-called crystal-poor rhyolite/obsidian flow. However, the clogging effect of bubbles is certainly not taken into consideration in these gas-driven filter-pressing scenarios [54–57]. Our results show that clogging by bubbles is important and must be considered in estimating the efficiency of gas-driven filter-pressing (particularly for low Bond number situations). In particular, the more the bubbles form, the slower the phase separation will be, ultimately limiting the efficiency of phase separation.


## 4. Conclusions

In summary, a series of lattice-Boltzmann simulations of multiphase flow has been performed in a porous medium that is generated via a crystal nucleation and growth algorithm. These simulations address the hydrodynamic response when phase exsolution/separation occurs in a flow through this medium driven by a hydraulic gradient. The focus of this study is limited on the effect of variations in surface tension and fluid volume fractions. For convenience, the phase separation process is simplified by taking the spinodal decomposition to describe a rapid demixing of a mixture of phases to form immiscible droplets out of a solution. Three stages have been identified during the evolution of phase separation, namely, transport of mixture in stage I, transport of growing droplets in stage II, and steady-state transport of two phases in stage III. Moreover, large droplets or clusters of droplets tend to settle in regions of low advection capacity, suggesting a transport-enhanced mechanism for droplet growth. Two distinguishable dynamic behaviors of droplet transport can be recognized, namely, droplet percolation at low non-wetting phase saturation and ganglion snap-off at high non-wetting phase saturation. The snap-off and formation of ganglion represent the most important features of droplet mobility especially at high non-wetting phase content. These large droplets or ganglia move slowly through or even are trapped in the pores due to capillary resistance. As surface tension increases, the breakup chance of the wetting phase increases and ultimately a completely clogged stage arises in the system. This completely clogged stage can be quantitatively defined using an averaged flux of the wetting phase calculated along continuous streamlines. The completely clogged stage occurs at higher wetting-fluid saturation for higher surface tension, due to higher resistance introduced by capillary stresses. Both the decrease in surface tension and the saturation increase in the non-wetting phase enhance the coalescence of droplets and thus the deviation of the volume-to-surface ratio to the one of spherical droplets. A regression analysis suggests that the normalized droplet volume is proportional to the power of an apparent Capillary Number. However, we tend to attribute a much higher weight of pore topology, rather than surface tension, on the exponent change of this power. Our results provide insights to the understanding of phase separation in natural systems, such as gassy soils [20–22], CO2 storage reservoirs [7–10], hydrothermal systems [11, 19], and water-rich magma reservoirs [16]. In all of these settings, the emergence, growth, and transport of the secondary phase introduce perturbations to the fluid flow driven by background pressure gradients, which are likely critical to the dynamic evolution of these settings.



---

*Source: 5176410-2019-05-14.xml*