Pore network modeling Variablysaturated flow and multicomponent reactive adsorptive transport 
Reactive Solute Transport
PoreFlow performs multicomponent 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 fluidfluid (i.e., nonwettingwetting) interfaces (shown by NW) in drained pores as well as to solidfluid (i.e., solidwetting) 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 c_{CU,}_{i} is the concentration in corner unit i with volume V_{CU,}_{i}. q_{ij,k} and q_{i,n} are volumetric flow rates in k^{th} 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 N^{ij}_{edge} edges of N^{throat}_{i}_{n} throats with flow towards the corner unit i. The second term on the r.h.s. accounts for the mass influx from N^{CU}^{;}^{I}_{in}_{;}_{ed}_{ge} 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 Q_{CU}_{;}_{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. s^{sw}_{CU}_{;}_{i }[ML^{}^{3}] and s^{nw}_{CU}_{;}_{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 s^{sw}_{CU}_{,}_{i }and s^{nw}_{CU}_{,}_{i }. Here, we assume linear equilibrium adsorption at both SW and NW interfaces, such that: where k^{sw}_{d,i }and k^{nw}_{d,i} [L] are pore scale adsorption distribution coefficients at SW and NW interfaces, respectively. The specific surface area, α^{α}^{w}_{CU,i}, is defined as where A^{α}^{w}_{CU,I }is the total area of the α interface within corner unit i (where α = s, n). Next, the mass balance equation for the k^{th} edge element of a drained pore throat may be written as (assuming that corner unit j is the upstream node): where V_{ij}_{,}_{k}_{, }q_{ij}_{,}_{k}_{,} and c_{ij}_{,}_{k} are the volume, volumetric flow rate, and concentration of k^{th} edge of pore throat ij, respectively. s^{sw}_{ij}_{,}_{k} and s^{nw}_{ij}_{,}_{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, k^{sw}_{d,ij,k }and k^{nw}_{d,i j,k} [L] are pore scale adsorption distribution coefficients at SW and NW interfaces, respectively. The specific surface area, a^{α}^{w}_{ij,k}, is defines as, where A^{α}^{w}_{ij,k} is the total area of the appropriate αw interface (where α = s, n) within the k^{th} 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 N^{ij}_{edge}_{ }pore throats as well as N^{CU,i}_{edge} corner units connected to it. D is the molecular diffusion coefficient and A_{ij,k} is the cross section of k^{th} edge of tube ij, which is a function of saturation. A_{i,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 R_{CU},_{i} and R_{ij,k} represent the reaction terms for corner units of drained pore body i and k^{th} edge of drained pore throat ij, respectively. Note that for multicomponent transport these parameters should be rewritten as R^{l}_{CU,i} and R^{l}_{ij,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 nonlinear 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 subsurface 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., R_{CU},_{i} and R_{ij,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 N_{reaction} is the number of reactions, r^{l} represents, for kinetic reactions, the rate of the reaction m and a^{l,m} is the stoichiometric coefficient of species l in the reaction m. The form of rate expression, r^{l}, 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 NewtonRaphson method to linearize the set of kinetic and equilibrium reaction equations. A sequential noniterative approach is implemented, in which, first the transport step, and then the reactive step is solved in each time step of the simulation.
References Regnier, P., O’Kane, J., Steefel, C., and Vanderborght, J. (2002). Modeling complex multicomponent reactivetransport systems: towards a simulation environment based on the concept of a knowledge base. Applied Mathematical Modelling, 26(9):913 – 927.vt35

Pore network modeling Variablysaturated flow and multicomponent reactive adsorptive transport 