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

Parmigiani, Andrea; Di Palma, Paolo Roberto; Leclaire, Sébastien; Habib, Faraz; Kong, Xiang-Zhao

doi:https://doi.org/10.1155/2019/5176410

Geofluids

On this page

Abstract Introduction Results Conclusions Appendix Data Availability Conflicts of Interest Acknowledgments Supplementary Materials References Copyright Related Articles

Special Issue

Geomechanics and Fluid Flow in Geothermal Systems

View this Special Issue

Research Article | Open Access

Volume 2019 | Article ID 5176410 | https://doi.org/10.1155/2019/5176410

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

Andrea Parmigiani,²Paolo Roberto Di Palma,³Sébastien Leclaire,⁴Faraz Habib,¹and Xiang-Zhao Kong¹

Guest Editor: Francesco Parisio

Received29 Nov 2018

Accepted06 Mar 2019

Published14 May 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.

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) CO₂ sequestration [4–10] where the dissolved CO₂ 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 populations

Compared 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 , 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, , is defined using the density fields of the wetting fluid, , and the non-wetting fluid, . Regions/lattice nodes occupied by the wetting and the non-wetting fluids are indicated by and , respectively (e.g., droplets in Figure 1(c) are the volume rendering of ). The flag is also used to retrieve the total momentum field, , and the individual momentum fields, and , of the wetting and non-wetting fluids, respectively. Integration of the vertical components of the momentum fields, namely, , , and , over the whole simulation domain yields the total discharge, , and its individual contributions, and , of the two immiscible fluids, respectively (Figure 1(a)).

(a)

(b)

(c)

Figure 1

(a) Three-stage (i.e., stages I, II, and III) temporal evolution of total discharge, , and non-wetting fluid discharge, , for three values and . (b) Temporal evolution of distribution (simulation time is normalized to the maximum lattice time of ) as a function of porous medium conductivity for null (middle panel) and non-null (bottom panel) hydraulic gradients with and . 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 , , and .

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 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, .

A set of simulations with a non-wetting phase saturation of and a surface tension of (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 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 (in lattice units), which corresponds to a relaxation time of 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 on the wetting phase only [40]. Therefore, the total driving force on the system is proportional to the wetting phase volume, . 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 (Figure 1): 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, , which depends on and the wetting fluid volume fraction (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, , 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 (due to the growth and transport of droplets) and the decrease in total discharge, , 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 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, , is assembled as a function of the normalized -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 -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), 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 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 velocity 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 is calculated as the ratio of momentum over fluid density, i.e., . Then, a seepage velocity of the wetting phase in the main flow direction (-direction) 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 -direction [41]. Mathematically, the -direction seepage velocity of the wetting phase is formulated as where is the length of the -th continuous streamline and is the total number of the continuous streamlines that can be identified and pass through the porous medium from the initial cross section in an individual simulation. A total velocity, , along all the continuous streamlines is calculated accordingly:

Starting from the initial cross section , 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 , indicated by arrows with different colors. Our results clearly indicate that the completely clogged stage occurs at higher wetting-fluid saturation for higher , due to higher resistance introduced by capillary stresses.

(a)

(b)

(c)

(d)

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) at and at a given time snapshot ( iterations). (b) Mean seepage velocity, , 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 (left) and (right) with at iterations. For better visualization, not all streamlines are reported.

Since the viscosity of both phases is fixed, the competition between capillary and viscous stresses is therefore explored by uniquely varying . This competition can be cast into a dimensionless capillary number , where is the characteristic velocity of the wetting phase. If is defined as the mean seepage velocity, , of the wetting phase, i.e., before the system is completely clogged, is yielded between and . If is defined as the mean velocity of from discharge , 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 Figure 3(a), where droplet distributions of three scenarios at at a later stage are reported with respect to , 0.005, and 0.05. When the capillary stresses dominate over the viscous stresses (at , 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.

(a)

(b)

Figure 3

(a) Snapshots of non-wetting phase (droplets) at steady state for different surface tensions (, , and ) and . (b) Characterization of droplet shape at different surface tensions using volume-to-surface-area, , versus “spherical radius” of individual droplets, , where and are the volume and surface area of individual droplets. For a spherical droplet, , indicated by the black solid line. Results are plotted for all simulations with saturations of . Inset: individual droplet volumes normalized by the averaged pore size of the medium are plotted against the apparent Capillary number, , 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.

By further varying , 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 (, depending on ). 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 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 (depending on ), ganglion snap-off and redevelopment represent the most important mechanism of droplet mobilities. Moreover, snap-off is more frequent at lower surface tension (low , see movies in the supplementary material). Our observation is in line with Deng et al. [45], stating that high capillary stresses inhibit snap-off.

Figure 3(b) reports the distributions of volume to surface 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 , 0.005, and 0.05 and for up to 0.2. This is because at higher , 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 or iterations, depending on ). 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 ), 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 Figure 3(b), where the normalized droplet volume (to the averaged pore volume) is plotted versus an apparent Capillary number which is formulated as . Here, is an analog of the dimensionless Bond number, ( is the density difference between fluid phases, is the gravitational acceleration, is the droplet radius, and is the surface tension), where is approximated by the droplet surface area and is replaced by 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 , where , 1.294, and 1.375, with respect to , , and . Note that the exponent is bounded by 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 clogging

It 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], CO₂ 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.

Appendix

Color Gradient Multiphase (CGM) Lattice Boltzmann Model

The core algorithm of our CGM approach follows the D3Q15 two-phase model of [26]. The two-phase fluid flows are resolved using two sets of distribution functions, one for each fluid, moving on the D3Q15 lattice with the velocity vectors as given in [26]. Note that a rescale between the lattice and the physical units is always needed, even though the lattice system is defined as , where is the lattice spacing and is the lattice time step.

The distribution functions for a fluid of color (with ) are denoted as , while is the color-blind distribution function. Excluding the boundary condition, the evolution algorithm of the fluid bulk follows the operators below: (1)Collision operator:(2)Streaming operator:where the symbol denotes the bra Dirac notation with an expansion with respect to the velocity space. The collision operator results in three main operations: the single-phase collision operator , the multiphase perturbation operator , and the multiphase recoloring operator

The first operator is based on the standard single-relaxation-time operator of the single-phase LB model

The density of the fluid is given by the first moment of the distribution functions where the superscript denotes the equilibrium. The total fluid density is given by , while the total momentum is defined as the second moment of the distribution functions: where is the velocity of the color-blind distribution functions. The general form of the D3Q15 equilibrium functions are defined in [26].

The term in Eq. (A.3) is a density distribution function modification designed to add external forces. In this study, this term is solely used to change momentum of the dispersed phase: where is an all-zero moment vector except for the indexes corresponding to the momentum such that

Here, , so that the constant forcing is applied on a single-component (here the wetting phase only) [35] pore and where the momentum indexes , , and are lattice-dependent and are given in [26]. Note that all these modifications are made in the moment space before being translated into the distribution space with the matrix multiplication . The matrix is defined in [26] for the D3Q15 lattice.

The effective relaxation parameter is defined so that the fluid viscosity is consistent with the macroscopic equations for a single-phase flow in the single-phase regions. When the viscosities of the fluids are different, an interpolation is applied to define the parameter at the interface. If is the kinematic viscosity of the fluid , we define the viscosity at the interface between the fluids

Then, the effective relaxation parameter is

The surface tension in the CGM is modeled by means of the perturbation operator [58, 59]. It takes the following form: where the color gradient approximates the normal to the interface and are lattice-dependent weights given in [26]. Reis and Phillips [59] and Liu et al. [60] have shown that this operator complies within the macroscopic limit, with the capillary stress tensor present in the macroscopic equations for two-phase flows if the weights are well chosen. The parameter is space- and time-dependent and is chosen to fit the surface tension value at the fluid interface: with being the interfacial surface tension. Although this operator generates the surface tension, it does not guarantee the fluid’s immiscibility. To minimize mixing and segregate the fluids, the recoloring operator needs to be properly selected.

This operator is used to maximize the amount of fluid at the interface that is sent to the fluid region, while remaining consistent with the laws of conservation of mass and total momentum. The recoloring operator presented here is based on Refs. [61, 62] and is as follows: where is a parameter controlling the thickness of the interface. The variable corresponds to the angle between the color gradient and the lattice velocity . A special relation for exists for which it is possible to easily control the spurious currents and the physical width of the interface with lattice refinement [63].

For more information on the CGM model, the reader may refer to [26].

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

X.-Z. Kong acknowledges support from ETH Research Grants ETH-12 15-2 and ETH-02 16-2 and Swiss National Science Foundation Grant 172760. A. Parmigiani was supported by a Swiss National Science Foundation Ambizione grant. We also thank Prof. Olivier Bachmann for his helpful suggestions on the magmatic systems. For a distribution of model, please contact Dr. A. Parmigiani (Email: [email protected]) and Dr. X.-Z. Kong (Email: [email protected]).

Supplementary Materials

There are three movies which describe the phase separation process in three different scenarios: (i) BF_2e-4_Sigma_5e-2_Snw_5e-2.gif shows the separation process of the non-wetting phase rendered in red with a hydraulic gradient of , a surface tension , and a non-wetting phase saturation . (ii) BF_2e-4_Sigma_5e-4_Snw_5e-2.gif shows the separation process of the non-wetting phase rendered in green with a hydraulic gradient of , a surface tension , and a non-wetting phase saturation . (iii) nullBF_Sigma_5e-2_Snw_5e-2.gif shows the separation process of the non-wetting phase rendered in blue with a hydraulic gradient of , a surface tension , and a non-wetting phase saturation . (Supplementary Materials)

References

A. Liebscher, “Aqueous fluids at elevated pressure and temperature,” Geofluids, vol. 10, no. 1-2, 19 pages, 2010.
View at: Publisher Site | Google Scholar
I. Leifer and R. K. Patro, “The bubble mechanism for methane transport from the shallow sea bed to the surface: a review and sensitivity study,” Continental Shelf Research, vol. 22, no. 16, pp. 2409–2428, 2002.
View at: Publisher Site | Google Scholar
N. Mahabadi, X. Zheng, T. S. Yun, L. van Paassen, and J. Jang, “Gas bubble migration and trapping in porous media: pore-scale simulation,” Journal of Geophysical Research: Solid Earth, vol. 123, no. 2, pp. 1060–1071, 2018.
View at: Publisher Site | Google Scholar
J. M. Hunt, “Generation and migration of petroleum from abnormally pressured fluid compartments (1),” AAPG Bulletin, vol. 74, no. 1, pp. 1–12, 1990.
View at: Publisher Site | Google Scholar
L. Zuo, S. Krevor, R. W. Falta, and S. M. Benson, “An experimental study of CO₂ exsolution and relative permeability measurements during CO₂ saturated water depressurization,” Transport in Porous Media, vol. 91, no. 2, pp. 459–478, 2012.
View at: Publisher Site | Google Scholar
M. Plampin, T. Illangasekare, T. Sakaki, and R. Pawar, “Experimental study of gas evolution in heterogeneous shallow subsurface formations during leakage of stored CO₂,” International Journal of Greenhouse Gas Control, vol. 22, pp. 47–62, 2014.
View at: Publisher Site | Google Scholar
L. Zuo and S. M. Benson, “CO₂ exsolution: challenges and opportunities in subsurface flow management,” Energy Procedia, vol. 63, pp. 5664–5670, 2014.
View at: Publisher Site | Google Scholar
B. M. Tutolo, A. J. Luhmann, X. Z. Kong, M. O. Saar, and W. E. Seyfried Jr., “Experimental observation of permeability changes in dolomite at CO₂ sequestration conditions,” Environmental Science & Technology, vol. 48, no. 4, pp. 2445–2452, 2014.
View at: Publisher Site | Google Scholar
B. M. Tutolo, X. Z. Kong, W. E. Seyfried Jr., and M. O. Saar, “High performance reactive transport simulations examining the effects of thermal, hydraulic, and chemical (THC) gradients on fluid injectivity at carbonate CCUS reservoir scales,” International Journal of Greenhouse Gas Control, vol. 39, pp. 285–301, 2015.
View at: Publisher Site | Google Scholar
L. Zuo, J. B. Ajo-Franklin, M. Voltolini, J. T. Geller, and S. M. Benson, “Pore-scale multiphase flow modeling and imaging of CO₂ exsolution in sandstone,” Journal of Petroleum Science and Engineering, vol. 155, pp. 63–77, 2017.
View at: Publisher Site | Google Scholar
S. Emmanuel and B. Berkowitz, “Phase separation and convection in heterogeneous porous media: implications for seafloor hydrothermal systems,” Journal of Geophysical Research, vol. 112, no. B5, 2007.
View at: Publisher Site | Google Scholar
D. Coumou, T. Driesner, P. Weis, and C. A. Heinrich, “Phase separation, brine formation, and salinity variation at black smoker hydrothermal systems,” Journal of Geophysical Research, vol. 114, no. B3, 2009.
View at: Publisher Site | Google Scholar
J. Choi and R. P. Lowell, “The response of two-phase hydrothermal systems to changing magmatic heat input at mid-ocean ridges,” Deep Sea Research Part II: Topical Studies in Oceanography, vol. 121, pp. 17–30, 2015.
View at: Publisher Site | Google Scholar
H. E. Huppert and A. W. Woods, “The role of volatiles in magma chamber dynamics,” Nature, vol. 420, no. 6915, pp. 493–495, 2002.
View at: Publisher Site | Google Scholar
W. Degruyter and C. Huber, “A model for eruption frequency of upper crustal silicic magma chambers,” Earth and Planetary Science Letters, vol. 403, pp. 117–130, 2014.
View at: Publisher Site | Google Scholar
A. Parmigiani, W. Degruyter, S. Leclaire, C. Huber, and O. Bachmann, “The mechanics of shallow magma reservoir outgassing,” Geochemistry, Geophysics, Geosystems, vol. 18, no. 8, pp. 2887–2905, 2017.
View at: Publisher Site | Google Scholar
J. Jiménez-Martínez, T. Le Borgne, H. Tabuteau, and Y. Méheust, “Impact of saturation on dispersion and mixing in porous media: photobleaching pulse injection experiments and shear-enhanced mixing model,” Water Resources Research, vol. 53, no. 2, pp. 1457–1472, 2017.
View at: Publisher Site | Google Scholar
T. Xia, E. Dontsov, Z. Chen, F. Zhang, M. Wei, and X. Kong, “Fluid flow in unconventional gas reservoirs,” Geofluids, vol. 2018, Article ID 2178582, 2 pages, 2018.
View at: Publisher Site | Google Scholar
J. L. Bischoff and K. S. Pitzer, “Phase relations and adiabats in boiling seafloor geothermal systems,” Earth and Planetary Science Letters, vol. 75, no. 4, pp. 327–338, 1985.
View at: Publisher Site | Google Scholar
G. C. Sills, S. J. Wheeler, S. D. Thomas, and T. N. Gardner, “Behaviour of offshore soils containing gas bubbles,” Geotechnique, vol. 41, no. 2, pp. 227–241, 1991.
View at: Publisher Site | Google Scholar
J. L. Grozic, P. K. Robertson, and N. R. Morgenstern, “The behavior of loose gassy sand,” Canadian Geotechnical Journal, vol. 36, no. 3, pp. 482–492, 1999.
View at: Publisher Site | Google Scholar
D. Ronen, B. Berkowitz, and M. Magaritz, “The development and influence of gas bubbles in phreatic aquifers under natural flow conditions,” Transport in Porous Media, vol. 4, no. 3, p. 295, 1989.
View at: Publisher Site | Google Scholar
J. K. Williams and R. A. Dawe, “Wetting phenomena and phase separation in porous media – qualitative features,” Journal of Colloid and Interface Science, vol. 125, no. 1, pp. 347–350, 1988.
View at: Publisher Site | Google Scholar
S. Leclaire, N. Pellerin, M. Reggio, and J. Y. Trépanier, “Multiphase flow modeling of spinodal decomposition based on the cascaded lattice Boltzmann method,” Physica A: Statistical Mechanics and its Applications, vol. 406, pp. 307–319, 2014.
View at: Publisher Site | Google Scholar
R. Shimizu and H. Tanaka, “Impact of complex topology of porous media on phase separation of binary mixtures,” Science Advances, vol. 3, article eaap9570, no. 12, 2017.
View at: Publisher Site | Google Scholar
S. Leclaire, A. Parmigiani, O. Malaspinas, B. Chopard, and J. Latt, “Generalized three-dimensional lattice Boltzmann color-gradient method for immiscible two-phase pore-scale imbibition and drainage in porous media,” Physical Review E, vol. 95, no. 3, 2017.
View at: Publisher Site | Google Scholar
B. Chopard, A. Dupuis, A. Masselot, and P. Luthi, “Cellular automata and lattice Boltzmann techniques: an approach to model and simulate complex systems,” Advances in Complex Systems, vol. 5, no. 02n03, pp. 103–246, 2002.
View at: Publisher Site | Google Scholar
X. Li, Y. Zhang, X. Wang, and W. Ge, “GPU-based numerical simulation of multi-phase flow in porous media using multiple-relaxation-time lattice Boltzmann method,” Chemical Engineering Science, vol. 102, pp. 209–219, 2013.
View at: Publisher Site | Google Scholar
H. Liu, Q. Kang, C. R. Leonardi et al., “Multiphase lattice Boltzmann simulations for porous media applications,” Computational Geosciences, vol. 20, no. 4, pp. 777–805, 2016.
View at: Publisher Site | Google Scholar
S. Leclaire, N. Pellerin, M. Reggio, and J.-Y. Trépanier, “Unsteady immiscible multiphase flow validation of a multiple-relaxation-time lattice Boltzmann method,” Journal of Physics A: Mathematical and Theoretical, vol. 47, no. 10, article 105501, 2014.
View at: Publisher Site | Google Scholar
S. Leclaire, N. Pellerin, M. Reggio, and J. Y. Trépanier, “A multiphase lattice Boltzmann method for simulating immiscible liquid-liquid interface dynamics,” Applied Mathematical Modelling, vol. 40, no. 13-14, pp. 6376–6394, 2016.
View at: Publisher Site | Google Scholar
S. Leclaire, A. Parmigiani, B. Chopard, and J. Latt, “Three-dimensional lattice Boltzmann method benchmarks between color-gradient and pseudo-potential immiscible multi-component models,” International Journal of Modern Physics C, vol. 28, no. 7, article 1750085, 2017.
View at: Publisher Site | Google Scholar
PALABOS, 2018, http://www.palabos.org.
S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond, Clarendon Press, 2001.
A. Parmigiani, C. Huber, O. Bachmann, and B. Chopard, “Pore-scale mass and reactant transport in multiphase porous media flows,” Journal of Fluid Mechanics, vol. 686, pp. 40–76, 2011.
View at: Publisher Site | Google Scholar
T. Akai, B. Bijeljic, and M. J. Blunt, “Wetting boundary condition for the color-gradient lattice Boltzmann method: validation with analytical and experimental data,” Advances in Water Resources, vol. 116, pp. 56–66, 2018.
View at: Publisher Site | Google Scholar
S. Saito, Y. Abe, and K. Koyama, “Lattice boltzmann modeling and simulation of liquid jet breakup,” Physical Review E, vol. 96, no. 1, article 013317, 2017.
View at: Publisher Site | Google Scholar
M. Avrami, “Kinetics of phase change. II Transformation‐time relations for random distribution of nuclei,” The Journal of Chemical Physics, vol. 8, no. 2, pp. 212–224, 1940.
View at: Publisher Site | Google Scholar
T. G. Hersum and B. D. Marsh, “Igneous microstructures from kinetic models of crystallization,” Journal of Volcanology and Geothermal Research, vol. 154, no. 1-2, pp. 34–47, 2006.
View at: Publisher Site | Google Scholar
E. S. Boek and M. Venturoli, “Lattice-Boltzmann studies of fluid flow in porous media with realistic rock geometries,” Computers & Mathematics with Applications, vol. 59, no. 7, pp. 2305–2314, 2010.
View at: Publisher Site | Google Scholar
P. R. Di Palma, N. Guyennon, F. Hesse, and E. Romano, “Porous media flux sensitivity to pore-scale geostatistics: a bottom-up approach,” Advances in Water Resources, vol. 102, pp. 99–110, 2017.
View at: Publisher Site | Google Scholar
V. Joekar-Niasar and S. M. Hassanizadeh, “Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: a review,” Critical Reviews in Environmental Science and Technology, vol. 42, no. 18, pp. 1895–1976, 2012.
View at: Publisher Site | Google Scholar
R. T. Armstrong, A. Georgiadis, H. Ott, D. Klemin, and S. Berg, “Critical capillary number: desaturation studied with fast x-ray computed microtomography,” Geophysical Research Letters, vol. 41, no. 1, pp. 55–60, 2014.
View at: Publisher Site | Google Scholar
W. R. Rossen, “A critical review of roof snap-off as a mechanism of steady-state foam generation in homogeneous porous media,” Colloids and Surfaces A: Physicochemical and Engineering Aspects, vol. 225, no. 1-3, pp. 1–24, 2003.
View at: Publisher Site | Google Scholar
W. Deng, M. Balhoff, and M. B. Cardenas, “Influence of dynamic factors on nonwetting fluid snap-off in pores,” Water Resources Research, vol. 51, no. 11, pp. 9182–9189, 2015.
View at: Publisher Site | Google Scholar
A. Parmigiani, S. Faroughi, C. Huber, O. Bachmann, and Y. Su, “Bubble accumulation and Its role in the evolution of magma reservoirs in the upper crust,” Nature, vol. 532, no. 7600, pp. 492–495, 2016.
View at: Publisher Site | Google Scholar
S. Arnórsson, A. Stefánsson, and J. O. Bjarnason, “Fluid-fluid interactions in geothermal systems,” Reviews in Mineralogy and Geochemistry, vol. 65, no. 1, pp. 259–312, 2007.
View at: Publisher Site | Google Scholar
S. Scott, I. Gunnarsson, S. Arnórsson, and A. Stefánsson, “Gas chemistry, boiling and phase segregation in a geothermal system, Hellisheidi, Iceland,” Geochimica et Cosmochimica Acta, vol. 124, pp. 170–189, 2014.
View at: Publisher Site | Google Scholar
S. Scott, T. Driesner, and P. Weis, “Boiling and condensation of saline geothermal fluids above magmatic intrusions,” Geophysical Research Letters, vol. 44, no. 4, pp. 1696–1705, 2017.
View at: Publisher Site | Google Scholar
J. L. Bischoff and R. J. Rosenbauer, “Phase separation in seafloor geothermal systems; an experimental study of the effects on metal transport,” American Journal of Science, vol. 287, no. 10, pp. 953–978, 1987.
View at: Publisher Site | Google Scholar
J. Cowan and J. Cann, “Supercritical two-phase separation of hydrothermal fluids in the troodos ophiolite,” Nature, vol. 333, no. 6170, pp. 259–261, 1988.
View at: Publisher Site | Google Scholar
D. I. Foustoukos and W. E. Seyfried, “Fluid phase separation processes in submarine hydrothermal systems,” Reviews in Mineralogy and Geochemistry, vol. 65, no. 1, pp. 213–239, 2007.
View at: Publisher Site | Google Scholar
F. K. Wilckens, E. P. Reeves, W. Bach et al., “The influence of magmatic fluids and phase separation on B systematics in submarine hydrothermal vent fluids from back-arc basins,” Geochimica et Cosmochimica Acta, vol. 232, pp. 140–162, 2018.
View at: Publisher Site | Google Scholar
T. W. Sisson and C. R. Bacon, “Gas-driven filter pressing in magmas,” Geology, vol. 27, no. 7, pp. 613–616, 1999.
View at: Publisher Site | Google Scholar
O. Bachmann and G. W. Bergantz, “Gas percolation in upper-crustal silicic crystal mushes as a mechanism for upward heat advection and rejuvenation of near-solidus magma bodies,” Journal of Volcanology and Geothermal Research, vol. 149, no. 1-2, pp. 85–102, 2006.
View at: Publisher Site | Google Scholar
M. Pistone, F. Arzilli, K. J. Dobson et al., “Gas-driven filter pressing in magmas: insights into in-situ melt segregation from crystal mushes,” Geology, vol. 43, no. 8, pp. 699–702, 2015.
View at: Publisher Site | Google Scholar
C. Huber and A. Parmigiani, “A physical model for three-phase compaction in silicic magma reservoirs,” Journal of Geophysical Research: Solid Earth, vol. 123, no. 4, pp. 2685–2705, 2018.
View at: Publisher Site | Google Scholar
A. K. Gunstensen, D. H. Rothman, S. Zaleski, and G. Zanetti, “Lattice Boltzmann model of immiscible fluids,” Physical Review A, vol. 43, no. 8, pp. 4320–4327, 1991.
View at: Publisher Site | Google Scholar
T. Reis and T. N. Phillips, “Lattice Boltzmann model for simulating immiscible two-phase flows,” Journal of Physics A: Mathematical and Theoretical, vol. 40, no. 14, pp. 4033–4053, 2007.
View at: Publisher Site | Google Scholar
H. Liu, A. J. Valocchi, and Q. Kang, “Three-dimensional lattice Boltzmann model for immiscible two-phase flow simulations,” Physical Review E, vol. 85, article 046309, no. 4, 2012.
View at: Publisher Site | Google Scholar
M. Latva-Kokko and D. H. Rothman, “Diffusion properties of gradient-based lattice Boltzmann models of immiscible fluids,” Physical Review E, vol. 71, article 056702, no. 5, 2005.
View at: Publisher Site | Google Scholar
S. Leclaire, M. Reggio, and J. Y. Trépanier, “Numerical evaluation of two recoloring operators for an immiscible two-phase flow lattice Boltzmann model,” Applied Mathematical Modelling, vol. 36, no. 5, pp. 2237–2252, 2012.
View at: Publisher Site | Google Scholar
S. Leclaire, N. Pellerin, M. Reggio, and J. Y. Trépanier, “An approach to control the spurious currents in a multiphase lattice Boltzmann method and to improve the implementation of initial condition,” International Journal for Numerical Methods in Fluids, vol. 77, no. 12, pp. 732–746, 2015.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2019 Andrea Parmigiani et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1163

Downloads

1481

Citations