🔬 condensed matter

A Conservative Darcy-Brinkman-Stokes Formulation for PHREEQC-Coupled Reactive Transport with Cell-Based Surface Complexation

This paper presents a mass-conservative Darcy-Brinkman-Stokes formulation implemented in GeoChemFoam and coupled with PHREEQC, enabling accurate reactive transport simulations with cell-based surface complexation in both resolved and unresolved porous media structures.

Original authors: Hannah P. Menke, Florian Doster, Julien Maes

Published 2026-10-15

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Hannah P. Menke, Florian Doster, Julien Maes

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Deep underground, water does not simply flow through empty tunnels like cars on a highway. Instead, it moves through a labyrinth of rock that is porous, meaning it is filled with tiny, interconnected holes and cracks. In many of these underground environments, the rock is not just a passive container; its surfaces are chemically active. As water moves, dissolved minerals and contaminants can stick to these surfaces, a process that slows their travel and changes their concentration. This interaction is critical for understanding how pollutants move through groundwater, how we might capture valuable minerals like lithium from brine, and how radioactive waste can be safely stored for thousands of years. The challenge for scientists is that the rock's internal structure is incredibly complex. Some parts have large, open channels where water rushes through, while other parts are made of fine, sponge-like material with microscopic pores that hold a vast amount of internal surface area. Traditional computer models often struggle to capture both the fast flow in the large channels and the slow, sticky chemistry happening inside the tiny, invisible pores simultaneously.

A team of researchers has developed a new way to simulate these complex interactions, allowing them to model both the flow of water and the chemical sticking of particles on a single, unified grid. Their approach treats the rock not as a collection of separate large holes and tiny sponges, but as a continuous material where the amount of water and the amount of reactive surface can vary from one tiny block of the model to the next. In this new method, the computer divides the rock image into small cells. For each cell, the model calculates how much water is actually moving through it and how much surface area is available for chemicals to stick to, even if that surface area is too small to be seen in the original image. By keeping track of the total mass of the chemicals—both the part floating in the water and the part stuck to the rock—the researchers ensure that nothing is lost or created by mistake during the calculation. This mass conservation is vital for making accurate predictions about how long it takes for a contaminant to travel or how much of a valuable mineral can be recovered.

The researchers tested their new method against a series of rigorous benchmarks to ensure it worked correctly. They started with simple, one-dimensional scenarios where the answer could be calculated by hand, confirming that their computer model matched the expected results for how fast chemicals move and how much they slow down due to sticking. They then moved to more complex, three-dimensional images of real rock and synthetic materials. In one test, they simulated a rock with a smooth, spherical shape to see if their method could accurately measure the surface area. They found that their approach could calculate the surface area with an error of less than one percent, even when the image resolution was relatively low. In another test, they compared their new method against a much more detailed, computationally expensive simulation that explicitly mapped every single tiny pore. Their new, faster method reproduced the results of the detailed simulation with a difference of less than one percent in the total amount of lithium captured, proving that it could handle the complexity of real-world rock without needing to resolve every microscopic detail.

The team then applied their method to three distinct, real-world scenarios to see how it performed in practice. First, they modeled a fractured rock system where lead ions moved through cracks coated with iron oxide. The simulation showed how the lead slowed down as it stuck to the coating, creating a clear picture of how the contaminant would spread through the network of fractures. Second, they simulated the capture of manganese from a salty brine solution containing a large amount of magnesium. Because manganese and magnesium compete for the same sticking spots on the rock surface, the model had to calculate which one would win. The results showed how the manganese gradually displaced the magnesium, accumulating on the surface of the rock in a pattern that depended on the flow of the water. Finally, they modeled the movement of strontium through a limestone formation, a scenario relevant to the safety of deep geological repositories for nuclear waste. In this case, the strontium had to compete with calcium, which was already present in the water. The simulation successfully tracked how the strontium pushed the calcium off the rock surfaces as it moved forward, creating a distinct wave of chemical exchange.

The key finding of this work is that it is possible to model reactive transport in complex, multi-scale rock formations without needing to resolve every single microscopic pore. By treating the reactive surface area as a property of the cell itself rather than just the visible boundary of the rock, the researchers can include the effect of invisible, internal surfaces that are crucial for chemical reactions. This allows for accurate simulations of how contaminants and valuable minerals move through the subsurface, even when the rock contains a mix of large open spaces and tiny, sponge-like regions. The method is conservative, meaning it strictly follows the laws of mass balance, ensuring that the total amount of material remains constant throughout the simulation. This reliability makes it a powerful tool for predicting the long-term behavior of underground systems, from cleaning up polluted aquifers to designing safer storage for radioactive waste. The researchers demonstrated that their approach works across a wide range of flow speeds and rock types, providing a robust framework for understanding the hidden chemistry of the deep earth.

Technical Summary: A Conservative Darcy-Brinkman-Stokes Formulation for PHREEQC-Coupled Reactive Transport

Problem Statement
Reactive transport in multi-scale porous media, such as carbonates, presents a significant modeling challenge. These systems often contain open macropores alongside microporous matrix regions. While image-based simulations can resolve macropores, the internal reactive surface area of microporous phases often falls below the image resolution (sub-voxel scale). Classical Darcy-scale models rely on constitutive upscaling assumptions that cannot capture spatial heterogeneity, while fully resolved pore-scale models (e.g., explicit-boundary Darcy-Brinkman-Stokes or DBS) fail when the dominant reactive area lies inside coarse computational cells rather than on the mesh boundary. Existing PHREEQC-coupled workflows typically impose chemistry as flux conditions on resolved interfaces, leaving unresolved internal surface areas either ignored or represented by ad hoc volumetric sinks that compromise mass conservation and the physical definition of cell-scale inventories.

Methodology
The authors present a mass-conservative DBS formulation implemented in GeoChemFoam (an OpenFOAM-based toolbox) and coupled to PHREEQC via the PhreeqcRM interface. The core innovation is a cell-based surface complexation approach where the reactive surface is a property of the cell volume rather than the mesh boundary.

Governing Equations: The flow is governed by the Darcy-Brinkman-Stokes equations, where a cell-centered porosity field ( $\epsilon$ ) varies smoothly from solid to open pore. Transport is solved for the physical mobile-water fraction ( $\theta$ ), which governs aqueous storage and diffusion.
Surface Area Definition: The reactive surface area ( $a_{surf}$ $a_{s u r f}$ ) is defined as a cell-centered bulk-volume density ( $m^2 m^{-3}_{bulk}$ $m^{2} m_{b u l k}^{- 3}$ ) composed of two additive terms:
1. Interface Term ( $a_{interface}^{surf}$ ): Reconstructed from the resolved image interface between pore and solid phases using smooth-interface algorithms (e.g., occupancyGaussian).
2. Intrinsic Term ( $a_{intrinsic}^{surf}$ ): A prescribed sub-voxel contribution derived from phase-specific surface densities ( $s_k$ ) weighted by the geometric phase fraction ( $f_k$ ) within the cell. This allows the representation of internal microporous surface area that is not resolved by the image segmentation.
Conservative Coupling: The local mass inventory is defined as $M = \theta C + a_{surf} \Gamma$ $M = θ C + a_{s u r f} Γ$ , where $C$ $C$ is aqueous concentration and $\Gamma$ $Γ$ is surface concentration. The solver uses a first-order Lie splitting scheme:
1. Transport Step: Advances aqueous species conservatively using finite-volume DBS transport.
2. Chemistry Step: Maps the updated aqueous state, mobile-water fraction, and cell-based surface inventory to PHREEQC. PHREEQC solves for local thermodynamic equilibrium (surface complexation) and returns updated aqueous and adsorbed species, ensuring strict mass conservation between the aqueous and surface phases within each cell.
Implementation: The reactive surface area and site densities are treated as fixed fields during the simulation. Cells with mobile-water fractions below a threshold ( $\theta_{mask}$ ) are excluded from chemistry.

Key Contributions

Decoupling of Flow and Chemistry Volumes: The formulation distinguishes between the numerical DBS porosity (governing flow resistance) and the physical mobile-water fraction (governing chemical solution volume), though they are set equal in the saturated cases presented.
Cell-Based Unresolved Surface: The method successfully integrates unresolved internal surface area as an intrinsic volumetric density, allowing reactive transport modeling in microporous matrix regions without requiring explicit mesh resolution of the internal pore network.
Conservative Mass Balance: The approach ensures that dissolved and surface-bound species enter the same local mass balance, avoiding heuristic redistribution of source terms.

Results and Validation
The formulation was verified through a comprehensive test suite:

1D Retardation Benchmark: The coupled advection-diffusion-sorption response matched analytical solutions and finite-domain references, demonstrating correct retardation factors and first-order convergence in space and time.
Closed-System Batch Balance: Mass balance tests confirmed that the transfer of mass between aqueous and surface phases within PHREEQC preserves total inventory with residuals below $10^{-9}$ .
Surface Area Reconstruction: Smooth-interface algorithms were tested against analytic geometries (spheres and circular posts). Errors were bounded within sub-percent levels (e.g., $<0.2\%$ for 2D posts, $<1\%$ for 3D spheres at sufficient resolution) compared to body-fitted snappyHexMesh references.
Resolved-to-DBS Two-Scale Benchmark: A lithium micromodel comparison between a fully resolved pore-scale simulation and the upscaled DBS formulation showed that the cell-based intrinsic area representation reproduced both aqueous and sorbed inventories within 0.6% across Péclet numbers ( $Pe$ ) of 0.1 to 10. While the total uptake was highly accurate, the sorbed streamwise profile showed shape discrepancies up to 4.75%, attributed to the smoothing of pore-scale concentration variations.
Application Cases: The model was applied to three complex image-based scenarios:
- 2D Fracture Network: Pb $^{2+}$ retardation by hydrous ferric oxide (HFO) coatings, capturing spatial heterogeneity in strong-site occupancy.
- 3D Kelvin Cell: Competitive Mn $^{2+}$ /Mg $^{2+}$ exchange, demonstrating the model's ability to handle multi-species competition in a microporous lattice.
- 3D Estaillades Limestone: Sr $^{2+}$ /Ca $^{2+}$ exchange simulating radionuclide migration, showing chromatographic exchange patterns where Sr $^{2+}$ displaces Ca $^{2+}$ in both aqueous and sorbed phases.

Significance and Claims
The paper claims that this framework supports conservative surface-complexation modeling in both resolved and unresolved pore structures. By treating the reactive surface as a cell property rather than a mesh boundary, the method overcomes the resolution limitations of traditional pore-scale models in multi-scale media. This is particularly significant for applications where selectivity and surface accessibility control performance, such as contaminant retention, radionuclide transport in geological repositories, and critical-mineral recovery (e.g., direct lithium extraction). The authors note that the accuracy of the intrinsic area representation relies on external characterization (e.g., BET adsorption) to constrain the sub-voxel surface density, and that the current implementation assumes local equilibrium and fixed surface areas, excluding kinetic reactions and mineral dissolution/precipitation.

Technical Summary: A Conservative Darcy-Brinkman-Stokes Formulation for PHREEQC-Coupled Reactive Transport

More like this