The Challenge
===========

Physical transport in a heterogeneous porous media and non-uniform distributions of chemical and biological species in the natural subsurface environment exert scaling behavior of reactive transport processes. In most cases the fully characterization of the system is impossible because of the inherent random distribution of the material properties. This has led to the development of reactive transport models across a wide range of spatial scales. Continuum scale models (e.g. Darcy equation for flow) and spatial averaging, smooth out any heterogeneity of a scale comparable or smaller than the representative elementary volume (grid-spacing). The advances in both computational power and the development of new numerical algorithms, have made possible pore-scale calculations including advection and diffusion of reactant in a solid matrix with a complex topology. Pore-scale models are capable of simulating complex interplay and non-linear coupling between surface reactions, flow and diffusive transport as well as ability to solve for the distribution of reactant inside pores. One of such numerical paradigm is the lattice Boltzmann method (LBM) that was developed as an extension to cellular automata models for fluid flows [1,2,3].

Lattice Boltzmann Method
===================
The LBM approaches mass and momentum conservation equations from a kinetic theory standpoint, where the macroscopic variables of interest, the local density (or pressure with an equation of state) and momentum of the fluid are retrieved from the statistical moments of particle distribution functions that follow the evolution of a discrete set of Boltzmann kinetic equations (Figure 1).

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/Lattice_Boltzmann_Phreeqc_1_6e78b3f9e4.png" alt=""/>
 <figcaption>Figure 1. Lattice topologies and lattice velocity indexing for the D2Q9 (a) and D2Q5 (b). (c) Schematic illustration of the bounce-back boundary condition at fluid–solid interfaces. The real and effective (numerical) boundaries are shifted by half a node for a straight interface (1st order accuracy) lattices. (d) Phase-field approach to surface reactions. Interface fluid nodes are considered for precipitation and interface solid nodes for dissolution. We use two flag variables, the first one is which corresponds to the type of nodes (solid or fluid dominated) and the second ζ which is null for all site but interface nodes where they can be -1 depending on there location with respect to the fluid–solid interface [2].</figcaption>
</figure>

A pore-scale reactive transport model based on LBM has been previously developed to study heterogeneous dissolution-precipitation and sorption reactions in a complex porous media [1,2] (Figure 2). The model was successfully validated with analytical solutions of various 1D and 2D reactive transport problems. It uses recent advances in the development of optimal advection-diffusion solvers within the lattice Boltzmann method framework. Another novelty of the model is that the treatment of surface reaction is independent of the fluid-rock interface topology and orientation.

<figure className="image">
<iframe width="560" height="315" src="https://www.youtube.com/embed/aDrscAjmYO8" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
 <figcaption>Figure 2. Simulating arsenic sorption on various iron minerals (i.e. ferrihydrite, goethite and magnetite). Panel (a) depicts 2D distribution of disolved arsenic in the pore space while panels (b), (c) and (d) illustrates adsorbed arsenic concentrations on ferrihydrite, goethite and magnetite, respectively. Concentrations have been normalized with respect to their inital values.</figcaption>
</figure>

The new coupled LBM-PhreeqcRM model (LAMP) uses flow and transport solver based on the Lattice Boltzmann method with the chemistry solver PhreeqcRM in the pore-scale system. PhreeqcRM is a newly-developed geochemical reaction module designed specifically to perform equilibrium and kinetic reaction calculations for reactive transport simulators that use an operator-splitting approach [4]. It comprises methods to transfer concentrations and other model properties, and retrieve selected outputs. Coupling with PhreeqcRM enables the existing multicomponent pore-scale transport model to be extended to simulate wide range of geochemical reactions.

The coupling procedure includes: initialization of the domain for initial and boundary conditions with PHREEQC; transferring the data to the LBM-based flow and transport solver; performing transport step; transferring the data to PHREEQC and performing the reactions; and update porous media (if dissolution-precipitations occurs) and sending back to transport model until the final convergence and time step is reached.

Results and Discussion
==================
To test the accuracy of the coupling procedure, two verification problems were solved using LAMP: 1) chloride tracer test and 2) ion-exchange. Both scenarios were run in a 2D micro channel and results were compared with the ones obtained from a 1D Darcy-scale PHREEQC model. In the first example a boundary water with a chloride concentration of 0.001 mol/kgw is injected at the inlet of the domain, progressively displacing the initial water with zero concentration of chloride (Figure 3). The cross-section chloride concentration profile resembles the classis parabolic velocity profile in a 2D channel i.e. higher Cl concentration along the centerline of the channel while decreasing as it approaches to the walls. As water continues infiltrating from the inlet, peak of the profile approaches to Cl concentration of 0.001 mol/kgw and once it reaches to the maximum value, the sharp parabolic front flattens as Cl concentrations gradually increase towards the vicinity of the walls.

In the second example, similar physical setup is used while it is assumed that ion-exchange reactions take place along the walls of the channel (Figure 3a). Initially the medium contains an equilibrated potassium-sodium-nitrate solution with cation exchanger. The initial water, then is flushed with boundary water composed of calcium-chloride solution. In order to assess the effect of cation exchanger on the composition of the effluent, two points along the cross section of the channel is chosen: one in the middle of the cross section of the channel and the other in the vicinity of boundary wall. Close to the longitudinal boundary, results show that sodium initially present in the column exchanges with the incoming calcium and its concentration reaches to a fixed value as long as the exchanger contains sodium. In the centerline point, however, sodium is flushed out of the system as cation exchange reaction doesn’t take place. Potassium, close to the wall, is released after sodium because of larger K value in the exchange reaction and reaches a constant concentration after an initial increase. Similar to sodium, potassium is completely displaced by the boundary water along the centerline.

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/Lattice_Boltzmann_Phreeqc_3_2db7e086d6.png" alt=""/>
 <figcaption>Figure 3. (a) Micro-channel geometry used for chloride transport and cation-exchange scenarios. In the cation-exchange problem it was assumed that exchanger exists along the longitudinal boundary walls. (b) Chloride front at the cross section of the channel (shown with the red line).</figcaption>
</figure>

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/Lattice_Boltzmann_Phreeqc_4_699a1424c8.png" alt=""/>
 <figcaption>Figure 4. Temporal evolution of the aqueous species concentrations at points close to the boundary walls (where cation exchangers exist) and at the centerline.</figcaption>
</figure>

Acknowledgments
==================
The LBM transport solver was developed at Georgia Institute of Technology in collaboration with Dr. Christian Huber, and Dr. Andrea Parmigiani. Their contributions are gratefully acknowledged.

References:
==================
[1] Shafei, B., Huber, C., Parmigiani, A., 2013. Pore-scale simulation of calcite dissolution and precipitation using lattice Boltzmann method. Mineralogical Magezine 77 (5).

 [2] Huber, C., Shafei, B., Parmigiani, A., 2014. A new pore-scale model for linear and non-linear heterogeneous dissolution and precipitation. Geochimica et Cosmochimica Acta 124, 109-130.

[3] Kang, Q.J., Chen, L., Valocchi, A.J., Viswanathan, H.S., 2014. Pore-scale study of dissolution-induced changes in permeability and porosity of porous media. Journal of Hydrology 517, 1049-1055.

 [4] Parkhurst, D., Wissmeier, L., 2015. PhreeqcRM: A reaction module for transport simulators based on the geochemical model PHREEQC. Advances in Water Resources 83, 176-189.

large_Lattice_Boltzmann_Phreeqc_1.png

medium_Lattice_Boltzmann_Phreeqc_1.png

small_Lattice_Boltzmann_Phreeqc_1.png

thumbnail_Lattice_Boltzmann_Phreeqc_1.png

Lattice_Boltzmann_Phreeqc_1.png

In this post, Reza Gooya (University of Copenhagen), Babak Shafei and Fidel Grandia, blog about a recently-developed coupled pore-scale lattice Boltzmann method with PhreeqcRM (LAMP).

Pore-scale

A new coupled pore-scale LAttice Boltzmann Method with PhreeqcRM (LAMP)