Pore network modeling

        Variably-saturated flow and multi-component reactive adsorptive transport

Reactive Solute Transport


PoreFlow performs multi-component reactive solute transport accounting for both advection and diffusion processes. The reactions can occur within the liquid phase, as well as between the liquid and solid phases which may result in an evolution of porosity and permeability. Similar to the flow simulations we use complex formulations for more accurate modeling of transport problems in presence of the nonwetting phase. This is done by refining the discretization within drained pores and calculating dissolved solute and adsorbed mass concentrations in edges of drained pore bodies and pore throats.

There are two options for simulating solute transport: i) advection only, ii) advection and diffusion. The advective transport module is also capable of simulating linear equilibrium and kinetic adsorption under both saturated and variably saturated conditions (without the need for coupling with reaction simulator). However, when including diffusive transport of components and simulation of system of chemical reactions, Biogeochemical Reaction Network Simulator, BRNS (Regnier et al., 2002) performs the adsorption/reaction part. This gives major advantages for simulation of complex reaction systems.


Adsorptive solute transport, without diffusion

Solutes may be adsorbed to the fluid-fluid (i.e., nonwetting-wetting) interfaces (shown by NW) in drained pores as well as to solid-fluid (i.e., solid-wetting) interfaces (shown by SW). Assuming flow from a given pore corner unit, j, towards corner unit i through corners of drained pore throat ij:

where cCU,i is the concentration in corner unit i with volume VCU,i. qij,k and qi,n are volumetric flow rates in kth edge of drain pore throat ij and drained pore body i, respectively. The first term on the r.h.s. is due to the mass influx via Nijedge edges of Nthroatin throats with flow towards the corner unit i. The second term on the r.h.s. accounts for the mass influx from NCU;Iin;edge neighboring corner units (within the same pore body) with flow towards corner unit i. The third term shows the mass leaving the corner unit, where QCU;i is the total flux leaving (or entering) the corner unit i. The last two terms account for the mass adsorbed onto the solid walls of the corner unit and on the NW interface within the corner unit, respectively. sswCU;i [ML-3] and snwCU;i are the mass adsorbed to SW and NW interfaces per unit volume of the corner unit.

To close the system, we need extra equations for sswCU,i and snwCU,i . Here, we assume linear equilibrium adsorption at both SW and NW interfaces, such that:

where kswd,i and knwd,i [L] are pore scale adsorption distribution coefficients at SW and NW interfaces, respectively. The specific surface area, ααwCU,i, is defined as

where AαwCU,I is the total area of the α interface within corner unit i (where α = s, n).

Next, the mass balance equation for the kth edge element of a drained pore throat may be written as (assuming that corner unit j is the upstream node):

where Vij,k, qij,k, and cij,k are the volume, volumetric flow rate, and concentration of kth edge of pore throat ij, respectively. sswij,k and snwij,k [ML-3] are the concentrations of masses adsorbed at SW and NW interfaces, respectively. Within the pore throats we assume linear equilibrium adsorption, so that we can write:

where, kswd,ij,k and knwd,i j,k [L] are pore scale adsorption distribution coefficients at SW and NW interfaces, respectively. The specific surface area, aαwij,k, is defines as,

where Aαwij,k is the total area of the appropriate  αw interface (where α = s, n) within the kth edge of pore throat i.


Solute transport, including diffusion

For the case of solute transport through the network including both advection and diffusion, the mass balance equation for a pore body corner unit is extended to the following form:

The fourth and fifth terms show the diffusive mass flux between corner unit i and Nijedge pore throats as well as NCU,iedge corner units connected to it. D is the molecular diffusion coefficient and Aij,k is the cross section of kth edge of tube ij, which is a function of saturation. Ai,n is the cross section of the pore body edge connecting corner unit i to corner unit n.

Similarly, the mass balance equation for an edge element of a drained pore throat (assuming that corner unit j is the upstream node) is modified to:

In above equations RCU,i and Rij,k represent the reaction terms for corner units of drained pore body i and kth edge of drained pore throat ij, respectively. Note that for multi-component transport these parameters should be rewritten as RlCU,i and Rlij,k for chemical species l.


Biogeochemical Reaction Network Simulator (BRNS)

PoreFlow uses Biogeochemical Reaction Network Simulator (BRNS) (Regnier et al., 2002) for solving a set of coupled non-linear equations modeling chemical reactions between components. BRNS has been used as a flexible numerical engine applicable for a variety of reactive transport problems in surface and sub-surface environments. One of the main features of BRNS is an automated procedure for code generation using a MAPLE interface. This provides high flexibility to include and combine alternative biogeochemical process descriptions in the model. Using this interface, kinetically controlled and equilibrium reactions, and combinations of both can be specified with arbitrary form of equations describing the transformations of the chemical species. A description on how symbolic programming can be used to create the model specific source code necessary to the numerical solution of the governing equations is explained in Regnier et al. [2002].

The reactive processes, rate parameters and equilibrium constants, species involved and stoichiometry are required to define a specific reaction network and to solve for the reactive terms in each pore body and pore throats (i.e., RCU,i and Rij,k terms in mass balance equations). The reaction term represents the sum of reactions affecting concentration of a given species l can be written as:

where Nreaction is the number of reactions, rl represents, for kinetic reactions, the rate of the reaction m and al,m is the stoichiometric coefficient of species l in the reaction m. The form of rate expression, rl, is arbitrary, even nonlinear, and can be a function of several concentrations within the system. The MAPLE interface is employed to generate Fortran files, such as the function residuals, which represent the reaction network, and the Jacobian matrix, which contains the partial derivatives of the function residuals with respect to the unknown concentrations. The generated Fortran files and the solver engine for solving a nonlinear set of equations ensued from kinetic and equilibrium reactions of multi components are embedded in PoreFlow. The reactive solver uses a Newton-Raphson method to linearize the set of kinetic and equilibrium reaction equations. A sequential non-iterative approach is implemented, in which, first the transport step, and then the reactive step is solved in each time step of the simulation.



Regnier, P., O’Kane, J., Steefel, C., and Vanderborght, J. (2002). Modeling complex multi-component reactive-transport systems: towards a simulation environment based on the concept of a knowledge base. Applied Mathematical Modelling, 26(9):913 – 927.vt35


 Pore network modeling

        Variably-saturated flow and multi-component reactive adsorptive transport