Pore network modeling

        Variably-saturated flow and multi-component reactive adsorptive transport

Application Examples


1- Non-reactive transport; saturated conditions

We consider the case of a pulse input of the solute into a 3D pore network domain. Averaging over the network cross section results in a 1D macro-scale concentration field. Solute transport through such a column is modeled using Advection Dispersion Equation (ADE),

where ϑ is porosity, v is average pore velocity, and DL is the longitudinal dispersion coefficient. In this equation, porosity and v values are determined directly from the pore network. To verify the pore network model, we compare the resulting breakthrough curve from the pore network model with the analytical solution for ADE for the same value of porosity, average velocity and the domain length as in the pore network. We use a network of size Ni = 40; Nj = Nk = 20 with an average velocity of 0.72 m/d.


Properties of Multi-Directional Pore Network

A pulse tracer with the concentration of 1.0 is injected through the inlet boundary for 1.8 pore volumes. At various time steps concentrations at the outlet face of the network are averaged over the cross section to obtain the breakthrough curve.

Comparison of BTC obtained from the pore network with analytical solution for a pulse input of solute to a column. The solid line shows average concentration at the outlet of the network and the symbols are analytical solution of continuum scale ADE.


2- Non-reactive transport; variably saturated conditions

An important parameter under variably saturated conditions is solute dispersion coefficient, which depends not only on the Darcy velocity, but also on degree of saturation (Toride et al., 2003). The degree of spreading of solutes is related to distribution of the pore velocities within the pore spaces. A decrease in saturation causes an increase in pore velocity fluctuations, which results in a higher solute dispersion coefficient. Using PoreFlow, we can simulate transport of solutes under different saturation levels. Since, resistance to the flow is included within drained pore bodies and pore throats, different fluxes are calculated for different edges of each drained pore, which are used in the calculation of solute transport. BTCs show a greater spreading of solute under lower saturations.

BTCs obtained from the pore network at different saturations: Sw = 1.0, 0.8, and 0.5. BTCs obtained under lower saturations show more spreading compared with BTC of tracer for saturated condition.


3- Adsorptive transport; variably saturated conditions

Adsorptive transport, such as transport of viruses and colloids, is of great importance in studies of porous media. Under variably saturated conditions, principal interactions usually occur not only at the solid-wetting (SW) interfaces, but also at the nonwetting-wetting (NW) interfaces. These interactions are greatly influenced by the water content. Concentration of adsorbed solutes depends on specific interfacial areas, which, in turn, are functions of saturation. This information will be used to calculate fluid flow and solute adsorption at interfaces and to determine the BTC of adsorptive transport. Averaging over the network, PoreFlow calculates interfacial areas as a function of the average saturation of the pore network together with capillary pressure-saturation and relative permeability-saturation curves.

(a) Capillary pressure-saturation curve, (b) relative permeability-saturation curve, (c) the normalized total interfacial area of solid-wetting phases, SW, and nonwetting-wetting phases, NW, as a function of average wetting saturation, Sw (the areas are normalized to the total area under saturated conditions, which only belongs to the solid-wetting surfaces), and (d) the resulting BTC from the network under variably saturated conditions. The BTC belong to kswd = knwd = 0.1 mm and Sw=0.5. Interfaces associated with the fluid films are neglected.


The slope of kr-Sw is larger at higher saturation. This is due to the fact that the nonwetting phase first occupies the larger pores (which have the most contribution to the conductivity of the phase), after which there is a considerable drop in conductivity of the wetting phase. Further decrease in saturation will only cause slight decrease in permeability of the wetting phase, since almost all larger and percolating pores are already drained. As shown, the SW interfacial area decrease during drainage (interfaces associated with the fluid films are neglected). This is due to the invasion of pores by the nonwetting phase. Initially, under saturated conditions, there is no NW interface present. Thus, the area of NW interface starts at zero, and increases with decrease in saturation. At any given average saturation, the fluids distribution within the pores and thus the amount of SW and NW interfaces in each pore are known. The system of mass balance equations for pore elements can be solved, and breakthrough curve can be obtained for a  given average saturation.


4– Reactive transport in wellbore cement under CO2 storage conditions

Wellbore integrity is widely regarded as one of the main potential risks regarding containment of carbon dioxide in geological storage systems. One of the primary concerns is that chemical degradation of the cement used to seal wellbores may provide or enhance preferential pathways for CO2 escape. Numerous laboratory studies have addressed the chemical interaction between cured cement and carbon dioxide. The process involves dissolution of CO2 in the pore fluid, which leads to the formation and dissociation of carbonic acid. This acidification is buffered primarily by reaction with portlandite, Ca(OH)2, and the calcium silicate hydrate, C-S-H, phases present in cement, thereby precipitating calcium carbonate, CaCO3, and eventually amorphous silica gel.

In order to integrate, explain or predict the highly variable results on cement degradation obtained in laboratory experiments or field situations, it is necessary to incorporate the key physical, chemical and microstructural processes that operate into numerical models that can be applied to any set of boundary conditions. PoreFlow is used to show the dynamics of reactive solute transport and pore structure evolution in wellbore cement under conditions of full hydro-chemical coupling. We study cement degradation reactions and the evolution of pore space when cement is exposed to CO2 at temperatures similar to those expected in CO2 storage sites. Solute concentration values for both pore bodies and pore throats within the pore network are calculated. To calculate concentrations changes due to chemical reactions Biogeochemical Reaction Network Simulator (BRNS) (Regnier et al., 2002) was used. BRNS is capable of handling a comprehensive suite of multi-component complexation and mineral precipitation and dissolution reactions as well as reaction networks characterised by multiple kinetic pathways.

A constant CO2 molar concentration equal to 2.0M equilibrated with water at 50°C is imposed at the inlet boundary of the model domain using a Dirichlet boundary condition.

Those reactions involving only aqueous species are treated as attaining equilibrium, whereas the dissolution-precipitation reactions of calcite and portlandite are incorporated in the model in the form of kinetic processes.


pH and molar chemical concentration profiles for different carbonate components at 4 different times (1, 6, 12 and 60 months) (after Raoof et al., 2012).


The dynamics of cement degradation involves interaction between the kinetics of reactions and the fluxes of aqueous species which determine their availability for reactions. Results show development of different regions: (1) a zone adjacent to the inlet face, which is characterized by an increase in porosity due to extensive dissolution, (2) a carbonation zone with decreased porosity, (3) the carbonation front which made a thin layer with the lowest porosity due to calcium carbonate precipitation, and (4) dissolution zone. Similar results have been reported through laboratory observations (Rimmelé et al., 2008).


Porosity, pH profiles and carbonate speciation after 1 year of CO2 invasion into the portlandite cement and creation of low and high permeability zones. The concentrations are normalized by their maximum values (after Raoof et al., 2012).


5– Porosity-permeability evolution

Progress of chemical reactions may cause change in porosity and flow filed in reacting porous media. PoreFlow can be used to obtain the corresponding relationship between porosity and permeability. We consider dissolution of carbonate calcium, CaCO3, due to the injection of a low pH solution. This problem has various applications such as injection of CO2 into a reservoir rock. The calcium carbonate equilibria parameters together with kinetic dissolution of calcite (occurring via three parallel reactions, Plummer et al., 1978) has been used.


Calcium carbonate equilibria at 25C [Koutsoukos and Kontoyannis, 1984]


Concentrations in equilibrium with CaCO3 is used as the initial conditions. Subsequent injection of solute with pH of 3.0 for 60 days caused the dissolution of calcite, and therefore, an increase in porosity. By averaging over the network domain at successive time steps, we calculate the evolution of porosity and permeability as chemical reactions progress. Dissolution of calcite causes an increase in pore throat sizes which have the most effect on permeability of the pore structure.


Initial distribution of pore throats together with pore throat size distribution after the dissolution reaction, b) relation between permeability increase with porosity. k0 = 1.2x10-12 m2 is the initial permeability




· P. Koutsoukos and C. Kontoyannis (1984). Precipitation of calcium carbonate in aqueous solutions. J. Chem. Soc., Faraday Trans. 1, 80(5):1181–1192.

· Plummer, L., Wigley, T., and Parkhurst, D. (1978). The kinetics of calcite dissolution in c 02-water systems at 5 deg to 60 deg c and 0. 0 to 1.0 atm c 02. American Journal of Science, 278.

· A. Raoof, H.M. Nick, T.K.T. Wolterbeek, C.J. Spiers (1984). Pore-scale modeling of reactive transport in wellbore cement under CO2 storage conditions. International Journal of Greenhouse Gas Control 11S x S67–S77

· G. Rimmel´e, V. Barlet-Gou´edard, O. Porcherie, B. Goff´e, and F. Brunet (2008). Heterogeneous porosity distribution in Portland cement exposed to co2-rich fluids. Cement and Concrete Research, 38(8-9):1038–1048

· N. Toride, M. Inoue, and F.J. Leij (2003). Hydrodynamic dispersion in an unsaturated dune sand. Soil Science Society of America Journal, 67(3):703–712

































































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