Review of pore network modelling of porous media: experimental characterisations, network constructions and applications to reactive transport

Pore network models have been applied widely for simulating a variety of different physical and chemical pro- cesses,includingphaseexchange,non-Newtoniandisplacement,non-Darcy ﬂ ow,reactivetransportandthermo-dynamicallyconsistentoillayers.Therealismofsuchmodelling,i.e.thecredibilityoftheirpredictions,dependsto a largeextenton the quality of thecorrespondencebetween thepore spaceofagiven medium and the porenet-work constructed as its representation. The main experimental techniques for pore space characterisation, in- cluding direct imaging, mercury intrusion porosimetry and gas adsorption, are ﬁ rstly summarised. A review of themainporenetworkconstructiontechniquesisthenpresented.Particularfocusisgivenonhowsuchconstruc- tions are adapted to the data from experimentally characterised pore systems. Current applications of pore net-workmodels are considered, with special emphasisontheeffects ofadsorption,dissolution and precipitation,as wellasbiomassgrowth,ontransportcoef ﬁ cients.Porenetworkmodelsarefoundtobeavaluabletoolforunder-standing and predicting meso-scalephenomena, linking single poreprocesses, whereother techniquesare more accurate, and the homogenised continuum porous media, used by engineering community.


Introduction
Flow and transport phenomena in porous media play an important role in many diverse fields of science and technology. For example, radioactive waste management is one of the most pressing problems facing the world today, because of the longevity of radionuclide and the possibility of their transport to the surface environment. For long-term performance assessment of nuclear repositories, knowledge concerning the transport of radionuclide in the back-fill material is required. One porous medium, bentonite, is presently considered as the best candidate for the high level waste disposal, due to its large specific surface area, high ion-exchange capacity, and sorption affinity for organic and inorganic ions. Sorption onto bentonite plays an important role in retarding the migration of radionuclide from a waste repository (Bourg et al., 2003(Bourg et al., , 2006Bradbury and Baeyens, 2003). Another field where such knowledge is highly demanded is the construction industry. Building materials such as bricks, concrete and sandstone all are porous. These materials may interact with their environment leading to degradation of the structures with time. A specific example is that the constituent ions in the hydrated cement paste matrix may leach out from the concrete and cause the concrete structure to become weak. Degradation processes like this strongly relate to the transport phenomena in concrete. Transport mechanisms involved in the degradation processes include: permeation of water or aqueous solutions, diffusion of gaseous compounds and sorption of ionic contaminated water. In such a case, transport properties in these porous media are usually considered as indicators for evaluating the durability and ultimately predetermining the service life of these structures (Wang and Jivkov, 2015;Wang et al., 2015a,b;Wang et al., 2016;. Another example is the contaminated underground water transport, which is controlled by a variety of transport mechanisms including permeation, diffusion and dispersion (McCarthy and Zachara, 1989).
There are several challenges in analysing transport in porous media. Firstly, the intricacy of the pore structure makes the transport processes in porous media very complex. Pores tend to have irregular surfaces and some of them even make dead ends. These factors influence the flow and transport behaviour significantly. Another difficulty of studying transport processes in porous media is the evolution of pore structure during their service or operation, etc. These changes may result from chemical, electrochemical or bacterial effects as well as from mechanical damage, such as micro cracking. The evolution of the pore space may lead to changes in the macroscopic transport properties, such as permeability and diffusivity. Furthermore, the transport properties in porous media are a function of types of species in solution and the phase (Appelo et al., 2010). Finally, many reactive pore network models (Li et al., 2006;Raoof et al., 2012) employ a reaction rate or mass transfer coefficient measured in a macroscopic experiment. For example, these authors determine reaction rate coefficients from spinning disk, flowthrough reactor, or continuous fluidized bed reactor experiments. These results may not be applicable to reaction rates in individual pores, as the smaller scale magnifies the effects of mass transfer resistance and concentration gradients. Raoof and Hassanizadeh (2010b) demonstrated the importance of incorporating kinetics developed for the pore scale, rather than relying on macroscopic averages.
Due to the complexity of the pore structure and the change of the environmental conditions, a great deal of experimental, theoretical and numerical approaches have been proposed and developed to study transport processes through porous media during the past decades (Abichou et al., 2004;Aytas et al., 2009;Boult et al., 1998;Ghassemi and Pak, 2011;Kohler et al., 1996;Li et al., 1992;Ren et al., 2009;Tsai and Chen, 1995;Wang et al., 2005). Field and reactive transport experiments as well as studies involving X-ray, SEM, TEM etc. are usually expensive and time consuming. The measurements are highly sensitive to the material composition, sample preparation, methodology and testing environment. Analytical solutions are typically restricted to problems with assumed homogeneous properties and specific boundary conditions, some of which have limited practical relevance or are complicated to evaluate.
Pore-scale simulations have improved the understanding of largescale natural processes and informed large-scale geotechnical applications. Their importance comes from the fact that they can produce rather cost-effective and accurate predictions for local transport (diffusion/ permeation), and at the same time allow for systematic variations of the system's parameters (pore space geometries, fluid properties, and boundary conditions) to assess their impact, which is much more difficult to achieve than with experiments (Meakin and Tartakovsky, 2009). With pore-scale models one can make improved assessments of macroscopic transport properties by varying the pore space structure parameters. This offers a way to understand the scale dependence of continuum transport parameters. Such scale dependence cannot be captured by an effective medium Darcy approach. The pore-scale modelling is dominated by particle-based methods. These include the promising lattice Boltzmann method (Hao and Cheng, 2010; and smoothed particle hydrodynamics (Tartakovsky and Meakin, 2006;Tartakovsky et al., 2007;Zhu and Fox, 2002).
The particle methods, while suitable for pore-scale analyses, become inefficient at the meso-scale, e.g. when the system requiring analysis has tens or hundreds interconnected pores in each direction (Tartakovsky et al., 2007). Further, these methods are time consuming and only very limited pore volumes can be addressed. In order to increase the pore volume with affordable computational resources and a reduced impact on the reliability of the results, the pore network model (PNM) approach has been used to study reactive transport phenomena. At the meso-scale the classical macroscopic equations, such as the Darcy law, are yet not needed, i.e. fluid flow and solute transport processes are simulated directly in the pores, but the precise particle dynamics, which can be analysed by particle methods, is not accounted for. This is done by creating a virtual representation of the porous medium consisting of pore bodies and pore throats of different sizes (the "geometry" of the porous medium), connected to each other as required (the "topology" of the porous medium). It is then possible to simulate the fluid flow and solute transport process of interest at the meso-scale through this network, with the relevant physics implemented on a pore-to-pore basis. The first network model was constructed by Fatt (1956) who exploited the analogy between flow in porous media and a random resistor network. Afterwards, the models have grown in sophistication and now can deal with irregular lattices, wetting layer flow, arbitrary wettability and any sequence of displacement in twoand three-phase flow, as well as a variety of different physical processes, including phase exchange, non-Newtonian displacement, non-Darcy flow, reactive transport and thermodynamically consistent oil layers, etc., (Balhoff and Wheeler, 2009;Blunt et al., 2002;Lopez et al., 2003;Ryazanov et al., 2009;Yiotis et al., 2006).
Compared to the pore-scale modelling methods, such as the lattice-Boltzmann and particle methods, pore-network models take less time to compute and require less computer capacity due to the inherent simplifications of the pore-space in the construction of the model. Pore network simulations are much less computationally demanding than direct methods. This allows researchers to incorporate more heterogeneity in modelling larger rock volumes. As natural geological media can be heterogeneous on all length scales, this is an important advantage of pore network modelling. Pore network models have been used extensively to simulate multiphase and single-phase fluid flow in porous media, and these models will continue to provide important insight and information in the future. However, pore network models are usually constructed based on simplified geometry due to the lack of experimental resolution in pore space characterisation. This could make the predictions of transport properties less credible. The greatest challenge in obtaining credible results is the identification of features/phenomena relevant to the modelled process and neglecting the rest to reduce the computational load. For example, to compute the single phase flow field, it is necessary that the principal connected pore space is modelled; while the very small pores below the image resolution can be neglected as they make little contribution to the overall behaviour (Blunt et al., 2013).
The suitability of the pore-scale modelling techniques for a given application depends on the governing equations, the underlying assumptions for the pore-scale flow and transport equations, as well as the length-scales of the (computational) domain, etc. While the lower scale limit of a pore-scale technique is determined by the scale of the governing equations, the upper scale limit is set by the computational power. In this review we focus on the pore network models with comprehensive account for the methods for obtaining pore space information, constructing pore networks and the application of pore network models to reactive transport in porous media.

Characterisation of pore space information
For pore network model construction, the geometry and topology of the pore space are required. There are several ways to characterise the pore space. One is imaging techniques such as producing 3D images by mapping the real interior structure of original samples (the destructive approach of cutting and stacking serial 2D sections, confocal laser scanning microscopy and non-destructive X-ray micro-tomography (micro-CT)), and constructing synthetic 3D rock images from high resolution 2D thin sections using statistical methods or geological process simulation. Other non-destructive approaches are mercury intrusion porosimerty and gas adsorption, which produce pore size distribution and surface area, but not the full connectivity. The main methods for space characterisations are reviewed in this section.
To image geological porous materials at the micro-scale, three types of micro-CT systems are in common use: medical CT, industrial X-ray generation tube and synchrotron micro-tomography. They primarily differ in X-ray energy and source, means of sample manipulation and detector geometry. Although the best image resolution reported in the literature is from synchrotron micro-CT, the samples need to be relatively small to achieve these resolutions. This may result in poor statistical representation of the bulk material. Typical spatial resolution that medical CT systems can achieve is between 200 and 500 μm, industrial systems range from 50 to 100 μm and synchrotron based CT systems can reach from 1 μm to 50 μm (Wildenschild et al., 2002a). Presently, laboratory systems with genuine submicron capabilities exist, some with voxel resolution of 20 nm, providing spatial resolution of 50-60 nm. Three main configurations are used in systems that seek submicron resolution; a good review can be found in (Schlüter et al., 2014;Withers, 2007). Specific details of these imaging techniques and their evolution can be found in previous reviews (Blunt et al., 2013;Ketcham and Carlson, 2001).
It should be noted that X-ray CT may obscure significant features or misinterpret attenuation values of a single material in different image sections resulting in complicated quantitative image analysis (Ketcham and Carlson, 2001). Commonly experienced problems include beam hardening, high-frequency noise, scattered X-rays, defective detector pixels or poorly cantered samples (Ketcham and Carlson, 2001;Wildenschild et al., 2002b). In addition, there might be errors and distortions arising from the CT reconstruction (Ketcham and Carlson, 2001). These problems cannot be completely avoided, but can be alleviated by careful detector calibration, metal filters and sample centring. Application of wedge calibration and dual energy scans (Rebuffel and Dinten, 2007) can also reduce beam hardening effects, but these techniques are rarely implemented in standard scanning procedures (Ketcham and Carlson, 2001).

Focused ion beams and scanning electron microscopy
Focused ion beams (FIB) and scanning electron microscopy (SEM) are destructive imaging techniques. SEM is a useful technique for extracting two-dimensional (2D) images of the microstructures but does not provide the third spatial component of the sample, which is important to find interconnected regions and pore volumes, shapes and sizes. Depending on the instrument, the resolution achieved can be between 1 and 20 nm. The world's highest resolution conventional SEM can reach a point resolution of 0.4 nm by utilising a secondary electron detector (http://www.nanotech-now.com/news.cgi?story_id= 42612). FIB is a well-established technique to acquire very high-resolution three-dimensional images, typically just a few micro-meters across (Curtis et al., 2010;Lemmens et al., 2010;Tomutsa et al., 2007). This technique can achieve less than 1 nm imaging resolution (http:// www.fibics.com/fib/tutorials/introduction-focused-ion-beam-systems/ 4/). Although this method provides a great potential for imaging hydrocarbon bearing rocks at high resolution and produces images with better quality than the electron beam imaging due to less charging of the surface, it is still significantly time consuming due to the refocusing between milling and imaging as well as the sample repositioning (Tomutsa et al., 2007). Due to these limitations, this technique exposes only small areas of observation and cannot provide adequate sampling to characterise the sample, e.g. from a shale reservoir.
The combination of FIB and SEM (FIB-SEM) is usually used to compute micro-structural properties of porous media. This combination allows for the observation of fine macro-pores and meso-pores within the porous medium (Keller et al., 2011;Tomutsa et al., 2007). This method can typically achieve voxel dimensions of tens of nanometers, thus allowing for analysing volumes of around (10-30) 3 μm 3 within practical measuring time (Michael et al., 2007).
Multi-scale imaging capabilities, such as the development and integration of a range of experimental and computational tools, can facilitate probing the structure of materials across various scales in an integrated fashion. Given their important function, multi-scale imaging capabilities have been an area of focus of several research groups in the recent years.
The combination of FIB, SEM and TEM (transmission electron microscopy) is used to analyse the structure of pore space (Wirth, 2009). For example, the backscattered scanning electron microscopy (BSEM) and focused ion beam SEM (FIBSEM) have been utilised to analyse rocks (Sok et al., 2010). These techniques are particularly suitable for examining carbonates due to their multi-modal pore structure, which can range from 10 nm to 10 cm. Sok et al. (2010) report on both plugs to pore scale registration in 3D, as well as pore-scale to submicron scale registration of features. The plug to pore scale approach requires one or more subsets of the sample originally imaged (at the scale of 4 cm and at 20 μm/pixel) which then have to be reimaged with Micro CT at a smaller sample volume (5 mm and at 2.5 μm/pixel). After that the data should be carefully registered and integrated. The pore-scale to submicron scale registration was done by cutting the subsamples and preparing thin sections from the original field of view of the Micro CT images. The thin sections were then imaged with BSEM and registered to the Micro CT image. Finally, in order to estimate reservoir properties for complex materials, such as carbonate and mudstones, FIBSEM images were also obtained at voxel resolutions of 50 nm. A longer term goal of Sok et al. group is to undertake complete multi-scale registration from the whole core and/or plug scale through the Micro CT length scale and down to the submicron (SEM or FIBSEM) scale.

Nuclear magnetic resonance
Recent NMR development has demonstrated the use of relaxation, cryoporometry, spectroscopy, diffusion, and imaging techniques to quantify pore structure (i.e. pore sized distribution, pore morphology, connectivity etc.), fluid properties, and rock heterogeneity. The NMR relaxometry poses the least demands on magnetic field quality, NMR imaging favours magnetic fields linear in space with a constant gradient, and NMR spectroscopy, places the highest demands on field homogeneity and field stability.
Magnetic resonance imaging allows the imaging of the interior of the rocks to obtain the spatial distribution across a much larger scale (Blümich et al., 2009). Callaghan (1993) discovered the diffusion diffraction phenomenon when the gradient wavelength approaches the characteristic pore size. Since the water molecules in the porous material will move randomly, they will probe the pore structure. This will influence the transverse relaxation time. Therefore, NMR can provide information on the pore-size distribution of a porous material by measuring this transverse relaxation time. The major advantage of NMR in comparison with classical methods is the short measurement time. This would allow the analysis of large quantities of samples as required to characterise a field or catchment-scale hydraulic properties, which are necessary for risk assessment (e.g., flood forecasting) and management (e.g., fertilization and pest control). The determination of pore size distributions by NMR relaxometry has its own drawbacks. First, the diffusion in induced magnetic field gradients can shorten transverse relaxation times (T 2 ). This must be checked and can be minimized by choosing sufficiently small echo times or by measuring longitudinal relaxation (T 1 ). Second, diffusion in internal gradients may affect different modes of a multimodal relaxation time distribution function in a different way. Third, one must be aware that the derivation of PSD is always a scaling procedure, which requires independent determination of the average specific surface area. The reason is that the relaxation times are controlled by the pore sizes and surface relaxivity, but the average S/V is controlled by the internal surface area of the porous media. Fourth, the assumption of the homogeneous distribution of pores and paramagnetic centres and the calculation of an average surface relaxivity parameter, especially for a multimodal relaxation time distribution function, leads to an overestimation of the large pores. All these issues, together with the necessary simplification of pore shape and geometry, can deviate the calculated pore size distribution from the real one (Stingaciu et al., 2010).
Nuclear magnetic resonance (NMR) relaxometry uses the random motion of molecules, whereas cryoporometry uses the melting-point depression of a confined liquid determining pore size distributions. It is suitable for measuring pore diameters in the range 2 nm-1 μm, depending on the absorbate. While NMR cryoporometry is a perturbative measurement, the results are independent of spin interactions at the pore surface and so can offer direct measurements of pore volume as a function of pore diameter (Mitchell et al., 2008). NMR cryoporometry (Strange et al., 1993) rely on the Gibbs-Thomson equation concerning the relationship between the characteristic pore length scale and the change in the freezing point of the liquid, or melting point of its solid crystal, due to confinement within the porous matrix. However, the freezing and melting behaviour of confined liquids/solids is often complex, with the thermodynamic properties of the confined material being modified from those of the bulk (Christenson, 2001). NMR cryoporometry offers the advantage of a more direct measure of the open pore volume.
The pore-size distribution obtained from cryoporometry is the derivative of the total liquid water content with respect to temperature. Because noise is always present in the NMR data, the derivative can have an irregular shape for pores having only a minor contribution to the total signal. However, no a priori shape of a distribution is imposed. On the other hand, the pore-size distribution obtained from relaxometry is the result of a numerical inverse Laplace transform. This is generally an ill-posed problem and hence some sort of regularization method is used. The numerical code that we use always yields a sum of log-normal Gaussian-shaped pore-size distributions. Therefore the shape of the relaxometry pore-size distribution does not accurately reflect the actual pore-size distribution. Only the position of the peaks and the total area of the distribution around a specific peak are relevant parameters in this pore-size distribution (Valckenborg et al., 2002).

Mercury intrusion porosimetry (MIP)
Mercury intrusion porosimetry is by far the most popular method for characterising porous materials with pores in the range of 500 μm down to 3 nm (Giesche, 2006;León y León, 1998;Rouquerol et al., 2012). Compared to alternative characterisation methods such as gas sorption, mercury porosimetry covers a much wider pore size range, while it is based on simpler physicochemical principles and it is much faster in operation. The most important limitations of the mercury intrusion when applied to extract pore size distribution of a porous material is that it is based on the assumption that the porous matrix can be represented by a bundle of cylindrical pores. In extracting pore size distribution during mercury intrusion, all pores are assumed to be equally accessible to the exterior mercury reservoir. This assumption can only be met if the pore structure is represented in the form of a bundle of capillaries or if pore connectivity is very high. In reality, however, pore network effects can be quite important resulting in the so-called pore shadowing or ink-bottle phenomenon. In this case, a large pore has smaller entrances connecting with the mercury reservoir and can only be filled at a pressure higher than that required by its actual dimensions. Capillary network models of varying complexity have been used over the past three decades and have been often quite successful in proving structural characteristics of the porous materials examined (Ioannidis and Chatzis, 1993a;Kikkinides and Politis, 2014). For obvious reasons this method could be applied to study only the open porosity features, since the closed pores are inaccessible to the mercury (Giesche, 2006).

Gas adsorption
Adsorption behaviour of porous materials is a function of their microstructural characteristics, such as the surface area, types of pores present in the material, topology of the porous network and the available pore volume. Characterisation of porous materials is, therefore, important in the development of adsorption applications. Mercury porosimetry and gas adsorption are determined from surface tension, capillary forces and pressure. The main difference is that with mercury porosimetry, large pores at the intrusion phase are determined first, while with gas adsorption, the smallest pores are measured first at the adsorption phase.
Physical gas adsorption has been used to study the pore characteristics of solid materials and the changes therein upon post-synthesis treatment. Gas adsorption is one of the most popular techniques used to characterise the pore space. It only allows determining the volume of open pores while closed porosity cannot be accessed. The advantage of this technique is that it allows assessing a wide range of pore sizes, covering essentially the completed micro-and meso-pore range in a timely and cost-effective fashion (Thommes, 2010).
Frequently used adsorptives are nitrogen (N 2 ), argon (A r ), and CO 2 , depending on the nature of the material (adsorbent) and the information required. Nitrogen at 77 K is considered to be a standard adsorptive for surface area and pore size analysis, but it is meanwhile generally accepted that nitrogen adsorption is not satisfactory with regard to a quantitative assessment of the micro-porosity, especially for micro-pores with widths smaller than 0.7 nm. Consequently, alternative probe molecules have been suggested, e.g., argon and carbon dioxide. For many micro porous systems, argon adsorption at 87.3 K appears to be very useful as N 2 adsorption in microspores occurs at lower p/p 0 values than Ar (Groen et al., 2003;Ravikovitch et al., 2000;Serrano et al., 2009). Despite this advantage, the low pressures induced restricted diffusion prevents argon molecules from entering the narrowest micropores, i.e., pores of widths b ca. 0.45 nm. Furthermore, Ar adsorption at 77 K shows limited application for meso-pore size determination, since the coolant temperature is below the bulk triple point. As a consequence, pore condensation vanishes in case the pore diameter exceeds approximately 12 nm (Thommes et al., 2002). When compared to nitrogen and carbon dioxide, it exhibits weaker attractive fluid-pore wall attractions for most adsorbents, which during adsorption does not give rise to specific interactions (like nitrogen and carbon dioxide because of their quadrupole moments) with most of surface functional groups and exposed ions. CO 2 is another usually preferred adsorptive, since these adsorption measurements are mostly performed at temperatures near ambient, which will enhance diffusion properties in the highly micro porous system compared to the low temperatures used in N 2 and Ar adsorption (Rouquerol et al., 2013). A drawback of CO 2 adsorption at ambient temperature is that in most commonly used equipments, which predominantly operate in the pressure range of vacuum to 1 bar, only a limited range of micropores can be measured, unless high-pressure CO 2 adsorption is used (Cazorla-Amorós et al., 1998;Ravikovitch et al., 2000). While CO 2 adsorption at 273 K is frequently used for the ultramicopore analysis of carbonaceous materials (Rodriguez-Reinoso, 2009), it is not a good choice for the pore size analysis of materials with polar sites, mainly because of the very specific interactions that CO 2 can have with functional groups on the surface.
After obtaining the physical adsorption isotherms, the interpretation of experimentally measured isotherms is required in order to arrive the specific characteristics of the material (Thommes, 2010). Classical theories include the Brunauer-Emmett-Teller (BET) method (Brunauer et al., 1938), Barrett-Joiner-Halenda (BJH) (Barrett et al., 1951), Broekhoff-de Boer (BdB) method (Broekhoff and De Boer, 1967), Horvath-Kawazoe (HK) model, (Horvath and Kawazoe, 1983), and the Saito-Foley (SF) model Foley, 1991, 1995). Typically, these theories involve certain assumptions, which are not universally applicable. For example, the notion of an adsorbed monolayer required by the BET theory to assess the surface area is clearly not applicable to materials with narrow pores and complex energy landscapes. The classical BJH, based on the Kelvin equation and corrected for multi-layer adsorption, is most widely used for calculations of the pore size distribution (PSD) over the mesopore and part of the macropore range. However, the BJH and BdB theories of extracting PSD, based on a description of a porous material as a collection of cylindrical or slit pores, are not appropriate for zeolites, metal organic frameworks, polymers etc. The conventional HK model for slit-shaped pores and SF model for cylindrical pore geometry are mainly applied for micro-pore size calculations (Groen et al., 2003). Such approaches allow for obtaining the pore size distribution in addition to the pore volume, but rely on similar macroscopic and thermodynamic assumptions concerning the nature of confined adsorbate. This leads to inaccurate determination of the pore size and volume. New models based on non-local density functional theory (NLDFT) and molecular simulations were developed (Lukens et al., 1999;Neimark and Ravikovitch, 2001;Neimark et al., 2000). It has been demonstrated that the application of these novel theoretical and molecular simulation based methods leads to: (i) a much more accurate pore size analysis (Neimark and Ravikovitch, 2001;Thommes et al., 2006), and (ii) allows performing pore size analysis over the complete micro/mesopore size range (Thommes et al., 2006). Currently there are methods for pore size quantification based on NLDFT and molecular simulations, which are commercially available and applicable to a range of adsorptive/adsorbent systems. They include hybrid methods that assume various pore geometries for the micro-and meso-pore size range, as it can be found for materials with hierarchical pore structures.
Previously published review on the use of gas adsorption for characterisation of porous materials hardly comments on these phenomena (Sing, 2001). One drawback is that the process could be very time consuming, but the determinable pore diameter is from 0.3 to 300 nm, a range not completely covered by mercury porosimetry.

Pore network model construction
The different experimental techniques for pore space characterisation present different limitations to the way representative pore networks can be constructed. For example, if a 3D image is available (experimental or synthetically generated), pore networks can be constructed directly from this image, assuming what is considered to be a pore and what to be a throat between pores. The need for up-scaling from the typically small imaged volume to larger domains leads to the need for construction of statistically representative networks. This requires analysis of the image to extract size distributions of pores and throats and their connectivity. However, if the experimental data comes from non-imaging techniques such as mercury intrusion porosimetry and gas adsorption where not all pore space characteristics are readily available, regular pore network construction approach is usually applied with assumed connectivity (Jivkov and Olele, 2012).
To a large extent the success of pore network models depends on the way they represent the real pore space in terms of its geometrical and topological characteristics for a given application. Previous works have clearly demonstrated the importance of the geometric properties of the porous media, such as the locations of pores and throats, the distributions of sizes and shapes of pores and throats Knackstedt et al., 1998;Oren, 1994). Most early works on pore-scale modelling assumed that the throats were cylinders with a circular cross-section or considered to be zero-volume connections between pore chambers. Pores were either not modelled explicitly at all (they simply represent throat junctions), or were spherical, cubical or cylindrical in shape. These methods have been applied to study the convection, transient and steady-state diffusion, permeability, etc. (Bryntesson, 2002;Jivkov and Xiong, 2014;Laudone et al., 2008;Meyers and Liapis, 1999;Meyers et al., 2001). Recent advances have allowed modelling a degree of irregularity in pore cross-sectional shape that was not available in earlier PNM (Al-Gharbi and . For some applications, such as multiphase flow, it may be necessary to pay specific attention to the local morphology of the throats (Fenwick and Blunt, 1998;Gao et al., 2012;Hui and Blunt, 2000;Man and Jing, 1999;Payatakes et al., 1973). For others, such as longer-scale permeability or diffusivity predictions, the throats could be considered as straight channels with locally averaged cross sections (De Josselin de Jong, 1958). In the latter case, the variability in throat morphologies becomes of secondary importance and the transport is controlled predominantly by the sizes and spatial positions of pores, and the connectivity of the pore set via throats with different permittivity. In addition, the models have to reflect basic topological properties, such as average pore coordination number and pore coordination spectrum. The effects of the average pore coordination on transport coefficients have been demonstrated in a number of works (Meyers and Liapis, 1999;Raoof and Hassanizadeh, 2010a). The pore spectrum represents the relative number of pores coordinated by different number of throats. Introducing the spectrum into the model poses a stronger constraint onto the construction than just utilising the average pore coordination number. This is due to the fact that totally different spectra could produce the same average coordination number and the effects of it have been demonstrated in a number of works (Jivkov et al., 2013;Meyers and Liapis, 1999).
Generally speaking there are three ways to construct a PNM representing a porous medium. The first method is to create a statistically equivalent network using distributions of basic morphologic parameters, while the second approach is to map a network structure directly onto a specific porous medium void space. The fundamental difference between the two methods is that the direct mapping provides a one-to-one spatial correspondence between the porous medium structure and the equivalent network structure, whereas the other type of network is equivalent only in a statistical sense to the modelled system. The last method is called the grain-based approach, which is based on the diagenesis of porous media.

Statistical reconstruction
2D pore space images are routinely available at high resolution. 3D images can be reconstructed using statistical methods with information obtained by analysing 2D thin sections. Methods based on a truncated Gaussian random field are often used in conjunction with the geometrical properties of the original pore space to reconstruct 3D images (Adler and Thovert, 1998). These geometrical properties include porosity, also called one-point correlation function, and two-point correlation function measuring the probability of finding two points separated by a certain distance within the same phase.
However, the one and two-point correlations are insufficient to adequately replicate the topology of the medium (Adler and Thovert, 1998;Ioannidis and Chatzis, 2000;Levitz, 1998;Roberts and Torquato, 1999;Yeong and Torquato, 1998a,b). Yeong and Torquato (1998a,b) used a combination of the two-point correlation function and the distribution of linear path, which is the probability of finding a line segment with certain length in the void space as a descriptor to characterise the pore geometry. In addition, Hilfer (Hilfer, 1991) introduced local porosity distribution and local percolation probability to improve the geometrical characterisation. Other descriptors, such as pore chord length (the length in the void between two solid voxels with a given direction) have proved useful in characterising the structure and generating 3D images (Levitz, 1998;Roberts and Torquato, 1999). The combination of oneand two-point correlation functions with these geometrical descriptors improves the reconstruction of connectedness and predictions of macroscopic properties such as permeability.
Despite this, these methods still fail to reproduce the long-range connectivity of the original pore space. On the other hand, Yeong and Torquato (1998a,b) developed a stochastic method based on simulated annealing, which was later extended by Manwart et al. (2000). Their approach is based on moving pore space voxels around to minimize the objective function and they obtained the correct porosity. This method should be able to match not only one-and two-point correlation functions but also multi-point correlation functions. Okabe and Blunt (2004) developed a multi-point statistical method aiming to reconstruct the 3D volume from thin section images. The approach preserves the typical void space patterns observed in 2D and consequently preserves the long-range connectivity. These statistical methods discussed above produce 3D representation form 2D images of the porous media with similar morphological statistics.

Grain-based model
Bryant and co-workers pioneered the use of geologically realistic networks (Bryant and Blunt, 1992;Bryant et al., 1993a,b). Their models are based on random close packing of equally-sized spheres. They represented diagenesis by swelling the spheres uniformly and allowing them to overlap. Compaction was modelled by moving the centres of the spheres closer together in the vertical direction, again allowing the spheres to overlap. Equivalent networks with a coordination number of four or less were then constructed. Single and multiphase flow was simulated through the pore space. They were able to predict the absolute and relative permeability, capillary pressure, electrical and elastic properties of water-wet sand packs, sphere packs and cemented quartz sandstone. They also reported a trend between permeability and porosity for Fontainebleu sandstone. This represented a major triumph in pore-scale modelling, since genuine predictions of transport and flow properties were made for the first time. They showed that spatial correlations in the pore size distribution were important for correct predictions: using the same pore size distribution, but assigning it at random to the throats in the network gave erroneous predictions of permeability (Bryant et al., 1993a). The main drawback with the work is its limited application-it could only be applied to media that is predominantly composed of spherical grains of the same size.
The next major advance came with the work of Oren, Bakke and coworkers at Statoil (Bakke and Oren, 1997;Lerdahl et al., 2000;Oren and Bakke, 2002;Oren et al., 1998). They developed a reconstruction method, where the packing of spheres of different size was simulated. The grain size distribution was derived directly from analysis of thin sections of the rock of interest. Compaction and diagenesis was modelled in a similar manner to Bryant et al. (1993a). Clays were also included in the model. Biswal et al. (1999) compared the pore space derived from this geological reconstruction of Fontainebleau sandstone with an image obtained from micro tomography. They also studied two stochastic models based on a correlation function representation. It was shown that the stochastic models differed strongly from the real sandstone in their connectivity properties. In contrast, the geological reconstruction gave a good representation of the connectivity of the rock and as a consequence could accurately predict transport properties (Oren and Bakke, 2002).
Oren and co-workers used the geological reconstructions to formulate topologically equivalent networks through which multiphase flow was simulated (Bakke and Oren, 1997;Oren et al., 1998). They predicted relative permeability for a variety of water-wet sandstones, showed promising results for a mixed-wet reservoir sample and matched three-phase water-wet data (Lerdahl et al., 2000;Oren and Bakke, 2002). The works of Pillotti, and Coehlo and co-workers have demonstrated a method to simulate the deposition of grains of non-spherical shapes, in particular, enabling the reconstruction technique to be applied more generally (Coelho et al., 1997;Pilotti, 2000).
There are three major concerns with the application of reconstruction methods to reservoir samples. Firstly, the reconstruction algorithm is based on explicit simulation of the geological processes by which the rock is formed. For many complex systems, involving microporosity and clays as well as a variety of different sedimentary processes, this may prove challenging. Furthermore, carbonate systems are not modelled at all. Statistical reconstruction methods are more general, since they require only a two-dimensional image of the system, and have been applied successfully to non-classic rocks. However, their application so far, has been mainly to predict single-phase flow simulated directly onto the pore-space reconstruction. The second problem, regardless of the approach, is that characterisation of the pore space requires detailed thin section analysis that might be unavailable or difficult to obtain. The final issue is that the appropriate characterisation of pore shape and wettability are not fully understood.

Direct mapping model
Direct mapping from a real sample will yield an irregular lattice (Piri and Blunt, 2005). Such irregular networks allow for validating physical assumptions for flow simulations by comparison with 4D imaging of mass transport. 4D (3D + time) imaging methodology enables visualization and quantitative assessment of dynamic pore scale processes in real time over variable experimental durations. More details can be found at http://www.solid-earth-discuss.net/se-2016-40. Based on 3D images, the following approaches have been used to construct irregular pore network models.

Medial axis algorithm
The medial axis based methods transform the pore space images into a medial axis that was the reduced representation of the pore space acting as a topological skeleton. The topological skeleton is constructed roughly along the middle of the pore channels either by a thinning algorithm (Baldwin et al., 1996;Liang et al., 1998) or a pore-space burning algorithm (Lindquist et al., 1996). The algorithms developed by Lindquist et al. (2000) have been primarily applied to consolidated systems but it is unknown whether they are applicable to unconsolidated porous media with higher porosity. Al-Raoush et al. (2003) developed algorithms which can be applied to unconsolidated porous media and demonstrated their application on a variety of unconsolidated porous media systems and the impacts of grain sizes, REV, and image resolution (Al-Raoush and Willson, 2005). The pore throats are decided by local minima along the branches whereas the pore bodies are determined at the nodes, which verify pore space partitioning (Jiang et al., 2007;Lindquist et al., 1996;Prodanović et al., 2006).
The medial axis mathematically preserves the topology of the pore space; it is difficult, however, to identify pores unambiguously. Furthermore, pores normally encompass more than one junction of the medial axis; therefore, various merging algorithms have to be developed to trim the skeleton and fuse the junctions together while avoiding unrealistically high coordination numbers (Jiang et al., 2007;Shin et al., 2005). The choice of threshold value for throat quality can be problematic without being examined by real flow simulations. In conclusion, medial axis algorithms readily capture the interconnectivity of the pore spaces but encounter the problem in identifying pores.

Maximal ball algorithm
The maximal ball algorithm (Al-Kharusi and Blunt, 2007;Silin and Patzek, 2006) finds the largest inscribed spheres centred on each voxel of the image that just touch the grain or the boundary. Then those spheres, which are included in other spheres, are removed, while the rest are called maximal balls. The largest maximal balls usually identify pores, which are connected by smaller balls between them (these balls are called throats). Maximal balls were first used by Silin and Patzek (2006) to find the dimensionless capillary pressure in drainage rather than to extract a pore network. The method was then adopted and developed by Al-Kharusi and Blunt (2007). A more comprehensive set of criteria was set to determine the maximal ball hierarchy including sphere clusters to handle equally sized balls. However, the algorithms of Al-Kharusi and Blunt use a tremendous amount of computer memory and time and consequently were limited to relatively small systems containing fewer than a thousand pores. Moreover, their method tended to form pores with very high coordination numbers. To overcome this issue, a two-step searching algorithm was developed by Dong and Blunt (2009) based on Al-Kharusi and Blunt's algorithm. The nearest solid was found to define a void bass instead of growing a ball layer by layer. This algorithm invented a clustering process to define pores and throats by affiliating the maximal balls into family trees according to their size and rank (Dong and Blunt, 2009). This method clearly identified the larger pores, but tended to indentify a cascade of smaller voids, which have sizes less than the image resolution.
One of the problems in network generation and modelling is the lack of a specific definition for what constitutes pores and pore throats in real materials with complicated pore geometries. Another problem is how the definition of pore morphology affects network topology and discretisation. Consequently, multiple network structures can be created for the same material or 3D data set and there is no straightforward way to define whether one network is more physically representative than another. Vogel (2000) has demonstrated that different latticebased topologies have similar water retention curves. Inversely, Ams et al. (2004) have demonstrated that different network topology induce different relative permeability.
The effects of pore density (number of pores per unit volume) on transport in porous media have been studied by Bhattad et al. (2011). The large variations in pore density are the result of two factors: different methods for seeding the search for pores (i.e. defined as maximal inscribed spheres) and different settings in the merging criteria (that merges overlapping inscribed spheres to created single pores). Generally, the fast algorithms create low-density networks with little or no overlap in inscribed pores and therefore, can also generate longer pore-throats. More complete pore searches are slow and result in higher pore densities and shorter pore-throats. The pore density affects secondary parameters, such as pore-size distribution and pore interconnectivity. However, the relationship is not always proportional or obvious. The most important results have shown that single-phase permeability is relative insensitive to pore density, while capillary pressure curves for quasi-static drainage are moderately sensitive to network structure (Bhattad et al., 2011).

Regular network model
Irregular PNM are sample-specific and potentially not statistically representative. Direct construction of irregular network is impossible (especially for micro-and meso-porous materials) in the absence of data for throat locations and sizes, which cannot be observed by current imaging techniques. A regular PNM constructed within a larger volume allows for capturing statistics from a number of imaged samples, i.e. improved statistical correspondence, and for calculating transport at distances closer to the engineering length scale.
Early three-dimensional network models are usually based on cubic lattice with a constant coordination number of six, Fig. 1(a) (Ioannidis and Chatzis, 1993b;Reeves and Celia, 1996). However, in actual porous media, the connectivity is distributed and can have coordination number greater than six (Dong and Blunt, 2009;Jivkov et al., 2013;Kwiecien et al., 1990). Therefore models with 26 pore coordination number have been proposed, Fig. 1(b) (Raoof and Hassanizadeh, 2010a). This does not appear to be physically realistic, however, because large numbers of throats intersect at points that are not pores. A pore network construction based on the Kelvin solid has been recently proposed and applied to permeability and diffusion analyses of porous media, Fig. 1(c) (Jivkov et al., 2013;Xiong et al., 2014). This is based on truncated octahedral cells providing variable pore coordination dependent on the allocation of pores and throats in the cell, with maximum coordination of 14. Another application of a non-cubic lattice support for pore-network construction uses rhombic dodecahedron as a unit cell, Fig. 1(d) (Vogel and Roth, 2001). This offers a maximum pore coordination of 12 with throats of equal lengths, i.e. sufficiently large physically admissible coordination. However, this cell is less representative for the volume around a pore than the truncated octahedral cell.
Historically, pores were related to the sites while pore throats were related to the bonds of regular PNM Meyers et al., 2001;Wilkinson and Willemsen, 1983). If such correspondence is to be statistically representative of the material modelled, sufficiently rich experimental information is required. This information includes shape and size distribution of pores and throats, as well as the pore coordination spectrum, i.e. percentages of pores coordinated by different numbers of throats (Gao et al., 2012;Jivkov et al., 2013). These can be obtained in structures with distinguishable pores and pore throats.
For porous media with macro-porosity, typically above 100 nm such as sandstone and limestone, the throat and pore sizes distribution can be resolved as well as the coordination spectra and the average coordination numbers can be calculated (Al-Raoush and Willson, 2005;Dong and Blunt, 2009). In such cases, topologically representative networks can be constructed by a system of pores and throats connecting all neighbouring pores, and subsequent elimination of throats to achieve the topological constraints. This has been illustrated for average coordination number in simple bases, such as cubic lattice (Raoof and Hassanizadeh, 2010a), as well as for full coordination spectra in a bi-regular lattice (Jivkov et al., 2013).
However, for micro-and meso-porosity, e.g. below 100 nm, the pore size distribution can usually be determined, but the resolution of the experimental techniques is not sufficient to segment all the throats and calculate their sizes. For such cases, different approaches to pore network construction are required. One possibility is to assume full connectivity between neighbouring pores in a selected lattice, but control the diffusivity of the throats through the sizes of the connected pores and the size of the solute molecules . This construction leads to a model length scale dictated by the prescribed total porosity and the assumption that pores are located at each lattice site. The approach was shown to provide insights into the effects of the structure on diffusivity and the results correlated well with experimentally measured diffusion coefficients . This approach, however, is not sufficient to explain the variability in reported experimental measurements of diffusion coefficients.
Clearly, the microstructure information used for the network construction is derived from samples, which are different from the ones used for measurements of diffusion coefficients. It is plausible, however, to assume that the statistical part of this information, namely the distribution of pore sizes, is not vastly different between samples of the same medium. This implies that there is an effect of an additional geometrical parameter on the measured coefficients. It could be possible that a network length scale that provides a constraint to the model rather than being determined from total porosity and location of pores in all lattice sites .
The two approaches that tackle incomplete pore space information, namely predefined connectivity for calculable length scale , and undetermined length scale for realistic connectivity , suffer from the lack of an additional constraint. This cannot be found from the pore space information and requires the consideration of the solid phase structure, e.g. the shape and size distribution of mineral grains. A methodology for incorporating the structure of the solid phase was then proposed which improves substantially the realism of the constructed pore network, both in terms of geometry and topology . The proposed methodology is conceptually different from previous works and it has the added benefit of constructing a PNM which can be paired directly to lattice models of the solid phase, developed for analysis of damage evolution via microcracking . The approach is conceptually similar to the dual models of pore spaces introduced by Hilpert (2007, 2008). This facilitates mechanistic investigations of combined mechanical-thermal-chemical-biological effects on diffusivity.  (Ioannidis and Chatzis, 1993b;Reeves and Celia, 1996); (b) cubic lattice with 26 pore coordination number (Raoof and Hassanizadeh, 2010a); (c) truncated octahedron (Jivkov et al., 2013;Xiong et al., 2014); (d) rhombic dodecahedron (Vogel and Roth, 2001).

Two-scale pore network models
Despite the advances of pore scale modelling, the simulation of transport properties in heterogeneous rocks still remains an open issue due to the very broad pore size distributions in some porous media. Such examples could be of considerable scientific and economic importance including many types of carbonates, clay-rich sandstones and clays. For example, carbonate reservoirs or tight gas sandstones are very important for hydrocarbon extraction and CO 2 sequestration. The failure of simulating most classical empirical relations (e.g. the Brooks-Corey parameterization of relative permeability and Archie's law for electrical behaviour) in these rocks due to the interaction of micro-and macro-porosity in turn has promoted the development of modelling methods where multiple pore scales are coupled (Mehmani and Prodanović, 2014). A short overview of existing two-scale pore network models is presented here and illustrated in Fig. 2. Jiang et al. (2013) developed a methodology to integrate networks extracted from images at distinct length scales. The pore network model was generated at each scale and then was integrated into a single two-scale network by characterising the cross-scale connection structure between the two networks. The shortcoming of this method is that it is computationally costly due to the number of network elements (Bultreys et al., 2015). Recognizing the computational problems when single micro-pores are taken into account, Mehmani and Prodanović (2014) proposed a two-scale pore network by the use of packing algorithms. The macro-network is constructed by Delaunay tessellation of the grain centres. Micro-porous networks are generated by downscaling existing networks extracted from macro-pores. This approach was capable of investigating fundamental two-phase flow properties of multi-scale porous media and micro-porosity was able to act in series (intergranular or pore-filling micro-porosity) and in parallel (intragranular or dissolution micro-porosity) to the macro-pores. However, in the construction process, distorted pores were produced when many small grains touched a large grain. A parameter, the ratio of the macro to micro length scale, needs to be determined for micro-porous regions from the measurements, which may be difficult to perform. Bultreys et al. (2015) developed a workflow to integrate networks of macro-pores and micro-porous regions extracted from micro-CT images. This methodology allowed micro-porosity to act both in parallel and in series with the macro-pore network. However, a representative network for the micro-porosity is necessary. In addition, the pore networks from Jiang et al. (2013) and Bultreys et al. (2015) were based on experimental data from micro-CT images which did not take into account the micro-pores that cannot be resolved by micro-CT. As the  (Jiang et al., 2013;Prodanović et al., 2015). (b) a lattice-based model with upscaled microporosity properties (Bekri et al., 2005). (c) a model which takes upscaled microporosity in parallel with macro-throats (dark grey) into account (Bauer et al., 2012). (d) a model with connectivity added by the (upscaled) microporosity (Bultreys et al., 2015). truncated cone shape is used to connect two neighbouring macro-pores, the tortuosity of the connection and geometric details about the bulk of the micro-porous cluster are neglected, which can lead to erroneous local conductivities.

Adsorption
PNM has been widely used to study adsorption in porous media which finds its relevance in many areas of science and engineering including radioactive waste disposal (Xiong et al., 2016;Xiong et al., 2014) and chromatographic separation processes (Câmara and Silva Neto, 2009;Meyers and Liapis, 1999).
There are two principal ways to simulate adsorption in PNM. The first one entails decoupling of adsorption and transport into two separate processes by the use of parameters associated with each process calculated separately. Firstly, the convection-diffusion problem is solved to obtain the concentration and flow field. This is governed by the following equations: where Eq. (1) is the Hagen-Poiseuille law, appropriate for describing flow in pores (Jacob, 1972) and valid for laminar flow; and Eq. (2) is the classical convection-diffusion law. The equations are written for transport through a throat connecting pores i and j, with length l ij [L], and cross section area A ij [L 2 ]. In Eq. (1) q ij [L 3 T −1 ] is the total volumetric flow rate through the throat, p i and p j [ML −1 T −2 ] are pressures at pores i and j, respectively, μ [ML −1 T −1 ] is the dynamic viscosity, and G ij is a cross-section shape coefficient: for circular, equilateral triangle and square tube, G ij is 0.5, 0.6 and 0.5623 respectively (Patzek and Silin, 2001). In Eq. (2) V ij [L 3 ] is the throat volume, c i and c j [ML −3 ] are concentrations at pores i and j, respectively, and D ij [L 2 T −1 ]is the throat diffusivity.
Then, a sorption isotherm, which is assumed to belong to some class of parameterized functions (Langmuir, Freundlich, Henry, etc.) with few tuning parameters, is used to describe the adsorption of species onto the walls. A coarse grained mathematical formula containing the sorption isotherm is used to reflect the effects of sorption on the change of pore sizes. For example, the adsorbate obstruction can be simulated by reducing the throats radii as a function of the concentration of the moving species (Meyers and Liapis, 1999;Xiong et al., 2016Xiong et al., , 2014 where t is the thickness of the adsorbed layer, R is the throat radius, r M is the radius of the diffusing species, and θ is the adsorption isotherm (e.g. θ ij = k d,ij c ij , θ ij = ac ij b ). Thus, the radius of a throat after adsorption at given concentration, C A , and hence radius after adsorption, becomes R ⁎ = R − t. Thus, the sorption is accounted for by modifying the pore geometries according to the changed pore sizes. This procedure is iterated until reaching a certain criteria of interest. Using this method, various types of adsorption reactions can be studied, including linear equilibrium, nonlinear equilibrium (Meyers and Liapis, 1999; and heterogeneous adsorption in which the adsorption parameters are spatially varying (Xiong et al., 2016).
Sorption isotherms in the published literature, however, are determined from bulk measurements and may not be appropriate for occluded regions within the pore space. Adsorption takes place at the watersolid interface and is thus controlled by solute concentration in the fluid next to the solid grain in individual pores. The concentration near the solid surface is different from the average concentration within the pore due to the development of concentration gradients within the pore space. Simulations with pore network models are usually based on the average concentration in individual pore throats. An alternative way is to simulate the adsorption in each pore element and then scale it up to the whole pore network. This has been demonstrated by Raoof and Hassanizadeh (2010b). They have performed pore-scale simulations for a wide range of local-scale distribution coefficients and Peclet numbers. These have provided relationships between the upscaled parameters and underlying pore-scale parameters. Such relations are useful for upscaling by means of a pore-network model. Raoof and Hassanizadeh have shown also that even if there were equilibrium adsorption at the pore wall (i.e. at the grain surface), one might need to employ a kinetic description at larger scales. They have also shown this kinetic behaviour through a volume averaging method, yielding very similar results for the upscaled kinetic parameters. From upscaling, they have found a negligible dependency of effective adsorption rate coefficients on pore water velocity. This finding is in agreement with Zhang et al. (2008) who has observed that upscaled sorption parameters are independent of pore-water velocity. The findings, however, contradict some experimental results, where kinetic adsorption-desorption coefficients have been shown to depended on velocity (Akratanakul et al., 1983;Brusseau, 1992;Maraqa et al., 1999;Pang et al., 2002).
The second principal method for simulating adsorption in PNM is by solving the advection-diffusion-sorption equation directly. The adsorption reactions include linear equilibrium (Raoof et al., 2013) and nonlinear equilibrium (Acharya et al., 2005;. The equation can be expressed as where V ij u ij c i −c j l ij is the advection term which represents the effects of velocity u ij on particle transport, F ij ¼ A ij D ij d 2 c ij dx 2 is the diffusive term, s ij = s ij (c ij ) [ML − 2 ] is mass adsorbed per unit area of throat wall, S ij [L 2 ] is the throat surface area.

Dissolution and precipitation
Modelling of dissolution/precipitation is of high interest to applications such as CO 2 geological storage, enhanced oil recovery and chemical degradation of the cement. For example, CO 2 storage is one possible method to reduce atmospheric emissions of CO 2 , and hence mitigates climate change. Some CO 2 may dissolve in the resident brine and generates a weak acid after injection into the low-permeability caprock. The acidic CO 2 -rich brine reacts with the host rock to produce solid carbonate which leads to precipitation process (Blunt et al., 2013). In enhanced oil recovery, strong acid has been injected near the well to dissolve the host rock and increase the reservoir permeability (Algive et al., 2010). Accurately computing dissolution and precipitation processes involve the development of advection, diffusion and reaction equations. In addition, for better estimation of the abovementioned process the following also need to be accounted for: the heterogeneities in the porous media (e.g. porosity, mineralogical and chemical compositions), geometrical changes (i.e. pore and throat size changes), topological changes (i.e. increase/decrease pore to pore connection) and the chemical gradient filed inside pores. This is due to the fact that mineral dissolution and precipitation reactions alter the geometry and topology of the porous media (Crandell et al., 2012). Li et al. (2006) used a pore-scale network model to simulate reaction of kaolinite and anorthite during carbon sequestration. The processes of diffusion, advection, aqueous reactions, dissolution and precipitation were simulated in individual pores. The information was then summarised and averaged over the network to obtain descriptions of processes at the continuum scale. As previously observed, they showed qualitatively that in situ reaction rates were significantly different from values found in batch lab reactors. Their network models, however, were 3D regular lattices and not physically representative of the real media, which made quantitative analysis challenging. Li et al. (2007) further investigated the influence of the spatial distribution of mineral species using hypothetical networks, while Kim et al. (2011) used network models mapped from X-ray computed tomography (XCMT) of real sandstones. They were able to determine the initial concentrations of minerals in the porous medium from the XCMT images, which allowed for more predictive modelling. Kim et al. (2011) studied the relationship of reaction rates with kinetic rate and bulk flow rate based on previous work (Li et al., 2007(Li et al., , 2006. The reaction rates were found dependent on flow rates with an approximately power law relationship. They noted that their models were limited by the size of the network model. Notably, the effects of decreased throat conductivity as a result of precipitation were not accounted for in their work. Thus the physical change in porosity is lacking in the network flow model due to surface minerals dissolving or precipitating. These models, therefore, do not capture the important topological changes that reactions produce and the impact that such changes may induce on upscaled reaction rates. Algorithms for incorporating physical change, due to dissolution and precipitation, need to be included in the models. Based on these limitations, Kim and Lindquist (2013) further extended the work to consider geometric changes (e.g. pore volume and reactive surface area) and topological changes, which opened new connections between porosity features. When geometrical changes are included, the steady state, core-scale reaction rates could not be directly used as input into larger-scale simulations, as reactions are time dependent. The Nernst-Planck (N-P) term in the governing transport equations, which ensures the electrical neutrality when species of different charges diffuse at different rates, was also taken into account. The agreement between simulations with and without the N-P term improved as the flow rate decreased. In these studies, reactions far from equilibrium, such as anorthite reaction, and reactions close to equilibrium conditions, such as the kaolinite reaction, were considered. Algive et al. (2010) used a pore-scale model to simulate reactive transport which caused precipitation in pores. The reactive transport and structural modifications were treated separately and updated step by step. The reactive transport was solved in each element of the pore network. The mean concentration in each pore was used to determine the effective parameters of the macroscopic reactive transport equation. The different pore geometries were also considered for reactive transport by using random walk techniques. Then the reactive transport equation was solved analytically between the nodes of the regular lattice network for the asymptotic regime and computed the node concentration by a matrix inversion. Moreover, the system was based on the assumption that it had nearly reached chemical equilibrium and mainly focused on regular or periodic networks in one or two dimensions or systems, where pore-to-pore heterogeneity was not taken into account. Algive et al. (2012) developed a methodology to use a reactive pore network model to extract upscaling factors to tie the pore-scale effects of reactive transport to the core-scale values of permeability and porosity. Their work simplified the geochemistry and transport by incorporating most geochemical dynamics and transport into dimensionless variables and in doing so prevented a full description of the different physical and chemical parameters. Mehmani et al. (2012) developed a novel approach that coupled several pore-scale models using mortar coupling domain decomposition to study the evolution of precipitation-induced cementation of calcite. They were able to study large changes in permeability and porosity by coupling of 64 porescale models (1 mm × 1 mm × 1 mm each) but with a limited description of the chemistry represented by two main parameters, the Damkohler number (D a ) and an "alpha" parameter which described the deviation from equilibrium of the precipitation reaction. Nogues et al. (2013) simulated carbonic acid-driven reactions in a 3-D network of pores to predict the changes in permeability and porosity at continuum scale. The pore network structure was based on a statistical characterisation of a synthetic microcomputed-tomography image of an oolithic dolostone. Only physical heterogeneities were considered in comparison to the work done by Kim et al. (Kim and Lindquist, 2013). To account for the changes caused by dissolution and precipitation, a mathematical construct was developed to modify the pore-pore conductivities by relating these changes to changes in pore volumes. It is important to note that while the surface area fractions assigned to different minerals do change, the model does not change the total surface area. Because the reactions are modelled at equilibrium, the reaction rate is independent of surface area. Thus the simulation results are negligibly affected in this situation. This model analysed the effects of different chemical, flow and mixing conditions on permeability. Raoof et al. (2013Raoof et al. ( , 2012) employed a coupled Complex Pore Network Model (CPNM) and Biogeochemical Reaction Network Simulator (BRNS) to solve reactive transport by pore networks. CPNM reproduced the microstructure of real porous media and calculated solute concentration values for both pore bodies and pore throats within the pore network. BRNS (Regnier et al., 2002) was responsible for handling a comprehensive suite of multi-component complexion, mineral precipitation and dissolution reactions, as well as reaction networks characterised by multiple kinetic pathways. The evolution of porosity and permeability was discussed in this work as well. In BRNS the processes of dissolution/precipitation were simulated by changes of radius of pores and throats, while the effects of opening/clogging pores/throats were not considered.
All modelling approaches mentioned above involve the assumption that mass transfer within the pore is not limited, and that the regular geometric solid-fluid interface and the aqueous concentrations are spatially uniform within a pore. Based on that, the pore network models do not resolve chemical gradients within pores and also for each single pore the average values for physical/chemical properties were used. The accuracy of this assumption in computing reactive flow requires further investigation.

Biomass growth
Biomass growth within porous media is of particular interest in a variety of industrial scenarios such as microbial enhanced oil recovery (Ezeuko et al., 2011), bioremediation of contaminated soil and water (Cunningham et al., 1991) and nuclear waste disposal. As the mass of the microorganisms grows, hydrodynamic properties of porous media change (Rolland du Roscoat et al., 2014;Stewart and Scott Fogler, 2002;Taylor and Jaffé, 1990). The relevance of PNM to biofilm growth modelling comes from the fact that the biomass is not homogeneously dispersed within the porous media and hence continuum models fail to adequately represent this pore-scale phenomena and its effect on the macroscopic properties of the media.
The morphology of biomass growth in porous media depends strongly on the structure of the porous medium, bacteria species and the prevailing hydrodynamic and nutritional conditions. Wimpenny et al. presented a good review of the factors affecting biofilm formation (Wimpenny et al., 2000). Three principal morphologies have been reported in the literature: flat and uniformly thick biofilms, mushroomshaped structures and streamers (Stoodley et al., 1999). Only the first two, however, have significant relevance to porous media since streamers generally occur at turbulent flow conditions. The most commonly modelled growth mode is a continuous layer on the surface of the soil grains (Suchomel et al., 1998). Sometimes, however, the biomass does not cover the entire surface of the soil grains, and thus microcolonies, i.e., a patchy biofilm, develop (Molz et al., 1986). In other cases, biomass grows in pores in the form of aggregates (Vandevivere and Baveye, 1992). When fungal growth occurs, mycelia can develop throughout the entire pore volume and envelop several soil grains (Dupin and McCarty, 2000). These forms of biomass growth can occur simultaneously (Dupin and McCarty, 2000). Biofilm morphology is important because it results in different clogging mechanisms. Rittmann (1993) noted that the difference between flat biofilms, and aggregate growth is crucial to model permeability reduction. It takes less biomass to plug a pore if it forms aggregates as opposed to continuous biofilms. Therefore, mushroom-shaped colonies have the potential to generate larger permeability reductions in comparison to the uniform biofilm morphology (Thullner, 2010).
The development of robust methods to engineer biofilms in porous media requires predictive mathematical models capable of determining the evolution of biofilms under different flow conditions. However, reported studies of biofilm evolution in micro models have revealed that the structural complexity of biofilms entrapped in a complex pore space makes accurate modelling difficult (Dupin and McCarty, 2000;Kim and Fogler, 2000). Several pore level models for biofilm growth under flow conditions have been developed by adopting different approaches. They include lattice Boltzmann-based models (Pintelon et al., 2009(Pintelon et al., , 2012, bacterial cell level individual-based models (Kreft et al., 2001;Picioreanu et al., 2004), and Pore Network Models (PNMs) (Kim and Fogler, 2000;Stewart and Kim, 2004;Suchomel et al., 1998;Thullner and Baveye, 2008). Ezeuko et al. modelled biofilm evolution in porous media using a new PNM that includes hydrodynamics and nutrient transport based on coupling of advective transport with Fickian diffusion, and a reaction term to account for nutrient consumption. Specifically, the PNM is used to examine the influence of biofilm formation on the hydraulic properties of porous media. In PNM models, pore geometries can range from a simple circular cylinder through to complicated voxel-by-voxel pore reconstructions from thin sections or micro-tomographical data. The PNM approach facilitates simulations of important processes by adopting simplified but realistic pore geometries while allowing implementations of other pore level physics and kinetics relevant to the process of interest (Ezeuko et al., 2011). Or and Tuller (2000) used the classic capillary tube model to represent the soil pores (i.e., the pores are simulated as a bundle of capillary tubes of different diameters) and studied flow in unsaturated fractured porous media. The analysis showed that the biofilm spatial distribution within the pore space has a significant effect on the hydraulic properties of the soil. Although the capillary model is widely used and is easy to obtain, it has a few shortcomings, primarily the over simplistic assumption regarding the binary nature of the pores (i.e., each pore is either completely water-filled or completely empty under a prescribed matric head). Moreover, the lack of pores connectivity and the unrealistic cylindrical geometry of the capillaries further reduce the model efficiency in determining soil hydraulic properties (Or and Tuller, 2000). Or and Tuller (2000) suggested to model the soil-pores by using channels having a triangular cross-section. The channels can be either water-filled (for saturated flow) or partially-filled (for unsaturated flow), where the water-filled fraction is related to the matric head through the curvature of the water surface. The lack of pore connectivity can be especially significant when dealing with biofilms. Thullner et al. (2004Thullner et al. ( , 2002 showed that under saturated conditions, modelling the pores as one-dimensional capillary tubes could not reproduce connectivity-related effects. These effects, such as flow bypasses generated from preferential biofilm growth, were demonstrated in experiments and were successfully reproduced using pore-network models. Rosenzweig et al. extended Or and Tuller (2000) concept and propose a framework in which the soil pores are simulated as a network of channels having a triangular cross-section. The representation of the soil pores by a network of triangular channels is superior to the capillary model in two aspects. It eliminates the binary (full/empty) behaviour and takes into account inter-pore connectivity. Furthermore, it provides basis on simulating the combined water flow, solute transport and biofilm growth under variably saturated conditions (Rosenzweig et al., 2013).
The coupled water flow, substrate transport, and biofilm growth under saturated conditions has been studied by several researchers by using pore network models (Dupin et al., 2001b;Kim and Fogler, 2000;Suchomel et al., 1998;Thullner and Baveye, 2008;Thullner et al., 2002). To study the effects of biofilm on porous media, the pore network models were built either of lattices of cylindrical tubes (Kim and Fogler, 2000;Suchomel et al., 1998;Thullner and Baveye, 2008), or two-dimensional channels (Dupin et al., 2001b). The biofilm was modelled as either a continuous layer coating pore walls (Suchomel et al., 1998;Thullner and Baveye, 2008;Thullner et al., 2002) or as isolated microcolonies (Dupin et al., 2001b;Thullner et al., 2002). As biomass growth depends on the environmental conditions and the type of bacteria, it is problematic to represent this accumulation appropriately in terms of its effect on porosity and permeability reductions. No universally applicable model has been developed for saturated porous media systems. Moreover, no attempts have been made to develop a model that is applicable to both saturated and unsaturated systems. This situation could be particularly relevant to conditions that trigger gas generation as a result of bacterial metabolic processes. Gas bubble formation and entrapment could also severely block the pore space. It is usually characteristic of denitrifying and methanogenic conditions where nitrogen (Soares et al., 1991) and methane gas (Beckwith and Baird, 2001) production is observed but not exclusively as CO 2 and other gases are produced as part of the metabolic processes of different species. Unsaturated conditions are crucial for understanding the flow and transport in environments such as groundwater recharge basins, bioremediation sites, and effluent irrigated fields. However, very few pore network models have been used to study the biofilm morphology under unsaturated conditions (Rosenzweig et al., 2013). Rosenzweig et al. proposed a pore network model of triangular channels to simulate variably saturated flows in biofilm-affected porous media. The effects of biofilm on the hydraulic properties of the network are demonstrated by examining three scenarios for the biofilm pore-scale morphology including plug (or microcolony) morphology and a uniform biofilm layer (Rosenzweig et al., 2013). The simulations demonstrated that the effect of biofilms on the hydraulic properties of the network is a complicated and nonlinear function that depends not only on the biofilm scenario but also on the saturation. A significant assumption in most existing biofilm simulation models is that biofilm-occupied regions are impermeable and that nutrient provision within the biofilm is purely as a result of diffusion. However, experimental studies have demonstrated that biofilm morphology is often extremely heterogeneous and can contain voids (e.g., (Flemming et al., 2000)). As a result, biofilms are not fully impermeable and can contain a significant amount of both static and dynamic water (Flemming et al., 2000). In the absence of biofilm detachment, the effect of permeability is only relevant when biofilm growth is nutrient mass transfer limited. Such conditions are prominent in subsurface conditions (Kim and Fogler, 2000) and are very relevant to bio-barrier scenarios where thick biofilms are desirable for permeability reduction. Thullner et al. studied the effects of permeability of bioflims on the hydraulic conductivity of porous media by pore network models (Thullner and Baveye, 2008). The computational pore network model consisted of cylindrical pores which were arranged in an orthogonal array in two or three dimensions. The nodes connecting adjacent pores are assumed to be volumeless. The biofilm is viewed as homogeneous with regard to its hydraulic and biological characteristics in the model. In each pore, microbial activity is assumed to take place exclusively within the biofilm. This model was able to simulate reductions of the hydraulic conductivity of the pore networks by two to three orders of magnitudes. These reductions are in good agreement with laboratory experiments with san columns or field situations (Thullner and Baveye, 2008). The amount of biomass in any given pore is controlled by biomass growth, decay processes inside the biofilm, and detachment of biomass due to shear forces. However, few pore network models simulate all the biomass growth or surface attachment mechanisms.

Discussion and conclusions
The paper reviewed three aspects of pore network modelling: main experimental techniques for characterisation of pore spaces; main methods for construction of pore networks reflecting available pore space characteristics; and application of pore network models to analysis of phenomena, which have not been covered in previous reviewssorption, dissolution and biofilm growthbut have increasing technological and scientific importance. The works reviewed demonstrate that pore network modelling is an efficient and useful approach to analysis of reactive transport at the meso-scale, less computationally demanding than direct methods and able to incorporate heterogeneity within larger rock volumes.
In many practical cases, the main source of uncertainty in pore space characterisation is the finite resolution of the imaging techniques. These contain a hierarchy of pores from nanometres to micrometres. The uncertainties could arise from hardware limitations or from compromising resolution to capture heterogeneities. The errors may arise from the relatively sparse numerical grid of the imaging techniques, whereby details of the structure, e.g. very narrow passages in the porous space and surface roughness of fibres, may not be reproduced by the images. Furthermore, the results for each individual sample may be affected by the methods used to minimize or remove noise and other artefacts present in the original reconstructed images. These are the main factors that currently limit the applicability of the imaging techniques especially for fine-structured materials. The limitation is particularly relevant to materials with tight pore spaces, such as clays and cement-based materials. For example in concrete, the so-called air porosity can be resolved by existing methods, but appears as disconnected set of pores. The transport pathways are formed by the gel porosity of the cement, as well as along the so-called interfacial transitions zones between cement and aggregates, the connectivity of which is not resolvable with the current techniques. Similar is the situation with compacted clays. In such cases, the connectivity of the transport pathways needs to be assumed and tested against macroscopic transport experiments until suitable connectivity is found (Xiong et al., , 2016. In addition to numerical uncertainties arising from discretisation and image quality discussed above, this deviation includes statistical variation arising from relatively small size and low number of samples used in computation. In principle, the accuracy can be improved by using higher resolution images and better statistics. The former condition may become feasible with the present rapid development of imaging techniques. The latter can be obtained by larger size or larger number of samples. Clearly, these improvements come with the cost of higher computational effort, which, however, does not appear critical from the point of view of applicability of the method. Another approach for overcoming the above mentioned problems is to use a combination of different techniques such as FIB/SEM, nitrogen adsorption and FIB (Keller et al., 2013).
It must be noted that pore network modelling is conceptually scale indifferent, i.e. the approach can be applied to any length interval, where the structure of the pore space has been experimentally observed and analysed. For example if a particular experimental technique allows for characterising pore features of sizes between 0.1 μm and 10 μm, then the corresponding PNM will take these into account. The elements of PNM, sites and bonds, are abstract and can be related to the measurable features in different ways depending on the available information. The individual pore networks can be extracted at each length scale and then integrated into a single multi-scale network by characterising the cross-scale connection structure between these networks (Bultreys et al., 2015;Jiang et al., 2013).
Pore network models have been proven to be an effective research tool to upscale reactive transport processes from the local scale to the continuum scale. The microscopic dynamics can be understood best at the scale of individual pores. It is challenging to directly correlate results obtained from pore scale simulations to improve quantitative predictions based on large-scale simulations. However, the results obtained from pore scale investigations contribute to the understanding of large-scale natural processes and improve large scale geotechnical applications. It is well known that some large scale reactive process may be difficult to predict or explain unless there is a well defined experiment that has been used to study it precisely. Pore network models allow for bridging the gap between laboratory simulations and larger scale transport. Although pore scale is difficult to evaluate with traditional experimental studies, detailed investigations allow for the identification of important pore-structure details that are essential and control the reactive transport processes at the engineering/continuum scale.
It should be mentioned that, in spite of the attractiveness and successful applications of pore network models, as any other numerical method they have certain limitations. First of all, successful simulation of transport requires adequate representation of the real porous media. Pore network models usually simplify the geometry of the pore space and sometimes cannot capture all characteristics due to lack of resolution in most simulations. Whether the pore networks could reasonably predict single and multi phase properties depends on whether all relevant details are identified and included. In many systems, the very small pores below the image resolution make a little contribution to the overall behaviour. In these cases, transport properties are well predicted when the principal connected pore space is modelled, while the smaller pores could be neglected. Otherwise, assumptions are required. Secondly, the detailed microscopic features, such as the exact molecular structures and intermolecular interactions cannot be represented explicitly in pore network models. Therefore, in such situations, many approaches have been developed to study upscaling reactive transport in porous media (Algive et al., 2010;Li et al., 2006;Varloteaux et al., 2013). Reactive transport processes are simulated at the pore scale, which allows for incorporation of physical and chemical effects' rate variations from pore to pore. The reactive transport is then upscaled from the pore scale to a heterogeneous porous medium with the aid of a sufficiently large three-dimensional pore network model. The PNM is used to determine the phenomenological transport coefficients and the porosity/permeability relationship, which can be used as inputs in reservoir scale simulation (Algive et al., 2010;Li et al., 2006;Varloteaux et al., 2013). The strength of the pore network models is that they can be developed elegantly to study the effect of pore space changing mechanisms. In most simulations of pore-scale mass transport the pore geometry is assumed to be constant, apart from the effects of reaction. However, in many important geo-systems the confining geometry changes slowly as a result of processes such as pressure solution (Tada and Siever, 1989) or rapidly because of fracturing and changes in pore pressure or matrix stress. It is important to include the coupling between diffusion and geo-mechanical processes for a number of important applications including oil recovery and carbon dioxide sequestration.
Modelling and simulation of mineral precipitation and dissolution on the pore-scale have not been subjects to extensive research in the geo-science community. In particular, most simulations have employed very simple models for dissolution and/or precipitation chemistry, and there has been little, if any, work on chemical processes coupled with multiphase diffusion. It is believed that this situation is likely to change rapidly. There is strong motivation for such simulations for both scientific discovery and practical applications. Furthermore, the simulated results are usually compared with experimental data, which will not give the same results realistically. For example, in the reactive diffusion, a small change in the pore space geometry can lead to a large change in the fluid-solid interface. Similarly, micro-scale roughness and trace impurities can also have a strong impact on the outcome of experiments. Consequently, experimental validation of computer models for reactive transport in porous media must be based primarily on the comparison of statistical measures such as correlation functions, saturations, etc.