Suitable Lattice Boltzmann Method for Hydrogen-Brine Multiphase flow?

Suresh Ahuja

Coupled Lattice Boltzmann Modeling Framework for Pore-Scale Fluid Flow and Reactive Transport

Siyan Liu
,
Reza Barati*
,
Chi Zhang
, and
Mohammad Kazemi

Cite this: ACS Omega 2023, 8, 15, 13649–13669

Publication Date:April 3, 2023

https://doi.org/10.1021/acsomega.2c07643

CC-BY-NC-ND 4.0.

📷

Article Views

1104

Altmetric

Citations

LEARN ABOUT THESE METRICS

Add to

ExportRIS

PDF (10 MB)

📷

ACS Omega

Abstract

📷

In this paper, we propose a modeling framework for pore-scale fluid flow and reactive transport based on a coupled lattice Boltzmann model (LBM). We develop a modeling interface to integrate the LBM modeling code parallel lattice Boltzmann solver and the PHREEQC reaction solver using multiple flow and reaction cell mapping schemes. The major advantage of the proposed workflow is the high modeling flexibility obtained by coupling the geochemical model with the LBM fluid flow model. Consequently, the model is capable of executing one or more complex reactions within desired cells while preserving the high data communication efficiency between the two codes. Meanwhile, the developed mapping mechanism enables the flow, diffusion, and reactions in complex pore-scale geometries. We validate the coupled code in a series of benchmark numerical experiments, including 2D single-phase Poiseuille flow and diffusion, 2D reactive transport with calcite dissolution, as well as surface complexation reactions. The simulation results show good agreement with analytical solutions, experimental data, and multiple other simulation codes. In addition, we design an AI-based optimization workflow and implement it on the surface complexation model to enable increased capacity of the coupled modeling framework. Compared to the manual tuning results proposed in the literature, our workflow demonstrates fast and reliable model optimization results without incorporating pre-existing domain knowledge.

This publication is licensed under

CC-BY-NC-ND 4.0.

📷
📷
📷
📷

1. Introduction

ARTICLE SECTIONS

Jump To

Pore-scale numerical modeling for flow and reactive transport is crucial to understand the underlying mechanisms of physical-chemical interactions of numerous subsurface natural processes. Advanced waterflooding enhanced oil recovery (EOR) (1) and CO2 injection and sequestration (2) are two typical subsurface reactive transport-associated scientific and engineering problems that involve complex physio-chemical processes. Meanwhile, they are by nature highly uncertain processes due to the high complexity level of geological settings, which brings more challenges to conducting high-resolution mechanistic modeling studies compared to the reactive transport problems that happen in an artificial structure or domain such as reaction beds or combustion engine cylinders. Take advanced waterflooding or low-salinity waterflooding EOR modeling as an example, multiple underlying mechanisms are potentially working together and affect the rock/mineral surface wettability, and hence improve the recovery, but the actual underlying mechanisms are not fully understood and still debatable. (3−5) Meanwhile, the modeling approaches to simulate this process have simplified the problems without incorporating all physio-chemical processes in the pore-scale high-resolution fashion. (6−9)

To investigate the reactive transport and the relevant multiphysics phenomenon associated with the abovementioned scientific and engineering problems, a modeling framework capable of handling high-resolution pore-scale single-/multiphase fluid flow with geochemical reactions and other known or potential new physics is crucial for achieving the desired research goals. The base of this type of modeling framework is the fluid flow model, which preferably can handle complex flow channel geometries, run in parallel CPUs/GPUs for performance, and incorporate other physio-chemical models.

Besides the conventional computational fluid dynamics (CFD) methods such as finite difference method, finite volume method (FVM), and finite element method (FEM), lattice Boltzmann method (LBM) is a unique modeling framework that is constructed with particle collision and streaming in regular lattice grids and solves macroscopic equations for the mesoscopic fluid flow kinetics. (10) LBM is proved to be a flexible and practical framework for modeling pore-scale fluid flow and multiphysics behaviors. (11−13) In addition to the single-phase fluid flow model, multicomponent multiphase LBM models are attributed to another major component of the LBM family as the equivalent of other multiphase CFD methods, for instance, the level set method, (14) front-tracking method, (15) and volume-of-fluid method. (16) The most popular multiphase LBM models used in this community are the Shan–Chen model, (11,17) the color-gradient model, (18) the free-energy model, (19,20) and the He–Chen–Zhang model. (21)

Regarding the pore-scale advection–diffusion and reactions, the LBM framework is flexible to adapt the conventional advection–diffusion equations and place the species concentration in scalar fields coupled with a fluid flow lattice. (22) Various studies have focused on the LBM reactive transport modeling, for example, LBM chemical transport, (23−25) reactive transport with mineral reactions, (25) mineral dissolution and dissolution-induced rock mechanical property changes, (26,27) CO2 reactive transport, (28) reactions between immiscible fluids, (29) and convection and heat transfer-associated problems. (30,31) More studies for the LBM multicomponent reactive transport models can be found in the review by Yoon et al. (32)

On the other hand, to enhance the flexibility and capabilities of LBM models in solving complex fluid flow and reactive transport problems, numerous studies have focused on coupling other or external models to different applications. For instance, the proposed applications include but are not limited to the coupled models that incorporate FVM with LBM; (33−35) multiscale coupling for heat, mass transfer, and reactions using a general reconstruction operator; (36) mass transport; (37) mineral dissolution and phase transition-associated multiphase reactive transport; (25,38−41) and hydrodynamic property evolution with mineral dissolution. (40,42)

Among those external solvers, PHREEQC (43) developed by USGS and its variants IPhreeqc (44) and PHREEQCRM (45) have been used extensively for reactive transport coupling due to their strong capabilities for solving complex reactions. (46) IPhreeqc utilizes the Microsoft component object model (COM) interface to communicate with the outside code. Numerous studies have used IPhreeqc for reactive transport modeling, including LBM-based models. (47−56) However, the data communication efficiency is limited for large-scale simulations. PHREEQCRM works as a C++ library at a lower level with a significantly improved data exchange speed, which is suitable for reactive transport coupling. There are a few notable applications coupled with PHREEQCRM used in reactive transport modeling. (57−60)

In recent years, with fast-paced progress in data analytics, artificial intelligence (AI), and machine learning (ML), AI has been extensively implemented in solving various scientific and engineering problems. For subsurface pore-scale modeling specifically, AI/ML models are used for different applications, for example, microstructure synthesis and pore space reconstruction using generative adversarial network-based models, (61−66) image segmentation for rock and flow, (67−70) petrophysical property and other physical property estimation, (71−79) prediction of fluid flow, (80) and effective reaction rates. (81) More applications of ML and geoscience-associated studies can be found in the thorough review studies by Tahmasebi et al. (82) and Wang et al. (83) According to their comprehensive review and many other studies, advanced AI/ML models have shown great capabilities in analyzing big data, complex pattern recognition and extraction, and knowledge discovery for physical science problems. AI/ML models are known for building surrogate models for property estimation. Meanwhile, the AI/ML models also demonstrate strong capabilities in assisting optimization of the numerical models.

However, the existing literature in pore-scale reactive transport usually suggests that, in the presence of multiple physical-chemical processes in porous media, the current models lack flexibility or extensibility to solve such complex conditions. For instance, in the advanced waterflooding EOR process, multiple reactions, including but not limited to equilibrium, kinetic, and surface complexation (SCM) reactions, happened simultaneously, which causes dynamic changing fluid–fluid, fluid–solid interaction, and electrokinetics on the rock surface. To comprehensively study and model such complex processes, it is critical to develop an effective framework to properly incorporate the aforementioned factors in the modeling process. In addition, the developed model should be able to scale up and go beyond the demonstration purposes in a parallel computing scheme to solve real-world problems. Furthermore, numerical fluid flow and reactive transport models usually involve tedious and time-consuming parameter tuning, which is computationally expensive. A good integrated model/parameter optimization tool that does not require too much human intervention would be greatly beneficial for obtaining optimum model and model parameters.

In this work, to address the issues mentioned above, a coupled modeling framework based on LBM fluid flow and advection–diffusion models using the PALABOS package is proposed. (84) The framework incorporates the geochemical reaction solver PHREEQC (85) to handle complex reactions within the porous medium. To gain access to the lower level of the reaction solver for better computational efficiency, the PHREEQC variant PHREEQCRM (45) is selected in the coupling model, and it works as a C++ library that can be compiled together with the LBM flow solver. A coupling interface is developed to connect the two solvers by adding other capabilities, for example, the cell mapping system that handles different physical-chemical property fields within the flow domain. The framework enables a flexible setup for various chemical reactions both in the aqueous phase and on a solid surface and utilizes the unique capabilities of the two codes as much as possible. In addition, an AI-assisted automatic model optimization workflow is developed that integrates and works with the coupled numerical model. The optimization method is applied to the SCM adjustment for investigating low salinity waterflooding mechanisms. In this paper, Section 2 (Methodologies) includes the detailed model descriptions and coupling mechanisms and the AI-assisted optimization workflow. Four model validation cases in Section 3 (Model Validation and Results) are presented, including simple fluid flow, diffusion validation against the analytical solutions, reactive transport benchmark with calcite dissolution, and AI-optimized SCM results. Section 4 provides conclusions of the work and future research and development directions.

2. Methodology

ARTICLE SECTIONS

Jump To

2.1. Overview of the Proposed Modeling Workflow

To investigate the pore-scale fluid flow and physiochemical processes effectively and comprehensively, a numerical modeling framework that couples fluid flow models (LBM-based) and geochemical reaction solvers (PHREEQCRM) is proposed in this study. A tailor-made coupling interface connects the two major components for data communications.

The major development goal of the coupled modeling framework is to efficiently and effectively study various pore-scale reactive transport-associated scientific problems, while preserving the flexibility of the future extension in terms of multiphase flow and new add-on physical or chemical models. Figure 1 illustrates the overall workflow of the modeling framework.

Figure 1 📷Figure 1. Schematic plot of the modeling framework. The framework has six major components, including initialization of the domain, fluid, and reactants based on both solvers; advection–diffusion process handled by PALABOS; coupling interface that connects the two solvers; chemical reactions handled by PHREEQCRM; electrical property calculations using PHREEQCRM; and model evaluations and visualizations. The simulation integrates all components dynamically in a loop shown in the figure.

2.2. LBM Fluid Flow Model

2.2.1. LBM Fundamentals

LBM represents a generic CFD modeling framework. LBM is not only capable of modeling single-phase and multiphase fluid flow but also of building complex flow-based multiphysics systems such as the advection–diffusion–reaction system and complex fluid–fluid and fluid–solid interaction scenarios with high parallel computational efficiency. (10)

One of the most widely used LBM models was proposed by Bhatnagar, Gross, and Krook (BGK) (86) with the approximation of the non-linear collisional dynamics by introducing the collision operator as

Ω𝑘=−1𝜏(𝑓𝑖−𝑓EQ𝑖)

where τ stands for the relaxation time used to reach the local equilibrium and fiEQ represents the particle distribution function at equilibrium. By assuming the existence of a local attractor and considering its impact on the collision process, the BGK model linearizes the collision term to reach local equilibrium using a fixed relaxation time and fluid viscosity. The BGK model utilizes single relaxation time τ in the evolution function (eq 9). Future research and model development will include the two-relaxation time or multi-relaxation time models for improved model stability purpose.(1)

Fluid flow macroscopic quantities are calculated in the BGK model based on the instantaneous status of particle distribution function within a specific grid block. For example, in a D2Q9 model (2D and 9 velocities), the lattice density ρLB(x,t) is obtained by summing the particle distribution function on all points within the given grid block

𝜌LB(𝒙,𝑡)=∑𝑖=08𝑓𝑖(𝒙,𝑡)

(2)

where fi(x,t) represents the particle distribution function for velocity i. Meanwhile, the lattice pressure PLB is recovered by

𝑃LB(𝑥,𝑡)=𝐶2s𝜌LB(𝑥,𝑡)

(3)

where Cs is a constant value for the lattice speed of sound, (87) which varies for different models. For instance, in the D2Q9 model, as mentioned, 𝐶s=1/3⎯⎯√, and the lattice velocity vector uLB(x,t) is obtained by

𝒖𝐋𝐁(𝒙,𝑡)=1𝜌LB(𝒙,𝑡)∑𝑖=08𝑓𝑖(𝒙,𝑡)𝒆𝒊

where ei represents the relative locations of the points within the specific grid block pointing to the adjacent grids, whose values are 0, 1, or −1; i is the index for the points within the grid block for the given model; x stands for the location of the grid block in the 2D or 3D simulation domain; and t describes the time step in the simulation for time-dependent particle distribution function evolution.(4)

The particle distribution function for the equilibrium status fieq is given based on the two local macroscopic lattice properties: lattice velocity u(x,t) and lattice density ρLB(x,t). To ensure the mass and momentum conservation during the collision and streaming process, specific equations and weights ωi are required for different models. (88,89) In the D2Q9 model, the equilibrium distribution functions are written as

𝑓eq0(𝒙,𝑡)=𝜔0𝜌LB(𝒙,𝑡)[1−32𝒖𝐋𝐁(𝒙,𝑡)2]

(5)

𝑓eq𝑖(𝒙,𝑡)=𝜔𝑖𝜌LB(𝒙,𝑡)[1+3𝒆𝒊·𝒖𝐋𝐁(𝒙,𝑡)+32[3(𝒆𝒊·𝒖𝐋𝐁(𝒙,𝑡))2−𝒖𝐋𝐁(𝒙,𝑡)2]]

where ωi represents a series of weights of a given lattice model to maintain the isotropy of the fourth-order tensor of velocities and the Galilean invariance. (88) The weight distribution of the points within a single lattice grid for three widely used lattice models is described in Table 1.(6)

Table 1. Weights ωi for Lattice Models D1Q3, D2Q9, and D3Q19

Modelcenterhorizontal/verticaldiagonalD1Q32/31/60D2Q94/91/91/36D3Q191/31/181/36

Take the D2Q9 model used in this work as an example; ω0 = 4/9 in eq 5 describes the equilibrium function for the center point within the grid block, where the rest of the 8 points are given by eq 6 with ωi = 1/36 or 1/9 depending on the geometrical position of the point with respect to the grid block.

Figure 2 demonstrates the lattice arrangement and velocity vectors for a D2Q9 model. The center point is labeled as 0, and indexes 1, 2, 3, and 4 indicate the velocity directions to the center of four sides, where 5, 6, 7, and 8 point to the four corners of the lattice grid.

Figure 2 📷Figure 2. Lattice arrangement and velocity vectors for the D2Q9 model.

The ei is determined by

📷

(7)

Equation 7 is converted to the following form

𝑒𝑖𝑒0=(0,0)𝑒1,3,𝑒2,4𝑒5,6,7,8==(±1,0),(0,±1)(±1,±1)

(8)

The overall BGK LBM model is given by the following evolution equation

𝑓𝑖(𝒙+𝒆𝒊,𝑡+1)=𝑓𝑖(𝒙,𝑡)−Ω𝜏[𝑓𝑖(𝒙,𝑡)−𝑓eq𝑖[𝜌LB(𝒙,𝑡),𝒖𝐋𝐁(𝒙,𝑡)]]

(9)

where Ωτ is the relaxation frequency defined by Ωτ = 1/τ; the relaxation time τ indicates the time used for the system to reach equilibrium during each streaming-collision cycle, which is associated with the lattice kinematic viscosity νLB (88)

𝜈LB=(𝜏−1/2)𝐶2s

(10)

A typical iteration cycle of an LBM BGK model starts with the initialization of the lattice velocity and lattice density to an equilibrium state based on the initial condition. Furthermore, the macroscopic quantities (e.g., lattice pressure, density, velocity) are obtained for the initial equilibrium state using eqs 2–4. Meanwhile, the calculated quantities for the current time step are used for computing the equilibrium function fieq based on eqs 5 and 6 (for the D2Q9 model in this case). Next, the collision and streaming processes based on the evolution eq 9 are executed. Specifically, the local particle distribution function fi(x,t) is altered by the current equilibrium function fieq within the relaxation time τ described in the evolution function. Followed by the collision step, the streaming step simply transfers the updated distribution functions to adjacent nodes, and a full iteration is completed. Boundary treatments are executed after the full cycle if applicable.

The general steps of the BGK LBM model can be described in Figure 3.

Figure 3 📷Figure 3. Schematic plot for the evolution of the BGK LBM model within a single time step. The orange box demonstrates the optional step, whereas the blue box indicates the necessary step during the evolution, modified from ref (90).

2.2.2. LBM Advection–Diffusion Model

A series of advection–diffusion types of problems can be formed within the LBM modeling framework, for instance, chemical transport, (24,25) heat transfer, (30,91) multicomponent reactive transport in porous media, (25) and reactive transport associated solid dissolution modeling. (92) Some of the studies coupled the LBM flow model with a third-party reaction solver to model the complex pore-scale reactive transport scenarios. (50,51,93−96) More related studies on LBM reactive transport-associated applications can be found in review papers. (46,97)

A typical advection–diffusion formulation is shown below

∂𝐶∂𝑡+∇·(𝐶𝒖)=∇·(𝐷∇𝐶)+𝑞

(11)

where C represents a scalar field for quantities that can be transported by the flow while diffuse based on the concentration or temperature gradient, u stands for the flow or advection velocity vector, D is the diffusion coefficient, and q is an optional source term which can be used to represent the locally destroyed or produced chemical species (if C is for chemical concentration field) or represents the consumption or production of heat (if C is the temperature field). These behaviors can be attributed to the mechanisms such as chemical reactions. In terms of the diffusion process, the associated LBM formulation can be written as the following evolution function

𝑓AD𝑖(𝑥+𝑐𝑖Δ𝑡,𝑡+Δ𝑡)−𝑓AD𝑖(𝑥,𝑡)=Ω𝑖(𝑥,𝑡)+𝑄𝑖(𝑥,𝑡)

(12)

The Qi(x,t) represents the same meaning as the source term q in eq 11 in a “discretized” form with respective to lattice space and time. For example, if we have an additional reaction solver to “produce” or “destroy” chemical species in each lattice grid at each lattice time step, this would be added in the “Q” source term of the LBM evolution equations to incorporate the advection–diffusion process.

The relaxation frequency is calculated as

Ω𝑖(𝑥,𝑡)=−1𝜏fAD[𝑓AD𝑖(𝑥,𝑡)−𝑓ADeq𝑖(𝑥,𝑡)]

(13)

The concentration or temperature term C is recovered by the equation below similar to the lattice density ρ in the BGK LBM flow model

𝐶=∑𝑖𝑓AD𝑖

(14)

Meanwhile, the diffusion coefficient D is obtained by the following equation in light of the lattice kinematic viscosity (ν) in the BGK LBM flow model

𝐷=𝐶2s(𝜏fAD−Δ𝑡2)

(15)

The associated lattice equilibrium distribution function is written as

𝑓ADeq𝑖=𝑤𝑖𝐶[1+𝑐𝑖·𝒖𝐶2s+(𝑐𝑖·𝒖)2𝐶4s−𝒖·𝒖2𝐶2s]

(16)

The advection and diffusion coupling can be achieved by explicitly linking the concentration or temperature field C with another fluid flow lattice with a unified solving time step.

2.2.3. PALABOS LBM Solver

LBM is chosen as the fluid-solving model due to its advantages, including but not limited to parallel, coupling, and complex boundary treating capabilities. The PALABOS (parallel lattice Boltzmann solver) package is selected as the base of the modeling framework for the fluid flow component.

PALABOS (98) is an open-source LBM-based fluid flow solver built with the C++ programming language, developed and maintained by the University of Geneva. It is capable of solving single/multiple-phase fluid flow and convection–diffusion type problems in open or closed geometry. It also features parallel computing capabilities by utilizing a few to thousands of CPUs to speed up the simulations. Meanwhile, PALABOS provides flexible application programming interfaces (APIs) to extend the existing capabilities, such as modeling fluid–solid interactions by coupling a large-scale atomic/molecular massively parallel simulator, (99) or simulating blood cell behaviors by coupling structural solver npFEM. (100,101)

Regarding the extensive APIs, PALABOS supports convenient non-local operations by “data processors” in addition to in-place LBM particle collision and streaming operations. Data processors enable robust outside algorithm implementation on either lattices or additional layers of data fields. Examples of this type of implementation include advection–diffusion coupling, multiphase flow, and solid–fluid interactions.

2.3. Geochemical Reaction Model

2.3.1. Geochemical Reaction Solver

Similar to the selection of the fluid flow solver, the widely used open-source geochemical reaction solver PHREEQC (pH-Redox-EQuilibrium) (102) developed and maintained by the United States Geological Survey (USGS) is selected in this work as the geochemical solver. PHREEQC is a scripting input-based reaction solver written in C and C++ languages for equilibrium, kinetic, ion-exchange, surface complexation reactions, solid solutions, and 1D reactive-transport solving capabilities.

In addition to the original version of PHREEQC, which is equipped with a graphical user interface, other variants are developed to be able to couple with thirty-party codes. For example, IPhreeqc (44) utilizes the Microsoft COM interface for communicating with outside models, which is widely used in modeling reactive transport. (47,51,54,55) Another variant PHREEQCRM (45) is used in this study, enabling most of the low-level APIs to be exposed to developers for deep code integrations while preserving a flexible scripting input system. Meanwhile, the reactions can be solved in parallel by utilizing mature interfaces such as OpenMP and MPI schemes. The implementation of PHREEQCRM library includes FEFLOW, (59) USGS PHAST, (57,58) and modeling ionic transport. (60)

2.3.2. PHREEQCRM Geochemical Reaction Models

One of the appealing features of PHREEQCRM is the “batch reactor” functionalities, which allow one or multiple geochemical reactions to be calculated within a single grid block at a time. Also, it is possible to enable different types of reactions flexibly by using their scripting input system. The available reaction types (keywords) include: (1) solutions, (2) Equilibrium_Phases, (3) exchange, (4) surface, (5) Gas_Phase, (6) Solid_Solutions, and (7) kinetics. Before executing different reactions, PHREEQC takes user-defined reactions and initializes the “batch reactor” to equilibrium with the initial chemical components. (45) In this study, three reaction types were considered: solutions, kinetics, and surface for aqueous phase equilibrium reaction, kinetic reactions, and surface complexation reactions (SCM) on the solid surface.

A discretized 2D/3D simulation domain contains thousands or millions of grid blocks that are considered parallel batch reactors in PHREEQCRM. They are independent of one another. The transport of chemical species is executed by a third-party flow model (PALABOS LBM model in this study) and moves in or out of species within a reaction cell. PHREEQCRM solves the equilibrium reactions and returns the updated species concentrations back to the flow model. (102)

2.4. Model Coupling Interface

A coupling interface is developed to connect the LBM flow solver PALABOS and the geochemical reaction solver PHREEQCRM in a dynamic fashion to solve pore-scale reactive transport problems. Both of the solvers are written in the C++ programming language, with APIs exposed to users to implement new functionalities, so the target of this task is to build an interface that is easy to use while maintaining high communication efficiency.

2.4.1. Data Pre-/Post-Processing and Script Input

The data preprocessing includes a series of functions for preparing the simulation data for fluid flow or chemical reaction domains. For instance, one of the functions is used to process the micro-CT scanned images to the simulation model readable ASCII/binary format data structure. Other functions are used for simulation parameters’ determination based on the model constraints. Also, a unified XML format-based script reading system is developed to preprocess the user inputs. The user input script includes various parameters for PALABOS and PHREEQCRM as well as other auxiliary features. Meanwhile, the code loads the “.pqi” PHREEQC input script with detailed chemical reaction definitions at the model initialization stage. The data post-processing includes various help functions to execute the main simulations, result analyses, and visualizations.

2.4.2. Model Initialization

The model initialization stage includes two steps: (1) simulation domain and reaction-type assignment; (2) initial equilibrium of the flow and reaction models. In the first step, the coupled model establishes the simulation domain from loaded porous media geometry information and creates cells mappings from the initial fluid–solid distribution. In addition to setting up the inlet/outlet flow and chemical flux boundary conditions, the “Bounce-Back” or “No-Flow” conditions are assigned within the “solid” grid blocks in LBM flow and diffusion lattices. Furthermore, the associated “solid” cells in PHREEQCRM are treated as “inactive” for geochemical reactions.

In the second step, the PHREEQCRM executes the initial sets of reaction calculations based on the user-defined equations in those reaction “active” cells, then obtains one or more scalar fields with the initial species concentration. It is worth noting that only the kinetic and surface complexation reactions on a type of cells called “solid interface cells” (see more details in Section 2.4.3, Cell Mapping Model) are enabled, where no flow exists but diffusion is enabled for mass transfer between the solid surface and outside aqueous phase cells.

On the other hand, the model initializes the fluid flow lattice from the initial and flow boundary conditions and obtains the initial velocity distribution. The initial velocity distribution fields were coupled with one or multiple concentration fields from PHREEQCRM, and the main simulation is ready to launch.

2.4.3. Cell Mapping Model

A unified cell mapping model was proposed to handle the fluid–fluid and fluid–solid interactions and geochemical reactions within different cell types. A dynamic mapping system is built on top of the porous media geometry, which initializes before the main simulation. The system includes four cell types:(1)

Fluid cells: cells represent the void space in the porous media

(2)

Fluid interface cells: cells belong to “fluid cells” and are positioned at the outer surface layer cells of the “fluid cells”

(3)

Solid cells: cells without any fluid flow and reaction enabled, defined from the initial input geometry

(4)

Solid interface cells: cells belong to and occupy the outer surface layer of the “solid cells”

The fluid cells host the LBM advection–diffusion and the PHREEQC equilibrium reactions. No flow or reactions happen within the solid cells. However, the equilibrium, kinetics, and surface complexation reactions on the outer layer (solid surface cells) are enabled. Meanwhile, the solid dissolution and precipitation processes are implemented within solid interface cells. The solid interface layer is still considered as “solid” due to the disabled fluid flow, but the diffusion is enabled for information exchange between them and aqueous cells since the solid might be partially dissolved. Another type of cells called fluid interface cells are also defined to handle future potential applications when complex fluid–fluid interactions are involved and in cases where the pore-scale geometry changes due to mineral precipitations.

Table 2 summarizes the advection–diffusion and reaction availability within the proposed four cell types.

Table 2. Compatibility of the Physics Processes and the Four Mapping Cellsa

fluid cellsfluid interface cellssolid cellssolid interface cellsfluid flowyesyesnonodiffusionyesyesnoyesequilibrium reactionyesyesnoyeskinetic reactionnononoyesdissolutionnononoyesprecipitationnoyesnoyessurface complexationnononoyes

The solid cells are inactive, whereas fluid cells only host advection–diffusion–reaction in the aqueous phase, and the interface cells are located at the boundary of the solid–fluid interfaces, hosting the fluid–solid interaction-associated processes except for the fluid flow.

2.4.4. PALABOS Data Processor

PALABOS has a flexible and powerful built-in data processor function, utilizing the C++ template programming features. The data processor allows the user to create flexible operations between lattices, scalar fields, and tensor fields in parallel. This function provides the foundation for coupling other models within PALABOS or from the third-party code. Various customized data processors are developed to couple the chemical concentration fields with the fluid flow lattice. In this work, a cross-lattice data processor is built to connect the fluid flow velocity field lattice and multiple advection–diffusion field lattices. Hence, in each simulation time step, the displacements caused by fluid flow are imposed on the concentration fields for all chemical species in the aqueous phase.

2.4.5. PHREEQCRM Basic Function Feature

PHREEQCRM shares a majority of the functionalities and features with the standalone version of PHREEQC, including a basic function input interpreter. Under specific formatting standards, the interpreter loads the user-defined input script (.pqi file). In addition to the pre-defined keywords followed by different reaction equations, basic mathematic operations and arithmetic functions can be defined in the script for a single reaction cell in PHREEQCRM. Furthermore, the defined operations are optimized in parallel similar to the reaction solver. Thus, customized analytical or numerical analysis on concentration fields can be done efficiently without spending a significant amount of data communication overhead due to data transfer. For instance, in this work, the calculation of the zeta potential (ζ) is achieved using a PHREEQCRM basic function executed after the general internal geochemical reactions are solved during each iteration. ζ is obtained by ref (103)

𝜁=𝜓𝑒−𝜅𝑑s

(17)

where ds is the distance between the slipping plane and the stern layer in the electrical double layer model, ψ stands for the surface potential, and κ represents the inverse of the Debye length which is obtained by

𝜅−1=𝜖r𝜖0𝑘b𝑇2𝑁A𝑒2𝐼⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯√

where NA is the Avogadro’s number; e is the charge of an electron; I represents the ionic strength; ϵr and ϵ0 are the temperature-dependent relative permittivity of the solution and the vacuum permittivity (8.8541878128–12 F m–1), respectively; kb = 1.38064852–23 J K–1 is the Boltzmann constant; and T is the temperature in Kelvin.(18)

2.5. AI-Assisted SCM Model Optimizations

The SCM describes solute adsorption from the aqueous phase to the mineral surface. (104) It is also considered one of the most explainable models for the fluid–solid interactions between the crude–oil–brine–rock (COBR) interfaces during advanced waterflooding EOR processes. One of the challenges in building the SCM models for COBR investigations is the adjustment of the SCM reaction constants (log K) based on the experimental measured zeta potentials (ζ). To validate the SCM model integration for the coupling code, an SCM model based on the experimental data by Tetteh et al. (105,106) is developed and optimized for the brine/calcite interface. Instead of estimating and calibrating the ζ within a bulk fluid solution, a more realistic 2D simulation domain is created with a centered spherical calcite grain surrounded by solutions with various salinities. Then, the PALABOS–PHREEQCRM coupled solver is used to obtain the surface properties such as ζ, surface potentials, and Debye length. The simulated ζ over the surface of the calcite grain is averaged and used as the target to compare with the experimentally measured ones and re-adjust the reaction constants if they do not match. Since the directions for adjusting the sets of reaction constants are unknown, and each one of the numerical simulations takes time, an optimization framework is developed to assist the parameter adjustment tasks. The framework connects the geochemically coupled PALABOS–PHREEQCRM simulator with a multi-layer perceptron (MLP) neural network (NN). A series of pre-simulated model parameters with the calculated ζ are prepared as the pre-train data for NN. The NN and simulator communicate with each other interactively to update the log K and match the target ζ as much as possible. The conventional mean squared error (MSE) loss function and the stochastic gradient descent (SGD) optimizer are used during the interactive training and optimization process.

Figure 4 depicts the overall workflow of the proposed optimization framework.

Figure 4 📷Figure 4. Illustration of the workflow for NN-assisted SCM model parameter optimization process.

3. Model Validation and Results

ARTICLE SECTIONS

Jump To

In this section, the model setup and validation results for the proposed framework are presented. The following model validation experiments are included: (1) an LBM single-phase fluid flow in a 2D channel used to verify basic fluid flow; (2) an LBM diffusion in a 2D channel experiment to verify basic Fick’s diffusion; (3) a reactive transport in a 2D channel experiment with a dissolving circular-shaped calcite grain (this case is used to verify the reactive transport coupling and surface reactions/calcite dissolution kinetics); and (4) a static SCM model in a 2D domain with artificial neural network (ANN) optimization (this case is used to verify the SCM modeling capabilities integrated in the framework) as well as the ANN optimization feature.

3.1. LBM Single-Phase Fluid Flow in a 2D Channel

The first validation experiment designed is a 2D Poiseuille flow scenario, as shown in Figure 5. A 2D rectangular-shaped simulation domain is created with a 300 × 100 lattice resolution, which represents a real-world physical dimension of 6 cm × 2 cm rectangle. An artificial pressure gradient (2 × 10–8 [Pa]) is applied across the channel with an inlet on the left and an outlet on the right. The no-flow boundary conditions are imposed on the top and bottom walls. Water is assumed as the fluid used in this experiment at a standard temperature of 25 °C, and it has a density of 997[kg/m3]

, viscosity of 8.891 × 10–4 [Pa·s], and kinematic viscosity of 8.917×10−7[m2/s]

. The fluid flow is assumed to be a steady-state laminar flow, and the viscous dissipation and gravity effects are neglected.

Figure 5 📷Figure 5. Model setup for the LBM single-phase fluid flow in a 2D channel.

To verify the solutions, the analytical solutions are derived based on the mass and momentum conservation equations

∇·V⇀=∂𝑢∂𝑥+∂𝑣∂𝑦=0

(19)

𝜌(∂V⇀∂𝑡+V⇀·∇V⇀)=𝜇∇2V⇀−∇𝑝+𝜌g⇀

(20)

After simplifying the momentum equation with the integration of boundary conditions, the below equation is obtained for the velocity profile u(y)

𝑢(𝑦)=−𝐻22𝜇d𝑝d𝑥[1−(𝑦𝐻)2]

(21)

where y is the coordinate across the width of the channel, H is the width of the channel, μ is the fluid viscosity, and dp/dx represents the assigned pressure gradient along the length of the channel. The average velocity is computed as

𝑢(avg)=−𝐻23𝜇d𝑝d𝑥

(22)

From eq 22, the maximum velocity calculated at the center is 1.87436 × 10–8 [m/s] and the average velocity is 1.24970 × 10–8 [m/s]. The same physical properties are applied and converted to lattice units, and the LBM simulation is conducted. The pressure gradient is represented by the lattice density gradient map shown in Figure 6.

Figure 6 📷Figure 6. Lattice density map generated from the given pressure gradient.

The 2D velocity profile obtained from LBM simulation after the steady state is reached, as shown in Figure 7.

Figure 7 📷Figure 7. Lattice velocity distribution induced by the pressure gradient of the fluid flow channel.

The 2D velocity distribution is evaluated from the LBM simulation, and the computed maximum (1.88450 × 10–8 [m/s]) and average 1.24371 × 10–8 [m/s] velocities after unit conversion agree with the analytical solutions with minor errors (0.48 and 0.54%, respectively). Furthermore, additional 10 locations across the width of the channel from the LBM simulated velocity profile are extracted, and they show good agreement with the analytical solutions, as shown in Figure 8.

Figure 8 📷Figure 8. LBM simulated velocities (red dots) against the analytical solutions (blue curve).

3.2. LBM Diffusion in a 2D Channel

A numerical experiment is designed to verify the capability of the LBM model in solving diffusion problems. The same simulation domain used in section 3.1 with the same resolution with some modifications (Figure 9) is re-used in this experiment. Instead of creating a fluid flow inlet from the left end, a single chemical species with a constant concentration c0 = 1.0 is imposed. The right-hand side of the channel is assumed to be infinitely long (cinf = 0), and the boundary conditions and boundary effects are not considered for top and bottom walls.

Figure 9 📷Figure 9. Model setup for the LBM diffusion validation in a 2D channel.

The following assumptions are made in our model setup: (1) there is no boundary effects such as surface adsorption to influence the diffusion; (2) the process is considered isothermal such that the rate of diffusion is not affected by the temperature; (3) the process is considered isotropic diffusion, which is the same in every direction and is characterized by a single diffusion coefficient D; and (4) the flow is one-dimensional Fick’s diffusion with a constant diffusion coefficient D = 10 × 10–6 [m/s]. The analytical solution of the concentration profiles is solved from Fick’s second law

∂𝑐∂𝑡=𝐷∂2𝑐∂𝑥2

(23)

The corresponding analytical solutions for the semi-infinite channel with a constant species concentration on one side is

𝑐(𝑥,𝑡)=𝑐oerfc(𝑥2𝐷𝑡⎯⎯⎯⎯√)

where co is the constant concentration (co = 1) on the left terminal, D represents the diffusion coefficient, t is the time, x is the location from the left terminal, erfc stands for the complementary error function, and c(x,t) describes the concentration of the distance x from the left terminal at time t. 7 time points (t = 1, 20, 50, 100, 150, 200, 250 s) are assigned as the testing time for the concentration profile comparison evaluations.(24)

A pure concentration field is created for single species in the LBM model without implementing any fluid flow in the channel. The lattice diffusivity is obtained from the physical diffusion coefficient and the domain lattice time. Figure 10 illustrates the concentration distribution at the initial stage (t = 1 s) from the LBM simulation. The 1D concentration profile is extracted from the center line across the domain as shown in Figure 10 at the given time steps.

Figure 10 📷Figure 10. Cross section (dashed line) at the center (y = (1/2)*H) of the 2D channel for concentration data extraction, where H is the width of the channel along the y-axis.

Figure 11 demonstrates the species concentration distribution at the given time points and the evaluation of the concentration profile at the center crossline at those time steps.

Figure 11 📷Figure 11. 2D concentration distribution from LBM simulation at t = 20, 50, 100, 150, 200, and 250 s.

The comparison results are shown in Figure 12. Good agreement is observed between the analytical solutions (curves) and LBM simulated results (dots).

Figure 12 📷Figure 12. Concentration profile comparison between LBM simulation and analytical solutions at t = 1, 20, 50, 100, 150, 200, and 250 s.

The errors are measured quantitively regarding three metrics: mean absolute error (MAE), MSE, and coefficient of determination (R2) which are defined in Appendix A. The evaluation results are shown in Table 3. The LBM simulated concentration profiles are in very good agreement with the analytical solutions at all time points except at t = 1 s. This can be attributed to the fluctuations at the beginning of the simulation, and the error becomes significantly smaller as the time steps go by.

Table 3. Errors between Analytical Solutions and LBM Simulation Results in Terms of MAE, MSE, and R2 Metrics

MAE errorMSE errorR2 scoreT = 1 s0.00300.00030.9772T = 20 s0.00062.1829 × 10–60.9999T = 50 s0.00045.4259 × 10–70.9999T = 100 s0.00031.9075 × 10–70.9999T = 150 s0.00031.1620 × 10–70.9999T = 200 s0.00044.2959 × 10–70.9999T = 250 s0.00093.1512 × 10–60.9999

3.3. Coupled Model for Advection–Diffusion–Reaction Transport in a 2D Channel with Calcite Dissolution Kinetics

The model validation results for single-phase fluid flow and pure diffusion from Sections 3.1 and 3.2 create a solid foundation for further reaction coupling. In this section, the advection–diffusion–reaction capabilities of the proposed framework are presented and verify the model performance by using the benchmark problem proposed by Molins et al. (107) This reactive transport benchmark problem is designed to evaluate the dynamic concentration fields in a 2D channel with a circular-shaped calcite grain located at the center. The average reaction rate is evaluated from the concentration differences between the inlet and the outlet. The simulation domain (Figure 13) is a rectangular-shaped 2D flow channel with a length of 0.1 cm and a width of 0.05 cm. A representing lattice grid is created with a resolution of 256 × 128. A circular-shaped calcite grain with a diameter of 0.02 cm is placed at the center of the channel. A constant velocity (0.12 cm/s) flow boundary condition is imposed at the inlet (left terminal). Meanwhile, an outflow condition is assigned to the outlet at the right terminal, indicating that no pressure gradient exists, or any driving force/resistance applied at the outlet. The top and bottom walls are non-reactive and impermeable and have no-flow boundaries.

Figure 13 📷Figure 13. 2D simulation configurations of this numerical experiment. The flow channel is filled with HCl solution, while it continuously flows through the channel with a constant flow rate from the inlet (left) to outlet (right). The calcite grain is placed at the center of the domain interacting with the solution. Upper and lower walls are impermeable and non-reactive boundaries. This figure is re-generated based on ref (107).

In the initial stage, the flow channel is assumed to be filled with hydrochloric acid (HCl) with a concentration of 10–2 [mol/L], which results in a pH = 2.0. After the injection starts, the solution with a concentration of HCl is pushed into the channel continuously. Meanwhile, the irreversible heterogeneous reaction starts to dissolve the calcite based on the stoichiometric equation

CaCO3(s)+H+→Ca2++HCO−3

(25)

and the kinetic reaction rate factor (rf) for the calcite dissolution is dependent on the H+ concentration

𝑟𝑓=𝑘H+𝛾H+𝑐H+

(26)

where kH+ is the rate constant with units of [mol·cm−2s−1], γH+ is the activity of H+ with units of [cm3·mol−1], and the cH+ is the concentration of H+ with units of [mol·cm−3] or [mol·L−1]

The detailed parameters of the benchmark problem are shown in Table 4.

Table 4. Simulation Parameters for the Reactive Transport Benchmark Problema

200cm–1rate constantkH+10–4.05mol cm–2 s–1activity coefficientγH+1000cm3 mol–1inlet concentration (pH = 2)C10–5mol cm–3Reynolds numberRe = uin ω ν–10.6 Péclet numberPe = uin ω D–1600parameterssymbolvalueunitsfluid densityρ1g cm–3kinematic viscosityν10–2cm2 s–1diffusion coefficientD10–5cm2 s–1inlet velocityuin0.12cm s–1width of the channelω0.05cmspecific grain reactive area📷

Simplified from ref (107).

The evolution of the concentration profile is evaluated according to the methods from Molins et al. (107) with the assumptions that several impact factors need to be considered: (1) the velocity distribution influenced by the inlet flow rate and the calcite grain; (2) the instantaneous equilibrium reactions in the fluid phase; (3) the kinetic reaction rate on the calcite surface; and (4) the multicomponent diffusion in the fluid phase as well as in the diffusive boundary layer defined on the calcite surface. In addition, the chemical reactions are assumed not to influence the fluid flow. The average reaction rate (R) is calculated as

𝑅=𝑄(𝑐out−𝑐in)𝜉𝐴

(27)

where ξ is the stoichiometric coefficient, A stands for the reacting surface area of the calcite, cin represents the universal inlet concentration from the boundary condition, and cout is evaluated as the flux-weighted-average concentration across the outlet boundary

𝑐out=∫𝑐𝒖·d𝑠/𝑄

(28)

where Q is obtained from the integration across the outlet

𝑄=∫𝒖·d𝑠

(29)

To calculate the reaction rate regarding the volumetric properties, a height of 1 cm is attributed to the 2D simulation domain for the benchmark problem.

Figure 14A demonstrates the steady-state velocity distribution from the LBM fluid flow lattice domain; the red color represents relatively higher velocity located at the top and bottom surface of the calcite grain.

Figure 14 📷Figure 14. (A) Velocity distribution under the steady-state condition; (B) 2D Ca2+ concentration distribution in the steady state; the higher Ca2+ concentration region surrounding the grain indicates the stable dissolution of CaCO3 on the grain surface.

Figure 14B shows the Ca2+ concentration distribution at the same moment with respect to Figure 14A. The formation of the Ca2+ concentration thin layer is observed around the grain as a result of the calcite dissolution, and the teardrop-shaped overall concentration pattern indicates the impacts of the flow velocity and direction.

Similarly, Figure 15A,B illustrates the H+ concentration evolution at the start of the simulation (t = 0.2 s) and after reaching a steady state at around t = 1.5 s. The H+ concentration distribution pattern demonstrates the H+ consumption during the calcite dissolution processes, which matches the Ca2+ concentration distribution pattern shown in Figure 14B.

Figure 15 📷Figure 15. (A) 2D H+ concentration distribution map; the color contour illustrates the beginning stage (t = 0.2 s) of the consumption of H+ due to calcite dissolution kinetic reaction on the solid surface. (B) 2D H+ concentration distribution map in the steady state (t ≥ 1.5 s).

The average effluent H+ concentration evolution at the right terminal of the channel directly from the concentration field from the LBM model is evaluated. Equation 27 is used to estimate the average reaction rate, and the results are compared according to Mollins et al. (107) The results are shown in Figure 16 for the outlet average effluent H+ concentration evolution from the beginning of the simulation until the steady state at 1.5 s. The overall trend from PALABOS–PHREEQCRM agrees with the other five codes despite the initial bump. The small bump is caused by the equilibrium reaction results from the initialization step before the fluid starts to flow during the PHREEQCRM initialization stage. Specifically, we believe that the PHREEQCRM generates slightly higher H+ concentrations (about 1.4%) at the assigned temperature of 40 °C; even we explicitly specified pH = 2.0. However, the effluent concentration reaches the desired equilibrium when the calcite dissolution kinetic reaction dominates the process, although the initial condition for the vortex method is different (pH = 7.0) compared to other methods, and the slightly different initial conditions do not affect the later and final solutions.

Figure 16 📷Figure 16. Averaged effluent H+ concentration evolution from this study and the other five codes. (107)

Meanwhile, Figure 17 shows the average dissolution rate over time for the PALABOS–PHREEQCRM framework and the other five codes. The dissolution rate curve and the trend of the proposed model demonstrate good agreement among the other five codes. In addition, PALABOS–PHREEQCRM has a slightly higher rate at the endpoint of this comparison with a rate of 4.59×10−8(mol·cm2·s−1)

; however, it is comparable to the other “Lattice–Boltzmann” approach 4.57×10−8(mol·cm2·s−1)

described in the benchmark problem (Table 5). This may be due to the slightly increased surface area attributed to the nature of the LBM-based methods that use “staircase”-shaped pixels/voxels to replace smooth boundary curves.

Figure 17 📷Figure 17. Averaged reaction/dissolution rate evolution from this study and the other five codes. (107)

Table 5. Average Steady-State Dissolution Rate from This Study Compared to Other Codes Described by Reference (107)

codesurface area (cm2)grain volume (cm3)average rate (mol cm–2 s–1)theoretical0.06283.14 × 10–4 this study0.06283.14 × 10–44.59 × 10–8OpenFOAM-DBS0.06283.14 × 10–44.18 × 10–8lattice–Boltzmann0.06283.14 × 10–44.57 × 10–8vortex0.06283.14 × 10–44.27 × 10–8

3.4. Static SCM Model and ANN Optimizations

In this section, the model validation for the SCM model and associated optimization using an ANN are presented. The SCM model is built based on the experimental work from Tetteh et al., (105) to investigate the electrokinetics on limestone surfaces when various brine compositions and salinities are applied. The rock sample is taken from Indiana limestone, which contains 99% calcite, (108) and the solid is assumed to be 100% calcite in the modeling process for simplicity. Tetteh et al. (109) synthesized the brines based on the water compositions from the Lansing Kansas City (LKC) group and then diluted the brine to create low salinity and seawater samples. In addition, they synthesized five single salt solutions to represent the low salinity brine with the same ionic strength. Table 6 shows the detailed brine composition.

Table 6. Brine Compositions Used in This Studya

brineCa2+Mg2+Na+K+Cl–SO42–IS (mol/L)TDS, ppm0.04-NaCl009200141800.0423380.04-KCl0001564141800.0429820.04-CaCl253400094500.0414800.04-MgCl203240094500.0412690.04-Na2SO4006130012810.041894FWS11,000280048,000500101,9132603.27164,473SWS2200560960010020,383520.6532,895LSW134345856124330.042006

The concentrations are in units of ppm prepared by Tetteh et al. (109) FWS: formation water; SWS: seawater; LSW: low salinity water.

The associated zeta potentials (ζ) are measured at the brine/calcite interfaces at both 25 °C (atmospheric conditions) and 40 °C (LKC reservoir conditions).

Extensive investigations have been done to study the calcite surface complexation reaction models and the reaction constants. (105,109−113) The proposed reactions and reaction constants are summarized in Table 7. To keep consistency with the experimental work by Tetteh et al., (105) the same calcite surface site density (4.95 sites/nm2) and calcite-specific surface area (1 m2/g) are deployed. Meanwhile, the PCO2 = 10–3.4 atm condition is used in the SCM model with calcite as equilibrated phases.

Table 7. Surface Complexation Reaction Constants for the Brine/Calcite Interface Proposed in the Literature

reactions (brine/calcite interface)Pokrovsky et al. (1999)Wolters et al. (2008)Hiorth et al. (2010)Brady et al. (2012)Brady et al. (2016)Tetteh et al. (2020)>CaOH + H+ ↔ >CaOH2+11.512.212.911.8511.8511.85>CaOH2+ + SO42– ↔ CaSO4– + H2O2.892.892.12.12.12.1>CaOH + HCO3– ↔ >CaCO3– + H2O5.64.93.325.84.285.8>CO3H ↔ >CO3– + H+–5.1–4.9–4.9–5.1–5.1–5.1>CO3H + Ca2+ ↔ >CO3Ca+ + H+–1.7–2.8–3.16–2.6–2.6–4.4>CO3H + Mg2+ ↔ >CO3Mg+ + H+–2.2–2.2–3.17–2.6–2.6–4.4

To validate the SCM model built on top of the PALABOS–PHREEQC framework, a 2D static model is designed instead of a simpler batch reactor shown in Figure 18. Brines with various compositions are filled in the closed system, and a pure calcite grain is placed at the center. To simplify the process, the equilibrium reactions are assumed to be instantaneous within the aqueous phase. The surface cells of the calcite grain (the “solid interface cells”) are set active to kinetic and calcite surface complexation reactions.

Figure 18 📷Figure 18. 2D PALABOS–PHREEQCRM SCM model validation scheme.

After the model initialization, PHREEQCRM calculates the SCM model-related properties at the grain surface cells, such as the surface charge, Debye length, and zeta potential (ζ). The average surface ζ is calculated and compared to the experimentally measured values. Furthermore, the SCM reaction constants defined in the PHREEQCRM input script file are adjusted to improve the matching. Instead of manual tuning, an ANN is developed to optimize the reaction constants automatically.

The ANN is built with a relatively simple MLP network for this work on top of the TensorFlow 2 framework. The network itself has three hidden layers with 50, 100, and 50 neurons, respectively. “ReLU” activation functions are assigned to each hidden layer, and a dropout layer with a dropout rate of 0.2 is attached after each hidden layer to alleviate potential overfitting issues. The weights of the layers are initialized with random normal distribution before training. The adaptive moment and learning rate estimation (Adam) optimizer and the standard MSE are used during the network training. The pre-trained data are shuffled within the data loader before being fed in the main training loop with a batch size of 1024. The network inputs contain 9 features including 8 simulated zeta potential values from the coupled numerical model plus the temperature. The outputs are 6 SCM reaction constants that we would like to optimize such that the resulted zeta potentials can match the experimental measurements. The input data are normalized using “standardization” scaling which scales the inputs separately by subtracting the mean and dividing by the standard deviation, such that the data distribution will have a mean of 0 and a standard deviation of 1 for better fitting performance.

At the beginning of the optimization, parallel simulations are performed on a wide array of parameter (six reaction constants) combinations and generate the corresponding ζSIM values for ANN pre-training. The data generated for pre-training containing about 20,000–40,000 samples “warm up” the ANN in the first few training epochs, which leads to a general desired direction to converge. Note that the ANN takes ζSIM as inputs and predicts reaction constants log K. The Ki is used to represent log K values for simplicity. Furthermore, the interactive optimization is performed by (1) obtaining a set of Ki from ANN predictions; (2) assigning the ANN predicted Ki to PALABOS–PHREEQCRM and running coupled numerical simulations for all 8 solution cases under given temperature conditions; the resulting 8 ζSIM values are used to compare the experimental measurements (ζEXP) to compute the loss function; (3) representing the differences between ζSIM and ζEXP by the MSE loss

MSE=1𝑛∑𝑖=1𝑛(𝜁EXP𝑖−𝜁SIM𝑖)2

(30)

where n is the number of ANN predictions, 8 is used since there are 8 Ki, 𝜁EXP𝑖 is the experimentally measured zeta potential and 𝜁SIM𝑖

represents the numerical model simulated zeta potentials; and (4) updating the ANN weights based on the MSE loss; (5) using the updated ANN to predict a new set of Ki and enter the next iteration. The iteration is terminated until it converges where the MSE loss stops decreasing and reaches a “steady-state” condition.

Table 8 demonstrates the Ki obtained from the interactive optimizations for both 25 and 40 °C conditions after convergence.

Table 8. Optimized SCM Reaction Constants for Brine/Calcite Interface for both 25 and 40 °C Conditions

reactions (brine/calcite interface)25 °C40 °CTetteh et al. (2020)>CaOH + H+ ↔ >CaOH2+11.08410.95711.85>CaOH2+ + SO42– ↔ CaSO4– + H2O1.9472.1072.1>CaOH + HCO3– ↔ >CaCO3– + H2O5.5675.4855.8>CO3H ↔ >CO3– + H+–4.578–4.527–5.1>CO3H + Ca2+ ↔ >CO3Ca+ + H+–3.614–3.454–4.4>CO3H + Mg2+ ↔ >CO3Mg+ + H+–3.352–3.157–4.4

Figure 19 shows the values of the MSE loss evolution during the optimizations for brine/calcite interface SCM models at 25 and 40 °C. A steady decline is observed initially and further converges after some oscillations.

Figure 19 📷Figure 19. MSE loss decay during the optimization for k values of the brine/calcite interface SCM under 25 °C (left) and 40 °C (right) conditions.

Figure 19 demonstrates the examples of the convergence processes of ζSIM for FWS and LSW calculated during the ANN optimizations. A similar convergence pattern is observed with the MSE loss decay curves, which indicates that the ANN optimization framework is capable of minimizing the difference between SCM predicted ζSCM and experimentally measured ζEXP, with the reaction constants Ki adjusted automatically (Figure 20).

Figure 20 📷Figure 20. Evolution of the optimized brine/calcite ζSCM for FWS (A) and LSW (B) under 25 °C conditions.

The optimized Ki are used to calculate the final ζSIM and compared to ζEXP as well as results from Tetteh et al. (105)Figures 21 and 22 summarize the comparison for five single salt solutions and FWS/SWS/LSW cases. The ANN optimized results show good agreement with the tuned SCM model by Tetteh et al. (105) with a similar pattern. The predicted ζSIM values are consistent with the experimental measurements for the five single salt solutions except for CaCl2 and MgCl2, where slightly reverse surface charges are obtained. This can be attributed to the increased adsorption of Ca2+ and increases >CO3Ca+ concentration due to calcite dissolution on the rock surface. (105) It is noticeable from the chart that the ζSIM from ANN optimized models are slightly closer to the ζEXP in most of the cases in Figure 21, indicating better Ki selection from ANN.

Figure 21 📷Figure 21. Comparison of ζEXP and ζSIM based on the reaction constants from ANN optimization and Tetteh et al. (105) regarding the five single salt solution cases.

Figure 22 📷Figure 22. Comparison of ζEXP and ζSIM based on the reaction constants from ANN optimization and Tetteh et al. (105) regarding the FWS/SWS/LSW cases.

Regarding the FWS, SWS, and LSW cases, all SCM models for FWS and SWS are found to be over-estimated ζ and under-estimated for LSW cases. Meanwhile, the ANN results show slightly more overestimation, while less under-estimated ζ in LSW cases. One possible explanation for the overestimation can be attributed to to the predicting accuracy issues of using the default PHREEQC thermodynamic database (phreeqc.dat) which is prone to not working well in high salinity scenarios. (102,114)

To quantitively evaluate the prediction performance, the MSE loss and coefficient of determination (R2) are calculated for three groups: five single salt solutions, FWS/SWS/LSW, and all cases; the two evaluation metrics are defined in Appendix A. The result of the comparison is shown in Table 9, and it shows that the ANN generates better overall prediction accuracy with smaller MSE loss. In the meantime, ANN predicts much less error in the single salt solution cases. However, it underperforms slightly in FWS/SWS/LSW cases compared to the SCM models from Tetteh et al. (105) It is worth noting that no additional weights are added to the constructed loss function for each individual case during the optimizations. Thus, the ANN tends to reduce the overall error of the fitting process instead of focusing on a single comparison.

Table 9. Performance Comparison of the SCM Predicted ζ Values

all casesfive single salt solutionFWS/SWS/LSWmetricsMSER2MSER2MSER2ANN71.4260.33464.8560.08449.4280.653Tetteh et al.80.8640.24695.580–0.35033.8040.763

The overall ANN-assisted interactive SCM model optimizations show promising performance without almost zero user intervention. Even though the results from ANN might not be the real global optimal (if one exists), it provides a user-friendly approach to automate the model tuning process, which saves a great amount of human work and time. It outperforms the manually tuned model in the experiments. Additionally, in cases where optimization is needed for more complex physical models, and there is no prior knowledge that exists to guide the manual tuning process, the proposed ANN-assisted optimization framework can offer significant advantages over the manual tuning method in an automatic fashion.

4. Conclusions, Discussion, and Future Work

ARTICLE SECTIONS

Jump To

In this paper, a numerical modeling framework is proposed that integrates LBM-based fluid flow and advection–diffusion models with geochemical reaction models to study pore-scale physiochemical processes. The framework and the demonstrated initial stage of numerical validations paved the way for modeling complex and realistic subsurface scientific and engineering problems such as advanced waterflooding and CO2 sequestration. Meanwhile, the framework is flexible and can be broadly applied in different pore-scale reactive transport simulations. In this work, PHREEQC is used as the geochemical reaction solver due to its capabilities of solving multiple reaction types as well as being able to couple third-party codes for multiphysics modeling. One of the variants of PHREEQC called PHREEQCRM is selected for the model building since it exposes low-level API for coupling purposes and works as a C++ library. Meanwhile, by using PHREEQCRM, much faster reaction solving and data communication between the two packages are obtained. In addition, a cell-mapping mechanism is developed within the framework to handle different reaction types at various cells, which allows the implementation of complex reactions such as SCM models on a solid surface. To validate the integrated modeling framework, a series of numerical experiments are designed to demonstrate the capabilities of this framework, ranging from simple 2D channel fluid flow and 2D diffusion to a 2D reactive transport benchmark problem. Furthermore, SCMs are added to the tests and validated against the SCM models with experimental measurements proposed in the literature. Meanwhile, a workflow using ANN to automatically assist the SCM model optimization is proposed. It successfully improves the SCM model tuning while illustrating the great potential to optimize more complex scenarios or models.

In Section 2 (Methodologies), the background and formulation of the LBM method in fluid flow and advection–diffusion problems are introduced at first, and the LBM-based PALABOS package is discussed briefly. Second, the geochemical reaction solver PHREEQC and its capabilities as well as the variant PHREEQCRM selected in the modeling framework are discussed. Furthermore, the development details of the coupling interface are demonstrated. Specifically, the model/data preprocessing and initialization that provide essential steps for the coupling are discussed. The cell-mapping mechanisms handle various reaction types within different mediums that allow modeling multiphysics in complex geometries. The PALABOS data processor and PHREEQCRM basic function features provide great capabilities and flexibility to transfer information between different physical or chemical fields while preserving parallel computing performance for built-in or customized functions. At last, the ANN-based automatic numerical model optimization workflow is illustrated and is applied to the SCM model tuning using the model framework.

In Section 3 (Model Validation and Results), the numerical validations are designed and performed for the model in various aspects. First, in the 2D LBM fluid flow-only validation experiment based on a channel Poiseuille flow problem, the proposed model shows good agreement with analytical solutions with a small error (around 0.5%). Furthermore, the code is verified for the 1D diffusion problem in a 2D channel, and the time-dependent concentration profiles are evaluated with analytical solutions. The results suggest that the proposed model has high accuracy (R2 ≥ 0.9997), although minimal performance loss is observed at the beginning of the simulation due to the initial oscillations. After validating flow or diffusion-only scenarios, the experiment was inspired and designed by the reactive transport benchmark problem proposed by Molins et al. (107) and used the same initial conditions for the simulation. The time-dependent averaged effluent H+ concentration is evaluated, and the averaged reaction/calcite dissolution rate as the metrics is estimated for comparison. The results demonstrate good agreement with the other five codes described in the benchmark problem. The last mode validation focuses on the SCM model verification for complex physical–chemical processes using an advanced waterflooding-associated mechanistic study as an example. In addition to the SCM model integration, an ANN-assisted automatic SCM model tuning workflow is proposed that dynamically interacts with the numerical model to optimize the reaction constants. The results show good zeta potential estimations from SCM reactions compared to experimental data. Meanwhile, the ANN optimized SCM model parameters improve the overall matching performance while significantly reducing human labor in the manual model adjustment.

It is worth noting that the proposed coupling framework has its limitations at this initial stage of study. For example, we made assumption in our initial modeling practice that the relatively high concentrations of chemical species (Ca2+ in this case) do not reversely affect the fluid transport process due to intermolecular interactions, such that the LBM fluid flow is purely driven by pressure differences and affected by the no-flow boundaries at the walls. We plan to include the intermolecular interaction effects in the future modeling approaches.

Another limitation of this modeling framework is that we have not fully explored and utilized the lower-level source code and functionalities of PHREEQCRM; thus, we do not have full control of some implementation details to meet our desired needs, for example, the initialization behavior with a small concentration bump described in Section 3.3. Future work would include the coupling improvements with in-depth integration of PHREEQCRM for customized modeling needs and adjustments.

Meanwhile, all the numerical experiments are designed and proceed in a 2D domain, and we target pore-scale ranging from sub-micrometers/micrometers to centimeters (rock core or a portion of rock core). The numerical validation in a 2D symmetrical channel indicates the applicability of a 3D symmetrical domain such as a cylindrical pipe. We reserve the 3D modeling scenarios for future studies. It is worth noting that we intended to validate the coupled model in a relatively simple pore geometry instead of a realistic pore network due to the limitation of the computational resources. We treated this work as a “proof of concept” type of study, which paved the way for further investigations on more complex and realistic pore structures by incorporating a computing cluster with abundant computational resources.

Author Information

ARTICLE SECTIONS

Jump To

Corresponding Author Reza Barati - Department of Chemical & Petroleum Engineering, University of Kansas, Lawrence, Kansas 66045, United States; 📷https://orcid.org/0000-0002-1064-9562; Email: [email protected]
Authors Siyan Liu - Department of Chemical & Petroleum Engineering, University of Kansas, Lawrence, Kansas 66045, United States; Computational Sciences and Engineering Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37830, United States; 📷https://orcid.org/0000-0003-2017-3251 Chi Zhang - Department of Meteorology and Geophysics, Institute of Meteorology and Geophysics, University of Vienna, Universität Wien, UZA II, Josef-Holaubek-Platz 2, Wien 1090, Austria Mohammad Kazemi - Department of Physics and Engineering, Slippery Rock University, Slippery Rock, Pennsylvania 16057, United States
NotesThe authors declare no competing financial interest.

Acknowledgments

ARTICLE SECTIONS

Jump To

This work was financially supported by the Department of Chemical and Petroleum Engineering, University of Kansas, and was partially supported by the Kansas Interdisciplinary Carbonates Consortium (KICC) as well as the ExxonMobil/GSA Student Geoscience Grant.

Appendix

ARTICLE SECTIONS

Jump To

Evaluation Metrics

In this study, we use MAE, MSE, and coefficient of determination (R2) to evaluate the model validation results and the ANN prediction performance. MAE is defined as

MAE=1𝑛∑𝑖=1𝑛|𝑌obs𝑖−𝑌pred𝑖|

(31)

MSE=1𝑛∑𝑖=1𝑛(𝑌obs𝑖−𝑌pred𝑖)2

(32)

𝑅2=1−∑𝑛𝑖=1(𝑌obs𝑖−𝑌¯pred)2∑𝑛𝑖=1(𝑌obs𝑖−𝑌¯obs)2

(33)

where Yiobs stands for the ith observation, 𝑌¯𝑜𝑏𝑠

is the mean value of the observations, and Yipred represents model or ANN predictions.

Nomenclature

ARTICLE SECTIONS

Jump To

CFD

computational fluid mechanics

LBM

lattice Boltzmann method

EOR

enhanced oil recovery

LBM collision operator

LBM relaxation time

lattice distribution function

fkEQ

lattice equilibrium distribution function

ρLB

lattice fluid density

PLB

lattice pressure

lattice speed of sound

uLB

lattice velocity

νLB

lattice kinematic viscosity

ωi

weight coefficient for i direction

SCM

surface complexation model

COBR

crude–oil–brine–rock

zeta potential

width of the flow channel

u(y)

flow velocity profile along y-axis

pressure

fluid viscosity

number of lattice nodes along the x-axis

number of lattice nodes along the y-axis

species concentration

constant species concentration at boundary

cinf

species concentration at the end of the assumed infinitely long channel

diffusion coefficient

reaction rate factor

rate constant

activity coefficient

Reynolds number

Péclet number

stoichiometric coefficient

reacting surface area

average reaction rate

neural network

ANN

artificial neural network

MLP

multi-layer perceptron

SGD

stochastic gradient descent

MAE

mean absolute error

MSE

mean squared error

coefficient of determination

Yiobs

observation values

Yipred

prediction values

reaction constant log K

FWS/SWS/LSW

formation water/seawater/low salinity water

References

ARTICLE SECTIONS

Jump To

This article references 114 other publications.

1Green, D. W.; Willhite, G. P.; Enhanced Oil Recovery; Henry L. Doherty Memorial Fund of AIME, Society of Petroleum Engineers, 1998; Vol. 6.Google Scholar

2Gluyas, J.; Mathias, S. Geological Storage of Carbon Dioxide (CO2): Geoscience, Technologies, Environmental Aspects and Legal Frameworks; Elsevier, 2013.Google Scholar

3Sheng, J. J. Critical Review of Low-Salinity Waterflooding. J. Pet. Sci. Eng. 2014, 120, 216– 224, DOI: 10.1016/j.petrol.2014.05.026[Crossref], [CAS], Google Scholar

4Lager, A.; Webb, K. J.; Black, C. J. J.; Singleton, M.; Sorbie, K. S. Low Salinity Oil Recovery - An Experimental Investigation. Petrophysics 2008, 49, 28– 35Google Scholar

5Willhite, G. P. Waterflooding; Society of Petroleum Engineers: Richardson, TX, 1986.[Crossref], Google Scholar

6Kazemi Nia Korrani, A.; Jerauld, G. R.; Sepehrnoori, K. Coupled Geochemical-Based Modeling of Low Salinity Waterflooding. SPE Improved Oil Recovery Symposium , 2014, No. 2008; Vol. 1–23.[Crossref], Google Scholar

7Shalabi, E. W. A. Modeling the Effect of Injecting Low Salinity Water on Oil Recovery from Carbonate Reservoirs, Dissertation, The University of Texas at Austin, 2014.Google Scholar

8Sanaei, A. Compositional Reactive-Transport Modeling of Engineered Waterflooding, Dissertation, University of Texas at Austin, 2019.Google Scholar

9Sanaei, A.; Tavassoli, S.; Sepehrnoori, K. Investigation of Modified Water Chemistry for Improved Oil Recovery: Application of DLVO Theory and Surface Complexation Model. Colloids Surf., A 2019, 574, 131– 145, DOI: 10.1016/j.colsurfa.2019.04.075[Crossref], [CAS], Google Scholar

10Chen, S.; Doolen, G. D. Lattice Boltzmann Method for Fluid Flows. Annu. Rev. Fluid. Mech. 1998, 30, 329– 364, DOI: 10.1146/annurev.fluid.30.1.329[Crossref], Google Scholar

11Shan, X.; Chen, H. Lattice Boltzmann Model for Simulating Flows with Multiple Pahses and Components. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1993, 47, 1815, DOI: 10.1103/physreve.47.1815[Crossref], [PubMed], [CAS], Google Scholar

12Huang, H.; Sukop, M. C.; Lu, X.-Y. Multiphase Lattice Boltzmann Methods: Theory and Application, 1st ed.; WILEY Blackwell, 2015.[Crossref], Google Scholar

13Yoon, H.; Kang, Q.; Valocchi, A. J. 12. Lattice Boltzmann-Based Approaches for Pore-Scale Reactive Transport. Rev. Mineral. Geochem. 2015, 80, 393– 432, DOI: 10.1515/9781501502071-012[Crossref], Google Scholar

14Sethian, J. A. A. Fast Marching Level Set Method for Monotonically Advancing Fronts. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 1591, DOI: 10.1073/pnas.93.4.1591[Crossref], [PubMed], [CAS], Google Scholar

15Tryggvason, G.; Bunner, B.; Esmaeeli, A.; Juric, D.; Al-Rawahi, N.; Tauber, W.; Han, J.; Nas, S.; Jan, Y. J. A Front-Tracking Method for the Computations of Multiphase Flow. J. Comput. Phys. 2001, 169, 708– 759, DOI: 10.1006/jcph.2001.6726[Crossref], [CAS], Google Scholar

16Hirt, C. W.; Nichols, B. D. Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries. J. Comput. Phys. 1981, 39, 201, DOI: 10.1016/0021-9991(81)90145-5[Crossref], Google Scholar

17Shan, X.; Chen, H. Simulation of Nonideal Gases and Liquid-Gas Phase Transitions by the Lattice Boltzmann Equation. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1994, 49, 2941– 2948, DOI: 10.1103/physreve.49.2941[Crossref], [PubMed], [CAS], Google Scholar

18Gunstensen, A. K.; Rothman, D. H.; Zaleski, S.; Zanetti, G. Lattice Boltzmann Model of Immiscible Fluids. Phys. Rev. A: At., Mol., Opt. Phys. 1991, 43, 4320– 4327, DOI: 10.1103/physreva.43.4320[Crossref], [PubMed], [CAS], Google Scholar

19Swift, M. R.; Osborn, W. R.; Yeomans, J. M. Lattice Boltzmann Simulation of Nonideal Fluids. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1995, 75, 830, DOI: 10.1103/physrevlett.75.830[Crossref], Google Scholar

20Swift, M. R.; Orlandini, E.; Osborn, W. R.; Yeomans, J. M. Lattice Boltzmann Simulations of Liquid-Gas and Binary Fluid Systems. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1996, 54, 5041– 5052, DOI: 10.1103/physreve.54.5041[Crossref], [PubMed], [CAS], Google Scholar

21He, X.; Chen, S.; Zhang, R. A Lattice Boltzmann Scheme for Incompressible Multiphase Flow and Its Application in Simulation of Rayleigh-Taylor Instability. J. Comput. Phys. 1999, 152, 642– 663, DOI: 10.1006/jcph.1999.6257[Crossref], [CAS], Google Scholar

22Parmigiani, A. Lattice Boltzmann Calculations of Reactive Multiphase Flows in Porous Media, Thesis, University of Geneva, 2011; Vol. 129.Google Scholar

23Kingdon, R. D.; Schofield, P. A Reaction-Flow Lattice Boltzmann Model. J. Phys. A: Math. Gen. 1992, 25, L907– L910, DOI: 10.1088/0305-4470/25/14/008[Crossref], Google Scholar

24Kang, Q.; Zhang, D.; Chen, S. Displacement of a Three-Dimensional Immiscible Droplet in a Duct. J. Fluid Mech. 2005, 545, 41– 66, DOI: 10.1017/s0022112005006956[Crossref], Google Scholar

25Kang, Q.; Lichtner, P. C.; Zhang, D. Lattice Boltzmann Pore-Scale Model for Multicomponent Reactive Transport in Porous Media. J. Geophys. Res.: Solid Earth 2006, 111, 1– 12, DOI: 10.1029/2005jb003951[Crossref], Google Scholar

26Liu, M.; Mostaghimi, P. Pore-Scale Simulation of Dissolution-Induced Variations in Rock Mechanical Properties. Int. J. Heat Mass Transfer 2017, 111, 842– 851, DOI: 10.1016/j.ijheatmasstransfer.2017.04.049[Crossref], [CAS], Google Scholar

27Liu, M.; Mostaghimi, P. High-Resolution Pore-Scale Simulation of Dissolution in Porous Media. Chem. Eng. Sci. 2017, 161, 360– 369, DOI: 10.1016/j.ces.2016.12.064[Crossref], [CAS], Google Scholar

28Tian, Z.; Wang, J. Lattice Boltzmann Simulation of CO2 Reactive Transport in Network Fractured Media. Water Resour. Res. 2017, 53, 7366– 7381, DOI: 10.1002/2017wr021063[Crossref], [CAS], Google Scholar

29Di Palma, P. R.; Huber, C.; Viotti, P. A New Lattice Boltzmann Model for Interface Reactions between Immiscible Fluids. Adv. Water Resour. 2015, 82, 139– 149, DOI: 10.1016/j.advwatres.2015.05.001[Crossref], [CAS], Google Scholar

30Shan, X. Simulation of Rayleigh-Bénard convection using a lattice Boltzmann method. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1997, 55, 2780– 2788, DOI: 10.1103/physreve.55.2780[Crossref], [CAS], Google Scholar

31Huber, C.; Parmigiani, A.; Chopard, B.; Manga, M.; Bachmann, O. Lattice Boltzmann Model for Melting with Natural Convection. Int. J. Heat Fluid Flow 2008, 29, 1469– 1480, DOI: 10.1016/j.ijheatfluidflow.2008.05.002[Crossref], [CAS], Google Scholar

32Yoon, H.; Kang, Q.; Valocchi, A. J. 12. Lattice Boltzmann-Based Approaches for Pore-Scale Reactive Transport. Rev. Mineral. Geochem. 2015, 80, 393– 432, DOI: 10.1515/9781501502071-012[Crossref], Google Scholar

33Luan, H. B.; Xu, H.; Chen, L.; Sun, D. L.; Tao, W. Q. Numerical Illustrations of the Coupling Between the Lattice Boltzmann Method and Finite-Type Macro-Numerical Methods. Numer. Heat Transfer, Part B 2010, 57, 147– 171, DOI: 10.1080/15421400903579929[Crossref], [CAS], Google Scholar

34Yu, D.; Mei, R.; Shyy, W. A Multi-Block Lattice Boltzmann Method for Viscous Fluid Flows. Int. J. Numer. Methods Fluids 2002, 39, 99– 120, DOI: 10.1002/fld.280[Crossref], Google Scholar

35Luan, H.-B.; Xu, H.; Chen, L.; Sun, D.-L.; He, Y.-L.; Tao, W.-Q. Evaluation of the Coupling Scheme of FVM and LBM for Fluid Flows around Complex Geometries. Int. J. Heat Mass Transfer 2011, 54, 1975– 1985, DOI: 10.1016/j.ijheatmasstransfer.2011.01.004[Crossref], Google Scholar

36Chen, L.; He, Y.-L.; Kang, Q.; Tao, W.-Q. Coupled Numerical Approach Combining Finite Volume and Lattice Boltzmann Methods for Multi-Scale Multi-Physicochemical Processes. J. Comput. Phys. 2013, 255, 83– 105, DOI: 10.1016/j.jcp.2013.07.034[Crossref], [CAS], Google Scholar

37Sullivan, S. P.; Sani, F. M.; Johns, M. L.; Gladden, L. F. Simulation of Packed Bed Reactors Using Lattice Boltzmann Methods. Chem. Eng. Sci. 2005, 60, 3405– 3418, DOI: 10.1016/j.ces.2005.01.038[Crossref], [CAS], Google Scholar

38Chen, L.; Kang, Q.; Robinson, B. A.; He, Y.-L.; Tao, W.-Q. Pore-Scale Modeling of Multiphase Reactive Transport with Phase Transitions and Dissolution-Precipitation Processes in Closed Systems. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2013, 87, 43306, DOI: 10.1103/physreve.87.043306[Crossref], [CAS], Google Scholar

39Kang, Q.; Zhang, D.; Lichtner, P. C.; Tsimpanogiannis, I. N. Lattice Boltzmann Model for Crystal Growth from Supersaturated Solution. Geophys. Res. Lett. 2004, 31, GL021107, DOI: 10.1029/2004gl021107[Crossref], Google Scholar

40Chen, L.; Kang, Q.; Carey, B.; Tao, W. Q. Pore-Scale Study of Diffusion-Reaction Processes Involving Dissolution and Precipitation Using the Lattice Boltzmann Method. Int. J. Heat Mass Transfer 2014, 75, 483– 496, DOI: 10.1016/j.ijheatmasstransfer.2014.03.074[Crossref], Google Scholar

41Kang, Q.; Chen, L.; Valocchi, A. J.; Viswanathan, H. S. Pore-Scale Study of Dissolution-Induced Changes in Permeability and Porosity of Porous Media. J. Hydrol. 2014, 517, 1049– 1055, DOI: 10.1016/j.jhydrol.2014.06.045[Crossref], Google Scholar

42Chen, L.; Kang, Q.; Viswanathan, H. S.; Tao, W.-Q. Pore-scale study of dissolution-induced changes in hydrologic properties of rocks with binary minerals. Water Resour. Res. 2014, 50, 9343– 9365, DOI: 10.1002/2014wr015646[Crossref], Google Scholar

43Parkhurst, B. D. L. User’s Guide to PHREEQC ─ a Computer Program for Inverse Geochemical Calculations; U.S. Geological Survey, 1995.Google Scholar

44Charlton, S. R.; Parkhurst, D. L. Modules Based on the Geochemical Model PHREEQC for Use in Scripting and Programming Languages. Comput. Geosci. 2011, 37, 1653– 1663, DOI: 10.1016/j.cageo.2011.02.005[Crossref], [CAS], Google Scholar

45Parkhurst, D. L.; Wissmeier, L. PhreeqcRM: A Reaction Module for Transport Simulators Based on the Geochemical Model PHREEQC. Adv. Water Resour. 2015, 83, 176– 189, DOI: 10.1016/j.advwatres.2015.06.001[Crossref], [CAS], Google Scholar

46Liu, S.; Zhang, C.; Ghahfarokhi, R. B. A Review of Lattice-Boltzmann Models Coupled with Geochemical Modeling Applied for Simulation of Advanced Waterflooding and Enhanced Oil Recovery Processes. Energy Fuels 2021, 35, 13535– 13549, DOI: 10.1021/acs.energyfuels.1c01347[ACS Full Text 📷], [CAS], Google Scholar

47Kazemi Nia Korrani, A.; Jerauld, G. R.; Sepehrnoori, K. Coupled Geochemical-Based Modeling of Low Salinity Waterflooding. SPE Improved Oil Recovery Symposium , 2014, No. 2008; Vol. 1–23.[Crossref], Google Scholar

48Patel, R.; Perko, J.; Jacques, D.; De Schutter, G.; Ye, G.; Van Breugel, K. Lattice Boltzmann Based Multicomponent Reactive Transport Model Coupled with Geochemical Solver for Scale Simulations. Computational Methods for Coupled Problems in Science and Engineering , 2013; pp 806– 817.Google Scholar

49Patel, R. A.; Perko, J.; Jacques, D.; De Schutter, G.; Van Breugel, K.; Ye, G. A Versatile Pore-Scale Multicomponent Reactive Transport Approach Based on Lattice Boltzmann Method: Application to Portlandite Dissolution. Phys. Chem. Earth 2014, 70–71, 127– 137, DOI: 10.1016/j.pce.2014.03.001[Crossref], Google Scholar

50Patel, R. A.; Perko, J.; Jacques, D.; De Schutter, G.; Ye, G.; Van Breugel, K. A Three-Dimensional Lattice Boltzmann Method Based Reactive Transport Model to Simulate Changes in Cement Paste Microstructure Due to Calcium Leaching. Constr. Build. Mater. 2018, 166, 158– 170, DOI: 10.1016/j.conbuildmat.2018.01.114[Crossref], [CAS], Google Scholar

51Fazeli, H.; Patel, R.; Hellevang, H. Effect of Pore-Scale Mineral Spatial Heterogeneity on Chemically Induced Alterations of Fractured Rock: A Lattice Boltzmann Study. Geofluids 2018, 2018, 1– 28, DOI: 10.1155/2018/6046182[Crossref], Google Scholar

52Fazeli, H.; Patel, R. A.; Ellis, B. R.; Hellevang, H. Three-Dimensional Pore-Scale Modeling of Fracture Evolution in Heterogeneous Carbonate Caprock Subjected to CO2-Enriched Brine. Environ. Sci. Technol. 2019, 53, 4630– 4639, DOI: 10.1021/acs.est.8b05653[ACS Full Text 📷], [CAS], Google Scholar

53Fazeli, H.; Masoudi, M.; Patel, R. A.; Aagaard, P.; Hellevang, H. Pore-Scale Modeling of Nucleation and Growth in Porous Media. ACS Earth Space Chem. 2020, 4, 249– 260, DOI: 10.1021/acsearthspacechem.9b00290[ACS Full Text 📷], [CAS], Google Scholar

54Kazemi Nia Korrani, A.; Sepehrnoori, K.; Delshad, M. Coupling IPhreeqc with UTCHEM to Model Reactive Flow and Transport. Comput. Geosci. 2015, 82, 152– 169, DOI: 10.1016/j.cageo.2015.06.004[Crossref], [CAS], Google Scholar

55Nardi, A.; Idiart, A.; Trinchero, P.; de Vries, L. M.; Molinero, J. Interface COMSOL-PHREEQC (ICP), an Efficient Numerical Framework for the Solution of Coupled Multiphysics and Geochemistry. Comput. Geosci. 2014, 69, 10– 21, DOI: 10.1016/j.cageo.2014.04.011[Crossref], [CAS], Google Scholar

56Muniruzzaman, M.; Rolle, M. Modeling Multicomponent Ionic Transport in Groundwater with IPhreeqc Coupling: Electrostatic Interactions and Geochemical Reactions in Homogeneous and Heterogeneous Domains. Adv. Water Resour. 2016, 98, 1– 15, DOI: 10.1016/j.advwatres.2016.10.013[Crossref], [CAS], Google Scholar

57Parkhurst, D. L.; Kipp, K. L.; Engesgaard, P.; Charlton, S. R. PHAST─A Program for Simulating Ground-Water Flow, Solute Transport, and Multicomponent Geochemical Reactions, U.S. Geological Survey Techniques and Methods 6-A8; U.S. Geological Survey, 2004.[Crossref], Google Scholar

58Parkhurst, D.; Kipp, K.; Charlton, S. PHAST Version 2 - A Program for Simulating Groundwater Flow, Solute Transport, and Multicomponent Geochemical Reactions, Modeling Techniques, Book 6; U.S. Geological Survey, 2010.[Crossref], Google Scholar

59Diersch, H. J. G. FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media; Springer Science & Business Media, 2014.[Crossref], Google Scholar

60Muniruzzaman, M.; Rolle, M. Multicomponent Ionic Transport Modeling in Physically and Electrostatically Heterogeneous Porous Media With PhreeqcRM Coupling for Geochemical Reactions. Water Resour. Res. 2019, 55, 11121– 11143, DOI: 10.1029/2019wr026373[Crossref], [CAS], Google Scholar

61Fokina, D.; Muravleva, E.; Ovchinnikov, G.; Oseledets, I. Microstructure Synthesis Using Style-Based Generative Adversarial Networks. Phys. Rev. E 2020, 101, 43308, DOI: 10.1103/physreve.101.043308[Crossref], [CAS], Google Scholar

62Mosser, L.; Dubrule, O.; Blunt, M. J. Reconstruction of Three-Dimensional Porous Media Using Generative Adversarial Neural Networks. Phys. Rev. E 2017, 96, 43309, DOI: 10.1103/physreve.96.043309[Crossref], [PubMed], [CAS], Google Scholar

63Mosser, L.; Dubrule, O.; Blunt, M. J. Stochastic Reconstruction of an Oolitic Limestone by Generative Adversarial Networks. Transp. Porous Media 2018, 125, 81– 103, DOI: 10.1007/s11242-018-1039-9[Crossref], [CAS], Google Scholar

64Feng, J.; He, X.; Teng, Q.; Ren, C.; Chen, H.; Li, Y. Reconstruction of Porous Media from Extremely Limited Information Using Conditional Generative Adversarial Networks. Phys. Rev. E 2019, 100, 33308, DOI: 10.1103/physreve.100.033308[Crossref], [CAS], Google Scholar

65Liu, S.; Zhong, Z.; Takbiri-Borujeni, A.; Kazemi, M.; Fu, Q.; Yang, Y. A Case Study on Homogeneous and Heterogeneous Reservoir Porous Media Reconstruction by Using Generative Adversarial Networks. Energy Procedia 2019, 158, 6164, DOI: 10.1016/j.egypro.2019.01.493[Crossref], Google Scholar

66Shams, R.; Masihi, M.; Boozarjomehry, R. B.; Blunt, M. J. Coupled Generative Adversarial and Auto-Encoder Neural Networks to Reconstruct Three-Dimensional Multi-Scale Porous Media. J. Pet. Sci. Eng. 2020, 186, 106794, DOI: 10.1016/j.petrol.2019.106794[Crossref], [CAS], Google Scholar

67Varfolomeev, I.; Yakimchuk, I.; Safonov, I. An Application of Deep Neural Networks for Segmentation of Microtomographic Images of Rock Samples. Computers 2019, 8, 72, DOI: 10.3390/computers8040072[Crossref], Google Scholar

68Niu, Y.; Mostaghimi, P.; Shabaninejad, M.; Swietojanski, P.; Armstrong, R. T. Digital Rock Segmentation for Petrophysical Analysis With Reduced User Bias Using Convolutional Neural Networks. Water Resour. Res. 2020, 56, e2019WR026597 DOI: 10.1029/2019wr026597[Crossref], Google Scholar

69Karimpouli, S.; Tahmasebi, P. Image-Based Velocity Estimation of Rock Using Convolutional Neural Networks. Neural Networks 2019, 111, 89– 97, DOI: 10.1016/j.neunet.2018.12.006[Crossref], [PubMed], [CAS], Google Scholar

70Wang, Y. D.; Shabaninejad, M.; Armstrong, R. T.; Mostaghimi, P. Deep Neural Networks for Improving Physical Accuracy of 2D and 3D Multi-Mineral Segmentation of Rock Micro-CT Images. Appl. Soft Comput. 2021, 104, 107185, DOI: 10.1016/j.asoc.2021.107185[Crossref], Google Scholar

71Kamrava, S.; Tahmasebi, P.; Sahimi, M. Enhancing Images of Shale Formations by a Hybrid Stochastic and Deep Learning Algorithm. Neural Networks 2019, 118, 310– 320, DOI: 10.1016/j.neunet.2019.07.009[Crossref], [PubMed], [CAS], Google Scholar

72Kamrava, S.; Tahmasebi, P.; Sahimi, M. Linking Morphology of Porous Media to Their Macroscopic Permeability by Deep Learning. Transp. Porous Media 2020, 131, 427– 448, DOI: 10.1007/s11242-019-01352-5[Crossref], Google Scholar

73Tembely, M.; AlSumaiti, A. Deep Learning for a Fast and Accurate Prediction of Complex Carbonate Rock Permeability From 3D Micro-CT Images; Abu Dhabi International Petroleum Exhibition and Conference, 2019.[Crossref], Google Scholar

74Wu, J.; Yin, X.; Xiao, H. Seeing Permeability from Images: Fast Prediction with Convolutional Neural Networks. Sci. Bull. 2018, 63, 1215– 1222, DOI: 10.1016/j.scib.2018.08.006[Crossref], [CAS], Google Scholar

75Liu, S.; Zolfaghari, A.; Sattarin, S.; Dahaghi, A. K.; Negahban, S. Application of Neural Networks in Multiphase Flow through Porous Media: Predicting Capillary Pressure and Relative Permeability Curves. J. Pet. Sci. Eng. 2019, 180, 445– 455, DOI: 10.1016/j.petrol.2019.05.041[Crossref], [CAS], Google Scholar

76Liu, S.; Barati, R.; Zhang, C. Fast Estimation of Permeability in Sandstones by 3D Convolutional Neural Networks. SEG International Exposition and Annual Meeting , September 15, 2019; p D033S046R002.[Crossref], Google Scholar

77Wu, H.; Fang, W.-Z.; Kang, Q.; Tao, W.-Q.; Qiao, R. Predicting Effective Diffusivity of Porous Media from Images by Deep Learning. Sci. Rep. 2019, 9, 20387, DOI: 10.1038/s41598-019-56309-x[Crossref], [PubMed], [CAS], Google Scholar

78Wei, H.; Zhao, S.; Rong, Q.; Bao, H. Predicting the Effective Thermal Conductivities of Composite Materials and Porous Media by Machine Learning Methods. Int. J. Heat Mass Transfer 2018, 127, 908– 916, DOI: 10.1016/j.ijheatmasstransfer.2018.08.082[Crossref], Google Scholar

79Zhang, Z.; Hong, Y.; Hou, B.; Zhang, Z.; Negahban, M.; Zhang, J. Accelerated Discoveries of Mechanical Properties of Graphene Using Machine Learning and High-Throughput Computation. Carbon 2019, 148, 115– 123, DOI: 10.1016/j.carbon.2019.03.046[Crossref], [CAS], Google Scholar

80Santos, J. E.; Xu, D.; Jo, H.; Landry, C. J.; Prodanović, M.; Pyrcz, M. J. PoreFlow-Net: A 3D Convolutional Neural Network to Predict Fluid Flow through Porous Media. Adv. Water Resour. 2020, 138, 103539, DOI: 10.1016/j.advwatres.2020.103539[Crossref], Google Scholar

81Liu, M.; Kwon, B.; Kang, P. K. Machine Learning to Predict Effective Reaction Rates in 3D Porous Media from Pore Structural Features. Sci. Rep. 2022, 12, 5486, DOI: 10.1038/s41598-022-09495-0[Crossref], [PubMed], [CAS], Google Scholar

82Tahmasebi, P.; Kamrava, S.; Bai, T.; Sahimi, M. Machine Learning in Geo- and Environmental Sciences: From Small to Large Scale. Adv. Water Resour. 2020, 142, 103619, DOI: 10.1016/j.advwatres.2020.103619[Crossref], Google Scholar

83Wang, Y. D.; Blunt, M. J.; Armstrong, R. T.; Mostaghimi, P. Deep Learning in Pore Scale Imaging and Modeling. Earth Sci. Rev. 2021, 215, 103555, DOI: 10.1016/j.earscirev.2021.103555[Crossref], Google Scholar

84Latt, J.; Malaspinas, O.; Kontaxakis, D.; Parmigiani, A.; Lagrava, D.; Brogi, F.; Belgacem, M. B.; Thorimbert, Y.; Leclaire, S.; Li, S.; Marson, F.; Lemus, J.; Kotsalos, C.; Conradin, R.; Coreixas, C.; Petkantchin, R.; Raynaud, F.; Beny, J.; Chopard, B. Palabos: Parallel Lattice Boltzmann Solver. Comput. Math. Appl. 2021, 81, 334– 350, DOI: 10.1016/j.camwa.2020.03.022[Crossref], Google Scholar

85Parkhurst, D. L. User’s Guide to PHREEQC, a Computer Program for Speciation, Reaction-Path, Advective-Transport, and Inverse Geochemical Calculations; U.S. Geological Survey, 1995.Google Scholar

86Bhatnagar, P. L.; Gross, E. P.; Krook, M. A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems. Phys. Rev. 1954, 94, 511– 525, DOI: 10.1103/physrev.94.511[Crossref], [CAS], Google Scholar

87Martys, N. N.; Douglas, J. J. Critical Properties and Phase Separation in Lattice Boltzmann Fluid Mixtures. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2001, 63, 031205, DOI: 10.1103/physreve.63.031205[Crossref], Google Scholar

88Qian, Y. H.; D’Humières, D.; Lallemand, P. Lattice Bgk Models for Navier-Stokes Equation. EPL 1992, 17, 479– 484, DOI: 10.1209/0295-5075/17/6/001[Crossref], Google Scholar

89Martys, N. S.; Chen, H. Simulation of Multicomponent Fluids in Complex Three-Dimensional Geometries by the Lattice Boltzmann Method. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1996, 53, 743– 750, DOI: 10.1103/physreve.53.743[Crossref], [PubMed], [CAS], Google Scholar

90Timm, K.; Halim, K.; Alexandr, K.; Orest, S.; Goncalo, S.; Erlend, M. V. The Lattice Boltzmann Method Principles and Practice; Springer, 2017.Google Scholar

91Huber, C.; Parmigiani, A.; Chopard, B.; Manga, M.; Bachmann, O. Lattice Boltzmann Model for Melting with Natural Convection. Int. J. Heat Fluid Flow 2008, 29, 1469– 1480, DOI: 10.1016/j.ijheatfluidflow.2008.05.002[Crossref], [CAS], Google Scholar

92Parmigiani, A. Lattice Boltzmann Calculations of Reactive Multiphase Flows in Porous Media, Thesis, University of Geneva, 2011; Vol. 129.Google Scholar

93Fazeli, H.; Patel, R. A.; Ellis, B. R.; Hellevang, H. Three-Dimensional Pore-Scale Modeling of Fracture Evolution in Heterogeneous Carbonate Caprock Subjected to CO2-Enriched Brine. Environ. Sci. Technol. 2019, 53, 4630– 4639, DOI: 10.1021/acs.est.8b05653[ACS Full Text 📷], [CAS], Google Scholar

94Fazeli, H.; Masoudi, M.; Patel, R. A.; Aagaard, P.; Hellevang, H. Pore-Scale Modeling of Nucleation and Growth in Porous Media. ACS Earth Space Chem. 2020, 4, 249– 260, DOI: 10.1021/acsearthspacechem.9b00290[ACS Full Text 📷], [CAS], Google Scholar

95Patel, R.; Perko, J.; Jacques, D.; de Schutter, G.; Ye, G.; van Breugel, K. Lattice Boltzmann Based Multicomponent Reactive Transport Model Coupled with Geochemical Solver for Scale Simulations. Computational Methods for Coupled Problems in Science and Engineering , 2013; pp 806– 817.Google Scholar

96Patel, R. A.; Perko, J.; Jacques, D.; de Schutter, G.; van Breugel, K.; Ye, G. A Versatile Pore-Scale Multicomponent Reactive Transport Approach Based on Lattice Boltzmann Method: Application to Portlandite Dissolution. Phys. Chem. Earth 2014, 70–71, 127– 137, DOI: 10.1016/j.pce.2014.03.001[Crossref], Google Scholar

97Yoon, H.; Kang, Q.; Valocchi, A. J. 12. Lattice Boltzmann-Based Approaches for Pore-Scale Reactive Transport. Rev. Mineral. Geochem. 2015, 80, 393– 432, DOI: 10.1515/9781501502071-012[Crossref], Google Scholar

98Latt, J.; Malaspinas, O.; Kontaxakis, D.; Parmigiani, A.; Lagrava, D.; Brogi, F.; Belgacem, M. B.; Thorimbert, Y.; Leclaire, S.; Li, S.; Marson, F.; Lemus, J.; Kotsalos, C.; Conradin, R.; Coreixas, C.; Petkantchin, R.; Raynaud, F.; Beny, J.; Chopard, B. Palabos: Parallel Lattice Boltzmann Solver. Computers & Mathematics with Applications , 2020.[Crossref], Google Scholar

99Tan, J.; Sinno, T. R.; Diamond, S. L. A parallel fluid-solid coupling model using LAMMPS and Palabos based on the immersed boundary method. J. Comput. Sci. 2018, 25, 89– 100, DOI: 10.1016/j.jocs.2018.02.006[Crossref], [PubMed], [CAS], Google Scholar

100Kotsalos, C.; Latt, J.; Chopard, B. Bridging the Computational Gap between Mesoscopic and Continuum Modeling of Red Blood Cells for Fully Resolved Blood Flow. J. Comput. Phys. 2019, 398, 108905, DOI: 10.1016/j.jcp.2019.108905[Crossref], Google Scholar

101Kotsalosa, C.; Latt, J.; Beny, J.; Chopard, B. Digital Blood in Massively Parallel CPU/GPU Systems for the Study of Platelet Transport. Interface Focus 2019, 11, 20190116, DOI: 10.1098/rsfs.2019.0116[Crossref], Google Scholar

102Parkhurst, D. L.; Appelo, C. A. J. Description of Input and Examples for PHREEQC Version 3 ─ A Computer Program for Speciation, Batch-Reaction , One-Dimensional Transport , and Inverse Geochemical Calculations, U.S. Geological Survey Techniques and Methods, Book 6, Chapter A43, 2013, 6-43A, p 497; U.S. Geological Survey, 2013[Crossref], Google Scholar

103Tetteh, J. T.; Alimoradi, S.; Brady, P. v.; Barati Ghahfarokhi, R. Electrokinetics at Calcite-Rich Limestone Surface: Understanding the Role of Ions in Modified Salinity Waterflooding. J. Mol. Liq. 2020, 297, 111868, DOI: 10.1016/j.molliq.2019.111868[Crossref], [CAS], Google Scholar

104Lutzenkirchen, J. Surface Complexation Modelling; Elsevier, 2006.[Crossref], Google Scholar

105Tetteh, J. T.; Alimoradi, S.; Brady, P. V.; Barati, R. G. Electrokinetics at Calcite-Rich Limestone Surface: Understanding the Role of Ions in Modified Salinity Waterflooding. J. Mol. Liq. 2020, 297, 111868, DOI: 10.1016/j.molliq.2019.111868[Crossref], [CAS], Google Scholar

106Tetteh, J. T.; Alimoradi, S.; Brady, P. v.; Barati, R. G. Electrokinetics at Calcite-Rich Limestone Surface: Understanding the Role of Ions in Modified Salinity Waterflooding. J. Mol. Liq. 2020, 297, 111868, DOI: 10.1016/j.molliq.2019.111868[Crossref], [CAS], Google Scholar

107Molins, S.; Soulaine, C.; Prasianakis, N. I.; Abbasi, A.; Poncet, P.; Ladd, A. J. C.; Starchenko, V.; Roman, S.; Trebotich, D.; Tchelepi, H. A.; Steefel, C. I. Simulation of Mineral Dissolution at the Pore Scale with Evolving Fluid-Solid Interfaces: Review of Approaches and Benchmark Problem Set. Comput. Geosci. 2020, 25, 1285, DOI: 10.1007/s10596-019-09903-x[Crossref], Google Scholar

108Tetteh, J. T.; Barati, R. Crude-Oil/Brine Interaction as a Recovery Mechanism for Low-Salinity Waterflooding of Carbonate Reservoirs. SPE Reservoir Eval. Eng. 2019, 22, 877, DOI: 10.2118/194006-pa[Crossref], [CAS], Google Scholar

109Tetteh, J. T.; Alimoradi, S.; Brady, P. V.; Barati, R. G. Electrokinetics at the Limestone Surface: Understanding the Role of Ions in Modified Salinity Waterflooding. J. Mol. Liq. 2020, 297, 111868, DOI: 10.1016/j.molliq.2019.111868[Crossref], [CAS], Google Scholar

110Hiorth, A.; Cathles, L. M.; Madland, M. V. The Impact of Pore Water Chemistry on Carbonate Surface Charge and Oil Wettability. Transp. Porous Media 2010, 85, 1– 21, DOI: 10.1007/s11242-010-9543-6[Crossref], [CAS], Google Scholar

111Brady, P. V.; Krumhansl, J. L.; Mariner, P. E. Surface Complexation Modeling for Improved Oil Recovery. In SPE Improved Oil Recovery Symposium; Society of Petroleum Engineers, 2012; pp 14– 18.[Crossref], Google Scholar

112Mahani, H.; Keya, A. L.; Berg, S.; Nasralla, R. Electrokinetics of Carbonate/Brine Interface in Low-Salinity Waterflooding: Effect of Brine Salinity, Composition, Rock Type, and PH on Zeta-Potential and a Surface-Complexation Model. SPE J. 2017, 22, 053– 068, DOI: 10.2118/181745-pa[Crossref], [CAS], Google Scholar

113Tagavifar, M.; Jang, S. H.; Sharma, H.; Wang, D.; Chang, L. Y.; Mohanty, K.; Pope, G. A. Effect of PH on Adsorption of Anionic Surfactants on Limestone: Experimental Study and Surface Complexation Modeling. Colloids Surf., A 2018, 538, 549– 558, DOI: 10.1016/j.colsurfa.2017.11.050[Crossref], [CAS], Google Scholar

114Sanaei, A.; Tavassoli, S.; Sepehrnoori, K. Investigation of Modified Water Chemistry for Improved Oil Recovery: Application of DLVO Theory and Surface Complexation Model. Colloids Surf., A 2019, 574, 131– 145, DOI: 10.1016/j.colsurfa.2019.04.075[Crossref], [CAS], Google Scholar

Cited By

ARTICLE SECTIONS

Jump To

This article has not yet been cited by other publications.

Download PDF

📷 Abstract Figure 1 📷Figure 1. Schematic plot of the modeling framework. The framework has six major components, including initialization of the domain, fluid, and reactants based on both solvers; advection–diffusion process handled by PALABOS; coupling interface that connects the two solvers; chemical reactions handled by PHREEQCRM; electrical property calculations using PHREEQCRM; and model evaluations and visualizations. The simulation integrates all components dynamically in a loop shown in the figure. Figure 2 📷Figure 2. Lattice arrangement and velocity vectors for the D2Q9 model. Figure 3 📷Figure 3. Schematic plot for the evolution of the BGK LBM model within a single time step. The orange box demonstrates the optional step, whereas the blue box indicates the necessary step during the evolution, modified from ref (90). Figure 4 📷Figure 4. Illustration of the workflow for NN-assisted SCM model parameter optimization process. Figure 5 📷Figure 5. Model setup for the LBM single-phase fluid flow in a 2D channel. Figure 6 📷Figure 6. Lattice density map generated from the given pressure gradient. Figure 7 📷Figure 7. Lattice velocity distribution induced by the pressure gradient of the fluid flow channel. Figure 8 📷Figure 8. LBM simulated velocities (red dots) against the analytical solutions (blue curve). Figure 9 📷Figure 9. Model setup for the LBM diffusion validation in a 2D channel. Figure 10 📷Figure 10. Cross section (dashed line) at the center (y = (1/2)*H) of the 2D channel for concentration data extraction, where H is the width of the channel along the y-axis. Figure 11 📷Figure 11. 2D concentration distribution from LBM simulation at t = 20, 50, 100, 150, 200, and 250 s. Figure 12 📷Figure 12. Concentration profile comparison between LBM simulation and analytical solutions at t = 1, 20, 50, 100, 150, 200, and 250 s. Figure 13 📷Figure 13. 2D simulation configurations of this numerical experiment. The flow channel is filled with HCl solution, while it continuously flows through the channel with a constant flow rate from the inlet (left) to outlet (right). The calcite grain is placed at the center of the domain interacting with the solution. Upper and lower walls are impermeable and non-reactive boundaries. This figure is re-generated based on ref (107). Figure 14 📷Figure 14. (A) Velocity distribution under the steady-state condition; (B) 2D Ca2+ concentration distribution in the steady state; the higher Ca2+ concentration region surrounding the grain indicates the stable dissolution of CaCO3 on the grain surface. Figure 15 📷Figure 15. (A) 2D H+ concentration distribution map; the color contour illustrates the beginning stage (t = 0.2 s) of the consumption of H+ due to calcite dissolution kinetic reaction on the solid surface. (B) 2D H+ concentration distribution map in the steady state (t ≥ 1.5 s). Figure 16 📷Figure 16. Averaged effluent H+ concentration evolution from this study and the other five codes. (107) Figure 17 📷Figure 17. Averaged reaction/dissolution rate evolution from this study and the other five codes. (107) Figure 18 📷Figure 18. 2D PALABOS–PHREEQCRM SCM model validation scheme. Figure 19 📷Figure 19. MSE loss decay during the optimization for k values of the brine/calcite interface SCM under 25 °C (left) and 40 °C (right) conditions. Figure 20 📷Figure 20. Evolution of the optimized brine/calcite ζSCM for FWS (A) and LSW (B) under 25 °C conditions. Figure 21 📷Figure 21. Comparison of ζEXP and ζSIM based on the reaction constants from ANN optimization and Tetteh et al. (105) regarding the five single salt solution cases. Figure 22 📷Figure 22. Comparison of ζEXP and ζSIM based on the reaction constants from ANN optimization and Tetteh et al. (105) regarding the FWS/SWS/LSW cases.

CAN Fuzzy Network- CANFIS be used to simulate groundwater quality?

What software can i use for mapping?

Can someone help me with the most appropriate tools to analyse a project on innovation s within tourism and Hospitality industry?

Can you measure ownership structure yearly?

How can a random sample of public bus transport users be generated when a sampling frame is not feasible?

How can I improve upon my data so as to get a normal distribution?

I would like to identify a promoter for a gene. Given a super-scaffold containing the gene. How would I go about this and what tools are available?

What is the best way to keep an anaerobic chamber anoxic in the absence of hydrogen?

How the parameters (such as Mw, Chemical bonds, polymer bland and so on) affect the mechanic property of polymer (such as storage, loss modulus) ?

What is the characteristic length to determine the Reynolds number for the scattered turbulent fluid flow?

Joint behaviour in Abaqus is not captured well?

Issues with Permeabilization during Flow Cytometry?

How to modify Johnson cook equation in ABAQUS?

How to find the right number of integration points in second order 3d finite elements?

Can Ethovision store the generated data and video recordings directly to NAS system?

What could be the reason for not having further low chi square value (better agreement) in my refinement result?

Why do the air flow results vary totally when adding scalar in CFD ?