ORIGINAL RESEARCH article

Front. Earth Sci., 27 June 2022

Sec. Structural Geology and Tectonics

Volume 10 - 2022 | https://doi.org/10.3389/feart.2022.947625

Discrete Bubble Flow in Granular Porous Media via Multiphase Computational Fluid Dynamic Simulation

  • 1. College of Environmental Science and Engineering, Dalian Maritime University, Dalian, China

  • 2. School of Civil Engineering, University of Queensland, Brisbane, QLD, Australia

Abstract

The coal seam gas industry has raised public concerns about the potential risk of groundwater contamination, where gases leaked from coal seams are thought to pollute groundwater. However, the basic principles and controlling parameters for gas seepage from deep ground formations to the ground surface have not been fully understood. As a possible mechanism for gas transport in the subsurface environment, discrete bubble flow was previously investigated using laboratory experiments by Ma et al. (Water Resour. Res, 2015, 51 (6), 4359–4373). This study developed a multiphase computational fluid dynamic (CFD) model to simulate discrete bubbly flow in a two-dimensional granular porous media at the pore scale. Following the experimental setup from Ma et al. (Water Resour. Res, 2015, 51 (6), 4359–4373), a “point source” with preset bubble fluxes was specified in a simulating domain representing the flume size in the earlier experiments. There were around 7,000 granular particles within this domain to model the porous media. This numerical model was validated by conserving the gas mass in the simulating domain. The simulation results provide more physical insights into complex bubble transport behaviour in porous media through specific plume parameters. The breakthrough time of the bubble plume and the cross-sectional averaged velocity of ambient pore water flow were manifested to be proportional to the gas release rates in the logarithmic scales. Also, the bubble plume width was also observed to be proportional to the gas release rates. Moreover, the gas distribution on the top boundary could be observed. The outcomes were further tested against the scaling solutions proposed by Ma et al. (Water Resour. Res, 2015, 51 (6), 4359–4373) with disagreements. The limitations of this multiphase computational fluid dynamic model were finally discussed.

Introduction

The coal seam gas industry has raised public concerns about the potentially high risk of potential groundwater pollution, where gases leaked from coal seam beds are prone to contaminate groundwater (). In addition, those light hydrocarbon gases in high concentrations, such as ethane, methane, propane, etc., can transport in hydrocarbon gas-storing rock mass, particularly in preferential flow paths comprised of fracture and fault zones (; ). Therefore, soil gas concentration anomalies in the topsoil are primarily attributable to those gas migration in the following mechanisms: 1) gas diffusion; 2) the dissolved gas carried by groundwater flow; 3) continuous gas phase driven by pressure gradient; and 4) bubbly flow transport (). Although many studies have comprehensively studied the former three mechanisms (, ; ), the basic principles and controlling parameters for bubbly flow from deep ground formations to the ground surface, proposed by , have not been comprehensively investigated and developed. This research objective is, so far, still at an initial stage with a few academic attempts (; ; ; ).

A conceptual model in Figure 1 illustrates subsurface microbubble transport processes, where bubbles were transported in the fracture, which acted as the point source of gas for the soil layer. The discrete bubble transport mechanism is not well understood, especially in the saturated shallow soils near the ground surface. This water-saturated porous media could act as a buffer zone for gas bubbles percolating through the basement rock and modifying bubble transport behaviour. Hence, soil gas concentration signals could be detected at the ground surface. has successfully carried out quasi-two-dimensional (2D) physical modelling to explore potential buffering effects to study microbubble migration and ambient groundwater circulation. A scaling solution, accounting for bubble injection rate, bubble plume tip velocity, width and breakthrough time of the plume front, was also established by . As a continuous work after this former study, we further performed 2D numerical modelling to simulate the transport process of discrete bubbles with a localized source of gas bubbles at the bottom boundary of the saturated soil later. This “point” source models the input condition provided by gas flow through fractures underlying the water-saturated soil layer. Based on these simulation outcomes, it is expected that more physical insights could be brought into complex bubble transport behaviour in porous media. Also, this work could adequately quantify both advective and dispersive bubble migrations for engineering practices.

FIGURE 1

Literature Review of Numerical Models at Multiscale

Numerical simulations of gas transport in porous media have been extensively studied previously. This section reviews two conventional gas transport models: the continuum model of two-phase flow in porous media () and the continuum “colloidal” discrete bubbles transport model () in order to identify their advantages and disadvantages. This review contributed to the development of a new model that could overcome those shortcomings.

The Model of Two-phase Flow in Porous Media

The two phases Darcy flow model treats both the water phase and gas phase as continuum phases (). Depending on how much space each phase occupies, the two phases interact dynamically by the relative permeability and effective saturation. The governing equations for water and gas flow are given ():

For the wetting phase (water):

For the nonwetting phase (gas):where subscript w is for the wetting phase (water phase); subscript nw is for the nonwetting phase (gas phase); z is the vertical coordinate, respectively; k is the media permeability [L2]; kr is relative permeability for the wetting or nonwetting phase [-]; μ is the dynamic viscosity [M/TL]; ρ is fluid density [M/L3]; g is the magnitude of gravitational acceleration [L/T2]; Se is effective saturation (varying from 0 to 1), and P is the pressure [M/LT2]. The capillary pressure could be expressed as Pc= Pnw–Pw.

Specific capacity is defined as:where C is the specific capacity [LT2/M] and is the porosity [L3/L3]. The , C, Se, k and Pc vary simultaneously. and applied the soil water retention and relative permeability functions to express changes in the , C, Se and k as functions of Pc. The capillary pressure was transformed to the equivalent water head as Hc= Pc/(ρwaterg).

The hydraulic properties relative to the wetting fluid are given below:

If Hc>0 (unsaturated zone):

If Hc≤0 (saturated zone):

For the nonwetting fluid:where ϕs is the total porosity; ϕr is the residual porosity to the total porosity; α, n, m and L are the parameters, which depend on the soil type. These parameters could be determined using unsaturated soil tests, including standard tests and soil column tests () with various unsaturated soil sensors ().

This approach described above can simulate continuous gas phase transport in the soil layer, driven by pressure. However, bubbles are driven by the buoyancy force and travel discretely in the porous media in terms of bubble migration. As a result, the gas relative permeability could not be applied to describe discontinuous bubble flow. Furthermore, this approach does not provide information on gas concentration, although bubble concentration is significant for interpreting the concentration mapping of soil gas. Therefore, different approaches to quantifying the microbubble flow are highly demanded and reviewed below.

The Continuum “Colloidal” Discrete Bubbles Transport

conducted qualitative experiments to study bubble flow in porous media, observing that large gas bubbles stuck to surfaces and were trapped at pore throats. These large bubbles then block the pathway of smaller rising bubbles. However, the bubbles considered in the present study are discrete bubbles in colloidal sizes. suggested that if bubbles are in colloidal sizes (in the range between 0.7 and 20 µm), Brownian motion will prevent them from sticking to the solid particle surfaces, and they could continue to move upward. tested this hypothesis using laboratory experiments to generate stable discrete bubbles and measure their transport properties in porous media. The results indicated that the surface interaction between sand grains and bubbles was not favourable for deposition. Under such surface interaction, bubbles can travel over significant distances in the subsurface environment. Furthermore, the individual gas bubbles can be seen as finite fluid particles within the liquid continuum. Therefore, microbubble transport in porous media may be treated as a particular case of colloid transport ():where C is the bubble concentration [M/L3]; R is the retardation factor for local sorption equilibrium; t is the time [T]; vo is the gas velocity [L/T]; D is the dispersion coefficient [L2/T] and is the trapping coefficient. As for applying the above colloidal transport equation, an important question arises regarding defining the dispersion coefficient and the microbubble flow velocity. Every single bubble has an individual dispersion ability and velocity. In an infinite viscous media without walls or other particles, the bubble movement in static water is assumed to follow Stoke’s law [22], which is applied to calculate the bubble velocity given bywhere ν is the bubble velocity relative to the water [L/T], d is the bubble diameter [L], ρw is the water density [M/L3], ρg is the gas density [M/L3], and µ is the dynamic viscosity of water [M/TL].

Stoke’s law is only applicable to small bubbles with an immobile surface. Eq 8 shows that the bubble velocity is directly related to the bubble diameter. The bubble radius grows when the hydrostatic pressure decreases, consequentially increasing bubble velocity. The bubble velocity given by Stokes’s law is only applicable to a single bubble motion. It does not consider interference by the other bubbles and the wall effect induced by the fracture/porous media.

When the bubbles come into close contact with walls, the terminal bubble velocities decrease. The viscous drag increases since the fluid, displaced by the rising bubble, is squeezed into a narrower cross-sectional area between the wall and particle, increasing the viscous resistance to flow. , who investigated the wall effects on the motion of a single particle, found that the wall effects are minor if the particle size is significantly smaller than the aperture. The bubbles velocity within a narrow fracture/pore space, including the wall effect, can be described by the following equation ():where r is the bubble radius [L], b is half the size of the fracture aperture [L], is the bubble velocity [L/T] along the fracture and is Stoke’s velocity without wall effect [L/T].

The interplay of bubble rebound and coalescence could be necessary when a group of bubbles migrate in the form of bubble trains. The bubble coalescence and rebound are controlled by three actions: attractive surface forces, hydrodynamic interaction, and dynamic bubble collisions (). When gas fluxes increase, bubbles can coalesce and produce vertically elongated bubbles, called “slugs”, creating continuous gas streams within the porous media (). In general, the maximum size of a bubble is constrained by the minimum cross-section area of the flow pathway, such as the throat of the fracture/granular porous media. established a model to estimate maximum bubble size as a function of the porosity and the grain size:where is the bubble diameter [L]; is the mean grain size [L], and is the media porosity.

This population colloidal microbubble transport model is efficient for modelling bubble flow in porous media. Nonetheless, this approach does not include the dynamic interaction between the bubble and ambient water flow. As prior experimental work by concluded, the ambient pore water flow may interact with bubble flow strongly.

Numerical Methodology

Multiphase Computational Fluid Dynamics

Considering the disadvantages of the aforementioned two models, we selected a two-phase bubbly flow model for this study. This pore-scale two-phase flow model described by the Navier-Stokes equation can simulate the buoyancy force and the dynamic interaction between bubbles and water. A numerical model using COMSOL Multiphysics engineering simulation software is developed. COMSOL is a finite element simulation software (). It originated from the MATLAB toolbox and has been well integrated with MATLAB, which provides convenience for pre-processing and post-processing. The program has been applied and tested extensively for a wide range of gas transport problems in saturated and vadose zones (; ). A pore-scale numerical model was conducted to simulate bubbly flow in a two-dimensional (2D) porous media, equivalent to the laboratory geophysical modelling published by . Compared to other two-phase flow models simulating the gas phase in continuous form (; ), this model can simulate discrete bubbles without coalescence.

Continuum models with mass and momentum conservation equations describe two-phase bubbly flow in porous media at the pore scale. Mass conservation (continuity) equations:

Momentum conservation equations (Navier–Stokes equation):where subscript denotes a liquid phase; subscript g denotes a gas phase; is the density; is the phase volume fraction; is the phase velocity; µ is the dynamic viscosity; is the phase pressure; is the local acceleration term; is the convective acceleration term; is the pressure gradient; is the gravitational force and is the viscosity force.

This numerical work follows the previous experimental setup from . The released bubbles remained spherical shape (the bubble size is around 300 µm) through the experiments and were about one order of magnitude smaller than hydrogel beads. In addition, no trapping and attachment of bubbles to solid matrixes were visually observed. Consequently, there is no direct contact between gas and solid matrixes throughout the experiments, leading to no capillary actions involving solid matrixes. Accordingly, capillary effects are not included in the bubble transport numerical model. Therefore, the Navier-Stokes equation is adequate and appropriate to simulate ambient pore water flow around discrete bubble transport. The discrete bubbles were simulated as spherical particles floating upward following a gas bubble model described by Eqs 1417.

According to the mass and momentum conservation equation, could be calculated aswhere is the velocity difference between the gas phase and liquid phase. The drag force exerted on each bubble is given by ).where is the drag force and is the diameter of the bubbles (; ; ). The Hadamard-Rybczynski formula (; ) can be applied to calculate :

The gas density is calculated from the ideal gas laws:where p is the pressure of the gas, M is the amount of substance of gas, R is the universal gas constant, and T is the absolute temperature of the gas. The relationship between the gas volume fraction and liquid volume fraction is:

Numerical Simulation Setup and Settings

As shown in Figure 2, a two-dimensional water-saturated domain was created with the dimensions of 700×900 mm2 (width height), the same size as the experimental flume from . There are approximately 7,000 spherical particles preset in the domain in order to simulate the transparent granular soil as a porous media (), which yields a porosity of 0.494. All solid particles were not movable and rotatable. Then, a sequential sphere packing method was applied to establish the pore-scale domain. This packing method employs trilateration equations to pack spheres with a predefined grain size distribution (GSD of hydrogel beads) and a given porosity (). In addition, the boundary conditions are specified, as shown in Figure 2.

FIGURE 2

A series of numerical experiments with different gas release rates of 5, 10, 20 and 30 ml/min were conducted. The gas inlet diffuser is a gas flux boundary. Thus, the gas release rates were converted into the gas flux of , , and kg/(m2·s) corresponding to the gas release rates and bubble diffuser setting applied in the laboratory experiments by .

Results and Discussion

Validation of Mass Conservation

In order to test the model performance, mass conservation is examined. As shown in Figure 3, the total gas flux at the top boundary of the domain was plotted over time with different gas release rates. The evolution of the gas flux over time exhibits a rapid increase and sharp overshoots in the initial stage. After the gas breakthrough, gas fluxes were stabilized and kept at a plateau under steady-state conditions. The overshoot could be caused by massive ambient gas transporting together with a bubble plume during a short time right after the gas breakthrough, but this local disturbance quickly fades away. Although the input gas fluxes are labelled for different gas release rates on each curve, they all match well with gas output fluxes at the steady-state conditions. Therefore, the model is manifested to be mass conservative.

FIGURE 3

Bubble Plume Width Development

In order to examine results simulated by the numerical models in comparison with the laboratory experiments carried out by , the simulation results are analyzed following the procedure for analyzing the experimental data in the previous work (). The width variations of the estimated bubble plume at the top of the domain against gas release rates are shown in Figure 3A after the bubble plume breakthrough and reaching the steady-state conditions. The slope of the trendline is 1.65, which is around 1.6 times that of the laboratory trendline (slope = 1). Figure 4.

FIGURE 4

Gas flux distributions at the top boundary of the domain are plotted to understand this inconsistency. Figure 3B shows that the gas flux variation at the top boundary corresponds with the increasing gas release rate. The modelled is narrower than the laboratory results provided by for a low gas release rate of 5 ml/min. However, for a high gas release rate of 30 ml/min, the modelled is wider than the laboratory results, which leads to a steeper trendline slope. These discrepancies could be attributed to differences between the 2D granular porous media in this simulation and the quasi-2D (actual 3D) transparent soil in the previous experiment by .

Gas Distribution on the Top Boundary

As previously discussed, there are 11 divisions in the upper part of the laboratory flume. Therefore, in order to compare against the laboratory results, modelled gas fluxes at the top boundary were assigned and normalized into 11 divisions corresponding to 11 cells at the top of the laboratory flume. As observed during the laboratory experiment, the gas discharge from the top of the domain concentrated over the plume area with only three middle cells receiving gas. The maximum rate occurred at the centre above the gas release point. However, the distributions are flatter regarding the modelled gas flux variations and exhibit a symmetrical spreading. Additionally, according to laboratory results, the gas discharge rate decreased with the gas release rate in the centre cell. Still, it increased in the neighbouring cells on both sides, whereas modelling results are more random without a consistent tendency (see Figure 5).

FIGURE 5

Bubble Plume Fronts Evolution

The evolution of the plume front was defined by a fixed gas volume fraction, which is equal to 20% of the maximum gas volume fraction. Note that the gas volume fraction is defined within a Representative Elementary Volume (REV) of mm2, consistent with the REV applied for analyzing the laboratory data by . As shown in Figure 6, Ht is the distance of the plume front from the release point before the bubble breakthrough. The evolution of the plume front over time exhibits an advancing speed that is analogous to 5 and 10 ml/min. The advancing speed also behaved similarly for the gas release rates of 20 and 30 ml/min, but both slightly increased compared with the results for the low gas release rates. Unlike the laboratory results given by , there is no consistent advancing speed among the results for the various gas release rates.

FIGURE 6

The plume breakthrough times were obtained for all cases with different gas release rates by tracking the bubble plume front. In order to test if the numerical results also follow the scaling solutions proposed for the laboratory experiments conditions by , versus was plotted. The slope for the trendline between and is -0.137, which is three times smaller than that of the scaling solution proposed for the laboratory experiment (–0.4) by . The overall is around 3–5 s smaller than the laboratory results of , as shown in Figure 7.

FIGURE 7

Bubble Plume Lateral Width Evolution

The mean lateral width of the bubble plume was tracked over time, as shown in Figure 8. The lateral width of the bubble plume gradually increased with the gas release rate. At the initial stage of the plume development (less than 15 s before the bubble plume front breakthrough), the higher gas release rates produced wider plume widths. At the later stage of the plume development, the plume width approached steady-state conditions with time. The temporal variation of the plume width standard deviation exhibited a gradual increase followed by a relatively constant level, especially after the plume breakthrough. Attempts were also given to apply a characteristic width to normalize/non-dimensionalize the numerical results. However, it failed to find a suitable characteristic width to converge all curves into a relatively uniform trend.

FIGURE 8

Bubble Plume-Induced Pore Water Movement

The numerical model also demonstrated pore water flows associated with gas ebullition. Circulation was clearly illustrated due to the lift-up of pore water inside the bubble plume and the draw-down of pore water outside the plume. As shown in Figure 9, the vertical pore water flow velocities were plotted along transects at various heights of the simulation domain. Pore water was also descending beside the ascending pore water in the domain centre (more likely within the bubble plume). This descending pore water could be caused by localized circulations inside the plume, which may confine the bubble plume size. Outside the bubble plume, almost all of the water flows downward, and the higher gas release rate leads to a faster pore water velocity.

FIGURE 9

The Inconsistency With Prior Experimental Outcomes

The average pore water velocity over the cross-sectional area outside the bubble plume was plotted versus the gas release rate in Figure 10. Based on the scaling solution described by , is proportional to . The numerical results were tested by plotting a trendline with an exponent of 0.6. Unfortunately, the slope of 0.6 underestimated the actual trend, which is better fitted by an incline of 0.8.

FIGURE 10

Compared with laboratory conditions from , the inconsistent modelled results of the bubble plume width and the breakthrough time could be influenced by the drag coefficients applied in the model. A simple numerical experiment was conducted to examine the influence on the bubble plume width. The drag coefficient in Figure 11B is 100 times larger than that of the drag coefficient in Figure 11A. Consequently, the lateral bubble plume width is almost doubled in Figure 11B. As a result, we exert more resistance force on bubbles by applying a more significant drag coefficient. Thus, the vertical velocity of the bubble plume would be reduced, prolonging the bubble plume breakthrough time and enhancing bubble dispersion. Those effects may lead to more realistic results than the laboratory data from . The bubble plume breakthrough time modelled in this work is less than the breakthrough time observed in the laboratory experiments from by approximately 3–5 s.

FIGURE 11

Further investigations are required to better understand the impacts of drag coefficients on the bubble transporting behaviour and ambient pore water flow and develop a proper drag coefficient function to replicate the laboratory experiments. A parametric and sensitivity study might also need to be carried out for those purposes in the future. However, this does not undermine the value of this numerical modelling work as the first attempt because it is still a good initiation in addition to the prior physical modelling.

Conclusion and Reflections

A pore-scale two-phase bubbly flow model was developed and tested in this study. The numerical model was applied to replicate the laboratory experiments and allow more detailed studies, which were not feasible in the laboratory experiments. The preliminary numerical results have revealed that the proposed continuum numerical method can simulate bubbly flow events and bubble-water dynamic interaction phenomena. This numerical model was validated by conserving the gas mass in the simulating domain. The simulation results provide more physical insights into complex bubble transport behaviour in porous media through specific plume parameters. The breakthrough time of the bubble plume and the cross-sectional averaged velocity of ambient pore water flow were manifested to be proportional to the gas release rates in the logarithmic scales. Also, the bubble plume width was also observed to be proportional to the gas release rates. Moreover, the gas distribution on the top boundary could be observed.

This study examined the numerical outcomes against laboratory results and the scaling solution described in our previous experimental works. Although the modelling results do not match the laboratory data quantitatively, the modelled results capture bubble plumes’ trends/tendencies/features under various gas release rates, such as the plume width, plume breakthrough time, and pore water flow velocity.

The proposed numerical method has a relatively high computing cost, with around one million elements needed to represent the domain. This requires approximately 1 week to complete one gas release rate simulation case. As it was demonstrated that the drag coefficient applied in the model impacts the bubble transport behaviour, further studies should be carried out to focus on the drag coefficient effects. In addition, the solid particles packing method should also be investigated further in the future.

Scaling solutions proposed in the prior experimental study may only be applicable for discrete bubble transport in porous media with a specific range of porosity and permeability. Therefore, further studies are required to examine the hydrodynamic interactions between bubble ebullition and pore water under different boundary conditions within different types of porous media. Furthermore, the effects of the water table on microbubble transport behaviour could also be studied.

Bubble transport in the porous media has been studied with a minimal ratio of bubble size to solid particle size without observing bubble trapping. By increasing this ratio, the porous media could trap the bubbles in the pore space or pore throat. This bubble trapping would certainly change bubble transport behaviour and the hydrodynamic interactions between bubbles and ambient pore water. Further studies are required to investigate various ratios between bubbles and solid particle sizes. Thus a bubble phase diagram could be produced to determine the correlation between the bubbles’ trapping and velocity.

Statements

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

Conceptual design: YM and GY; Modelling method: YM; Writing: YM and GY; Editing: GY; Supervision: AS.

Funding

The work was funded by the National Key Research and Development Program of China (2021YFB2601100) and Fundamental Research Funds for DMU 3132022168. This work was also funded by the Linkage Project (LP120100662), Multiscale two-phase flow in complex coal seam systems, funded by the Australian Research Council and supported by a Queensland Science Fellowship awarded to AS.

Acknowledgments

The authors would like to acknowledge support from the Geotechnical Engineering Centre (GEC), School of Civil Engineering, the University of Queensland, Australia.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Summary

Keywords

multiphase computational fluid dynamic, bubble plumes migration, pore water flow, subsurface gas distribution, porous media

Citation

Ma Y, Yan G and Scheuermann A (2022) Discrete Bubble Flow in Granular Porous Media via Multiphase Computational Fluid Dynamic Simulation. Front. Earth Sci. 10:947625. doi: 10.3389/feart.2022.947625

Received

18 May 2022

Accepted

07 June 2022

Published

27 June 2022

Volume

10 - 2022

Edited by

Xiaofei Jing, Chongqing University of Science and Technology, China

Reviewed by

Pei Zhang, Westlake University, China

Tao Wang, Hohai University, China

Updates

Copyright

*Correspondence: Ye Ma, ; Guanxi Yan,

† ORCID: Ye Ma, orcid.org/0000-0001-5999-4650; Guanxi Yan, orcid.org/0000-0002-1466-7805

This article was submitted to Structural Geology and Tectonics, a section of the journal Frontiers in Earth Science

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics