Efficient flow and transport simulations in reconstructed 3D pore geometries

Advances in Water Resources 33 (2010) 1508–1516

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Efficient flow and transport simulations in reconstructed 3D pore geometries Yan Zaretskiy a,⇑, Sebastian Geiger a, Ken Sorbie a, Malte Förster b a b

Institute of Petroleum Engineering, Heriot-Watt University, EH14 4AS Edinburgh, UK Fraunhofer Institute for Algorithms and Scientific Computing, Schloss Birlinghoven, D-53754 Sankt Augustin, Germany

a r t i c l e

i n f o

Article history: Received 21 April 2010 Received in revised form 18 August 2010 Accepted 22 August 2010 Available online 27 October 2010 Keywords: Pore-scale modeling Finite element Finite volume Solute transport Navier–Stokes equation Algebraic multigrid

a b s t r a c t Upscaling pore-scale processes into macroscopic quantities such as hydrodynamic dispersion is still not a straightforward matter for porous media with complex pore space geometries. Recently it has become possible to obtain very realistic 3D geometries for the pore system of real rocks using either numerical reconstruction or micro-CT measurements. In this work, we present a finite element–finite volume simulation method for modeling single-phase fluid flow and solute transport in experimentally obtained 3D pore geometries. Algebraic multigrid techniques and parallelization allow us to solve the Stokes and advection–diffusion equations on large meshes with several millions of elements. We apply this method in a proof-of-concept study of a digitized Fontainebleau sandstone sample. We use the calculated velocity to simulate pore-scale solute transport and diffusion. From this, we are able to calculate the a priori emergent macroscopic hydrodynamic dispersion coefficient of the porous medium for a given molecular diffusion Dm of the solute species. By performing this calculation at a range of flow rates, we can correctly predict all of the observed flow regimes from diffusion dominated to convection dominated. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Understanding the fundamental physics of fluid flow and transport in porous media is a crucial topic for all fields dealing with mass and heat transfer in subsurface systems, such as groundwater reservoirs, saline aquifers targeted for CO2 storage, oil and gas fields or geothermal fields. Continuum-scale models, which are used to describe flow at the scale of a reservoir, rely on various macroscopic parameters, such as the effective permeability, porosity or dispersivity. However, these macroscopic parameters in turn depend on the underlying microscopic properties of the system. They can be either measured experimentally or in some cases computed by means of numerical modeling of physical phenomena at the pore scale. The advantages of the second approach lie in the fact that, once the numerical algorithm is established, it can complement experimental data by providing an additional information and understanding about the phenomena under study. For instance, numerical modeling can show the details of fluid flow within separate pores of a rock sample and give insight into why a macroscopic parameter such as dispersivity or the extent of reacted grain surface area behaves in a certain way. Therefore, it is important to develop pore-scale models that allow us to predict the fluid flow properties directly from the information about the ⇑ Corresponding author. Tel.: +44 131 451 8299; fax: +44 131 451 3127. E-mail addresses: [email protected] (Y. Zaretskiy), sebastian.geiger@ pet.hw.ac.uk (S. Geiger), [email protected] (K. Sorbie), malte.foerster@scai. fraunhofer.de (M. Förster). 0309-1708/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2010.08.008

pore space geometry. One impediment to this upscaling process is the complexity of the geometry and topology of the pore system in realistic rocks. The field of pore-scale simulation has grown over the last decade [1] and applications include calculating the spatial distribution of the velocity field [2], solute transport [3,4], mixing and chemical reactions [5–7] and multi-phase flow [8–11]. Another important application of pore-scale modeling is related to upscaling and averaging in porous media. One example of this is deriving the relations between the underlying microscale variables and their averaged macroscale counterparts, such as the true pressure in the Navier– Stokes equation versus the average pressure in Darcy’s law [12]. Prediction of the physical properties of the porous media from their microscopic origins involves an accurate characterization of the detailed structure of the pore space. Direct non-destructive pore space measurements are now possible by the means of (synchrotron) X-ray computed microtomography [13,14]. In addition to this, statistical [15–19] and process-based [20] models for reconstructing 3D structure from 2D thin section images have been proposed. The two most wide-spread types of methods for the numerical simulations of single-phase transport processes at the pore scale are pore-network simulations [10] and Lagrangian particle-based methods, among which lattice Boltzmann methods [21] are the most popular ones. The lattice Boltzmann (LB) method is usually applied directly on the voxelized regular grid, but extensions to unstructured meshes are generally possible [22]. Pore-network modeling involves additional conversion of the measured pore

Y. Zaretskiy et al. / Advances in Water Resources 33 (2010) 1508–1516

space structure into a network of pore bodies interconnected by narrow pore throats of various cross-section shapes, which makes it possible to analyse samples containing a large amount of pores. These methods can be used to study various flow properties but suffer from certain drawbacks. As already mentioned, pore-network models do not preserve the original sample geometry which is important for some applications such as the reactive transport with the fluid–rock interactions where a particular goal is to understand how the reacted rock surface evolves over time. Standard LB methods with a single relaxation time and bounce-back boundary conditions are known to manifest a non-physical viscosity dependence of the permeability [23]. This shortage can be overcome by employing multiple-relaxation schemes but at the expense of extra computational costs. As stated above, LB method is not the only particle-based approach and other examples include Dissipative Particle Dynamic (DPD) [24], Smoothed Particle Hydrodynamic (SPH) [25] and Moving Particle Semi-implicit (MPS) [26]. Another approach is to directly solve the Navier–Stokes equation on a domain that represents a porous structure. For the single-phase applications this is the preferred approach because of its better numerical efficiency compared to all particle-based methods. Having said that, the main obstacle here is still the amount of computational work, which grows rapidly with the complexity of the domain geometry (i.e., the number of nodes required to correctly represent the domain boundaries increases) and with the dimensionality of the problem. This is because the inversion of the matrix representing the discretized equations for the geometry of interest is very costly, prohibiting the use of large pore samples or complex geometries. Therefore, this method is usually applied either in 2D [27] or in 3D with a simplified domain geometry, such as packs of spheres [28,29]. Previous work that deals with Navier–Stokes solution on experimentally measured pore geometries uses relatively small sample sizes [30–32]. An interesting approach was demonstrated in [33], where grain surfaces are introduced not through mesh conformance but with an immersed boundary method. Other methods were suggested that approximate the solution of the Navier–Stokes problem [34], but their comparison against the exact solution requires further analysis. The large body of research that deals with various methods for single-phase flow modeling at the pore-scale suggests an on-going interest in methods that link 3D images of porous rocks an numerical modeling of their various transport properties. In this work, we present a computationally efficient approach to simulate single-phase fluid flow and solute transport in porous media. It is based on the unstructured finite element–finite volume (FE–FV) method that resolves pore space geometry in great detail while solving the fundamental Navier–Stokes and advection–diffusion equations without invoking additional physical assumptions. Our FE–FV realization is made using the Complex System Modeling Platform, CSMP++ [35], our in-house object-oriented finite element–finite volume based C++ library tailored to simulating a wide range of fluid flow processes in geometrically complex porous media. The computational domain is discretized by an unstructured FE method and contains a complementary FV grid. Unstructured meshes with varying element sizes result in accurate capturing of grain boundary shapes while at the same time avoiding the overrefinement of internal regions of large pores. The meshes are constructed from X-ray computer tomography scans of real 3D porous samples, but could equally well be applied to any pore system model, e.g. from a 3D numerical reconstruction method [18,19]. In order to increase the efficiency of the calculations, we use CSMP++ in conjunction with an algebraic multigrid solver [36]. AMG methods have the distinctive advantage that they scale linearly as a function of the number of the unknowns.

1509

We solve the Stokes equation at the pore scale for the full FE grid to obtain the velocity field for single-phase flow, from which the permeability of the sample can be calculated (see below). With this velocity field we then solve the tracer solute transport equation for a species with molecular diffusion coefficient, Dm, in order to calculate the upscaled dispersion coefficient D. This procedure is repeated at a range of flow velocities to obtain D as a function of flow rate; more precisely we construct two plots: D/Ua versus Peclet number Pe = Ua/Dm and D/Dm versus Pe, where U is an average flow velocity and a is an average grain size. The first plot highlights the transition between various types of flow regime while the second one allows a quantitative comparison of sample dispersivities with other published results. Both results are compared with the previously reported experimental and numerical findings.

2. Methodology and implementation 2.1. Stokes equations The partial differential equations that govern the incompressible Newtonian fluid flow within an arbitrary domain X are the Navier–Stokes equation

q

@u ¼ lr2 u  qðu  rÞu  rp; @t

ð1Þ

and the law of mass conservation for an incompressible fluid

r  u ¼ 0:

ð2Þ

Here p, q, l, u are pressure, fluid density, fluid viscosity and velocity vector, respectively. In what follows, we consider the case of creeping flow, characterized by a very low Reynolds number (101). This condition occurs in most porous media flow problems and is well studied [37]. In this case the non-linear term in Eq. (1) vanishes. Additionally, in this work we focus on the steady-state flows and therefore we solve Eq. (1) without considering the time-derivative term. Overall, the governing Stokes equation describing laminar steady-state flow of a single fluid phase in a porous medium is

lr2 u ¼ rp:

ð3Þ

The boundary conditions, necessary to close Eqs. (2) and (3), include the no-slip condition at the rigid boundaries, i.e. ujCrigid ¼ 0, as well as the proper conditions at the inflow and outflow boundaries. Instead of assigning the inflow velocity profile, which can be considered as a more traditional type of Dirichlet boundary conditions, we follow the procedure described in [38], in order to impose the prescribed pressure drop across the porous sample. Namely, we set that

p  l@ n ujCin=out ¼ Pin=out

ð4Þ

for both the inflow and outflow boundaries Cin/out, where Pin/out are the desired pressure values and @ n denotes a spatial derivative in a direction normal to the corresponding boundary. The discrepancy between the intended and realized pressure drop due to the second term of the left-hand side of Eq. (4) is proportional to 1/r3 in a threedimensional case, where r is a curvature radius of the boundary surface. Hence, in our case of flat boundaries, this mismatch is cancelled. The steps necessary to construct the Galerkin finite element formulation of the Stokes equations are well-known and can be found, for instance, in [39]. We outline them in Appendix A and include an example benchmark showing the convergence of our numerical method.

1510

Y. Zaretskiy et al. / Advances in Water Resources 33 (2010) 1508–1516

2.2. Advection–diffusion equation Transport of an inert solute, denoted as C(x, t), within a fluid with an absence of volumetric sources and sinks is described by the general advection–diffusion equation

@Cðx; tÞ ¼ r  ðDm rCÞ  r  ðuCÞ; @t

ð5Þ

where x is a coordinate vector and Dm is a molecular diffusivity. The velocity u in each pore is directly given by the solution of the Stokes equation as described above. As a boundary condition, we keep C at a constant value at the inflow boundary and allow it to undergo pure advection at the outflow boundary,

Cðx; tÞjCin ¼ C 0 ;

rCðx; tÞ  njCout ¼ 0:

ð6Þ

Eq. (5) is solved by Godunov operator splitting, which assumes that the total time derivative can be calculated as a summation of the diffusion step and advection step. Applying the operator splitting technique has the advantage that we can use two different numerical methods to solve the diffusion and advection components, respectively. This numerical procedure is described in Appendix B. 2.3. Algebraic multigrid method and parallel discretization The 3D volumetric meshes, reconstructed from the porous sample CT measurements, are characterized by a large amount of nodes (of the order of millions) and result in an excessive size of a linear system corresponding to an employed numerical method. In our approach we use the algebraic multigrid library SAMG [40], developed at the Fraunhofer Institute for Algorithms and Scientific Computing. AMG solvers have the advantage that the CPU time and memory requirements scale with OðnÞ, where n is the number of unknowns, thus mitigating the difficulties associated with the domain size. This is a significant improvement over direct, ILU, or conjugate gradient solvers for which the numerical work increases with Oðna Þ where a > 1. The idea of AMG is to generalize the two main principles of geometric multigrid [41], error smoothing and coarse grid correction, directly to certain classes of sparse matrix equations. For the construction of the multilevel hierarchy no geometric information needs to be available. Coarse levels, inter-grid transfer, and coarse level operators are built automatically, based purely on algebraic information explicitly or implicitly contained in the underlying matrix (such as size and/or sign of the matrix entries). This makes

AMG especially suited for applications on unstructured grids where the implementation of geometric multilevel approaches is too complicated if possible at all. The SAMG package is not just an AMG implementation for one particular problem type but rather a complete multilevel framework. Several setup strategies are available that cover various types of coupled systems of PDEs. We found that in our particular case the solver works best if we use ILU factorization with no fill-in (ILU(0)) as a smoother in combination with a point-based coarsening approach and fluid pressure as a primary unknown. The latter means that the construction of coarser levels of the FE linear system from Stokes equation is performed based on the connectivity structure of the fluid pressure diagonal matrix blocks. Transfer operators are constructed for each unknown separately, following the standard AMG procedure for single-unknown problems [42,43]. Finally, SAMG is not used as a stand-alone solver. Each V-cycle rather serves as a preconditioner for GMRES(20) (GMRES restarted in every 20 steps) iterative solver. A second technique that we employ in our code to gain further speedup and overcome possible memory limitations for large-scale simulations is parallelization [44]. We use the node-level graph partitioner METIS [45] to automatically decompose the domain into a number of adjacent subdomains and each of them is then treated on a separate processor. The resulting partitions overlap by exactly one layer of elements that represent the communication volume between processors. Since every CPU contains only a portion of the whole computational domain, this procedure mitigates the memory requirements, particularly on 32-bit machines, by a factor which is almost linearly proportional to a number of partitions. 2.4. Mesh generation In order to investigate the flow through a realistic pore geometry, one has to transform it into a proper finite element mesh. The first stage of this process is performing an X-ray computer tomography scanning of an actual porous sample. This experiment results in a spatial X-ray absorption profile of the sample. Based on the absorption contrast between a fluid and a solid one can infer the location of pores and grains with an error depending on the scanner resolution and create a 3D binary image of the sample (see Fig. 1a) [14]. The next step is construction of a watertight isosurface (Fig. 1b), which isolates the region of interest, i.e. the pores from the solids. Finally the interior of the isosurface is filled with volumetric elements and this provides the FE mesh.

Fig. 1. (a) Voxelized representation of a Fontainebleau sandstone sample obtained by thresholding the corresponding X-ray CT data. Blue voxels denote solid space, brown ones denote pore space. Each voxel is a cube with the edges of 7.5 lm. The overall size of the sample is 1.5 mm in each direction. (b) FE-2 mesh reconstructed from the voxelized data. The outer surface contains 763 thousand of triangles. The surface interior is filled with 3.3 million of tetrahedral elements (not shown). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Y. Zaretskiy et al. / Advances in Water Resources 33 (2010) 1508–1516

There are numerous algorithms for isosurface extraction, with the standard approach being the so-called marching cubes algorithm [46]. However, this method is of limited applicability for our particular problem. Firstly, it produces topologically inconsistent non-manifold surfaces. A manifold surface has a property that a small enough sphere around any point belonging to it will be divided into two parts that lie on the opposite sides of the surface. Non-manifold surfaces do not share this property, which in our case means that their volumetric interior is not continuous everywhere. Another limitation of the marching cubes algorithm is that it produces surfaces with an excessively large amount of triangles, which can pose an obstacle for subsequent data processing by a volumetric mesher. In this work we use two different software packages to study the implications of a mesh construction step. The first one is the proprietary MimicsÒ software by Materialise, which is a complete package with various algorithms to translate raw data from X-ray CT scans into a format suitable for finite element meshing. MimicsÒ provides its own implementation of a modified marching cubes algorithm, which is guaranteed to produce manifold meshes with a smaller degree of fragmentation. It also has a mesh optimization functionality to improve the isosurface quality in terms of the aspect ratio of its triangles. The second is a 3D mesh generation module from the open source CGAL library [47]. The CGAL mesh generator avoids the isosurface stage and construct the tetrahedral elements complex directly from the voxelized data, allowing for almost entirely automated mesh generation. 3. Application We will now demonstrate the application of our numerical procedure to investigate the flow properties and emerging transport characteristics in the 3D Fontainebleau sample. 3.1. Reconstructed 3D meshes Fig. 1 demonstrates the application of a mesh generation workflow to a Fontainebleau sandstone sample. The sample has the dimensions of 1.5  1.5  1.5 mm3, a permeability of 1.19  1012 m2 and a porosity of 13.5% [48]. The CT scanner resolution is 7.5 lm, therefore the resulting 3D image of the sample has the size of 200 voxels in all three directions. The mesh generated with the CGAL library contains 3.1 million tetrahedra and 932 thousand nodes. The calculated porosity of this FE mesh is 13.2%, which is slightly lower because the CGAL mesher discards the regions of isolated pore space. We will refer to this mesh as FE-1.

1511

The isosurface from MimicsÒ contains 474 k triangles and 236 k vertices. The surface was consequently processed using the ANSYS ICEM CFDÒ software to generate the volumetric mesh, representing the pores, with 3.3 million tetrahedra and 763 k nodes. This mesh will be refered to as FE-2. The calculated porosity of the FE-2 is 17%. This discrepancy can be explained by the isosurface extraction process which used an option that reduces experimental noise but at the same time increases the diameter of the pores slightly. However, due to limited access to the MimicsÒ software, we were not able to modify the mesh once this difference has been found. Moreover, we use this mismatch as an instructive example of the sensitivity of pore-scale fluid flow predictions to the tiny changes of pore geometry. 3.2. Pressure and velocity profiles For calculating the velocity distribution, a pressure gradient of 1 Pa was applied across the sample in the x direction. The fluid had a uniform viscosity of 103 Pa s. No-slip conditions were assigned to the grain boundaries, as well as the domain boundaries parallel to the x-axis. The mesh was partitioned into 16 subdomains. The computing facility, used to run the calculation, was the Heriot-Watt University Cluster, which consists of 40 AMD Opteron 875 nodes, each supplied with 4 GB of RAM. Every computing node can run up to four jobs simultaneously. The peak memory requirement for our simulation was 6.5 GB. Hence we had to run the simulations in parallel, although they could be performed on a 64-bit desktop PC with 8 GB RAM equally well. Since the Heriot-Watt Cluster only allows for non-dedicated simulations, i.e. it is not possible to reserve entire CPU nodes and their RAM for a simulation, we used 16 processors to ensure that we always have sufficient memory available. The required CPU time was 2 min to solve the linear problem (A.2) with three million degrees of freedom. The scaling behaviour of SAMGp [44,49] suggests that the CPU time will increase linearly with the decrease in processors, i.e. approximately 4 min will be needed to solve the linear problem with eight processors, etc. This is still significantly faster compared to the LB method because we do not need to solve a transient problem. The results for pressure and velocity distribution within the pore space are shown in Fig. 2. Velocity is represented along the streamlines emerging from the inlet boundary. Streamlines effectively represent the preferential flow paths of an imaginary tracer, transported by the fluid. These paths do not cover the whole pore space, as the intricacy of pore geometry inevitably leads to some pores being stagnant, i.e. no flow occurs in them. Another consequence of the pore space structure is that the flow with the maxi-

Fig. 2. Fluid pressure (left) and velocity (right) distributions along the porous sample, obtained from the solution of the Stokes equations. Velocity is represented as streamlines.

1512

Y. Zaretskiy et al. / Advances in Water Resources 33 (2010) 1508–1516

mum velocity magnitude takes place not inside the narrowest throats, which is the case, for instance, in a steady flow through a single tube of variable radius. Instead the fluid favors wider and less resistant throats. In order to demonstrate these effects more clearly, Fig. 3 presents a complementary 2D velocity distribution, which corresponds to a photomicrograph of a grainstone sample with the dimensions of 10  6 mm2. Based on velocity distribution, one can estimate the resulting permeability of the sample by integrating the velocity vector over the outflow boundary Cout

R

k¼

lL u  n dCout Aðpout  pin Þ

ð7Þ

;

where A is the area of the sample cross-section perpendicular to the flow, L is sample length and pout and pin are the pressure values at the outflow and inflow boundaries, respectively. In our case these calculations yield a permeability value of 0.71  1012 m2 for the FE-1 mesh and 3  1012 m2 for the FE-2 mesh. Since these values differ by a factor of 4, we also ran a LB simulation on the voxelized representation of the same sample using the SHIFT scheme [50] and the same lateral boundary conditions. A relaxation time was taken equal to unity, which makes the permeability to be least affected by viscosity [23]. The resulting LB permeability is 0.66  1012 m2, which is very close to the result obtained with the FE-1 mesh. Permeability values are summarized in Table 1. This result allows us to make two conclusions. Firstly, the discrepancy between numerical values and an experimental permeability of 1.19  1012 m2 is not inherent to the FE modeling approach, but is rather a consequence of a limited sample size. Indeed, as demonstrated in [48], the REV size for this particular sample is VREV = 6.8 mm3, while we had access to a portion of data with a volume of Vp = 3.38 mm3. Since the sample is smaller than the REV scale, a no-slip boundary condition on the lateral sides of the domain that introduces an additional drag, may further reduce the permeability value, as demonstrated in [32]. Finally, as shown in [48], where the same sandstone sample was studied, when a domain exceeding the REV size is taken, the LB

method yields a permeability of 1.28  1012 m2, a result much closer to the experimental value. Secondly, the permeability mismatch demonstrates that the isosurface reconstruction step has to be treated with care. Though an isosurface generation can be seen as an attempt to restore the geometry information from the jagged voxelized CT data, it is a post-processing procedure and can introduce geometrical errors. In our case the noise reduction procedure of MimicsÒ effectively reduced the CT resolution while the CGAL isosurface was a byproduct of the 3D mesh generation and as such was conformant to the CT data with a minimal degree of smoothing. It should be noted that the voxelized X-ray CT data is already an approximate representation of the original pore geometry and the accuracy of the computed permeability will not only depend on the accuracy of the numerical method but also on the quality of the CT data, the segmentation and smoothing algorithms and other factors. 3.3. Solute transport The modeled single-phase velocity distribution can be used to calculate the hydrodynamic dispersivity of the porous sample, which describes the effective parameter to model solute transport at the macro, i.e. continuum, scale. Dispersivity in porous media results from the interplay of Brownian motion of the solute molecules, i.e. molecular diffusion, and the velocity differences among different streamlines, i.e. mechanical dispersion. In this section we will demonstrate how our model can capture the effects of this interplay consistently with experimental findings [51]. In order to estimate the sample dispersivity, we first solve the advection–dispersion equation, using the procedure outlined in Section 2.2. From this we obtain a spatial distribution of a solute concentration C(x, t) at each time step. We can therefore calculate the average breakthrough concentration CðtÞ by integrating C(x, t) across the outflow boundary Cout, weighted with a local flux value

R CðtÞ ¼

Cðx; tÞ uðxÞ dCout R : uðxÞ dCout

ð8Þ

We then fit the resulting breakthrough curve against the Ogata– Banks analytical solution of a 1D advection–dispersion equation [52]

C OB ðtÞ ¼

Fig. 3. Velocity magnitude distribution across a 2D model of a grainstone sample. Region 1 represents a stagnant area where velocity is low. Region 2 demonstrates the development of a preferential flow path.

Table 1 Permeability values obtained from different numerical models. The experimentally measured permeability is 1.19  1012 m2. Model

Permeability value (m2)

FE-1 FE-2 LB

0.71  1012 3.0  1012 0.66  1012

       Co L  Ut LU L þ Ut erfc pffiffiffiffiffiffiffiffi ; erfc pffiffiffiffiffiffiffiffi þ exp D 2 4Dt 4Dt

ð9Þ

where Co is a fixed inflow concentration, U is the mean flow velocity and D is the (longitudinal) hydrodynamic dispersivity. The abovementioned interplay between convection and molecular diffusion is expressed in terms of a throat Peclet number Pe = Ua/Dm, where a is a characteristic length of the sample. In this particular case we took a equivalent to a mean Fontainebleau sandstone grain size of 250 lm. Low Peclet numbers correspond to the diffusion-dominated flow regime, while high ones are characteristic of convection-dominated regime. Experimental studies [51] suggest that in the first regime dispersivity decreases with velocity if all other parameters are equal. In the latter regime it increases slowly with velocity. By changing the imposed pressure gradient and the molecular diffusivity we can verify this dependency between the Peclet number and an upscaled sample dispersivity. Fig. 4 demonstrates two breakthrough curves that correspond to the diffusive (Pe = 0.18) and convective (Pe = 55) flow regimes, respectively. Though both curves are well approximated by an analytical solution (9), a close inspection of the second (Pe = 55) curve reveals that the concentration does not reach its maximum inflow level within the time span of the modeling. As indicated in [53], the size necessary for a breakthrough curve to have no tailing, is linearly proportional to the Peclet number. When the solute diffuses

Y. Zaretskiy et al. / Advances in Water Resources 33 (2010) 1508–1516

Fig. 4. Flux-weighted breakthrough curves, corresponding to Peclet numbers Pe = 0.18 (a) and Pe = 55 (b). High Peclet numbers characterize a convectiondominated regime in which case long tailing occurs. That is, the concentration takes a longer time to reach the maximum value than predicted by the analytical solution.

away from the major flow paths it becomes trapped or delayed within stagnant zones [54]. If the solute diffusivity is low (i.e. Peclet number is high), it takes a significant amount of time for the solute to diffuse back to high-velocity streamlines and hence causes the tailing effect. In case of higher (molecular) diffusivity, e.g. thermal diffusion, this heterogeneity effect becomes less dominant and the tailing is reduced. Fig. 5a shows the calculated FE-2 sample dispersivities for a wide range of throat Peclet numbers. The results are plotted in a dimensionless fashion as log10(D/Ua) versus log10(Pe) [3,4,51,55]. This way points corresponding to different values of molecular diffusivity collapse onto a single curve. The overall form of this curve compares very well with the experimental results [51], which demonstrates that our code captures this well-known emergent behavior. At values of Pe < 1, the curve has a negative slope because dispersion is almost constant in the diffusion-dominated regime. A minimum in the value of D/Ua is predicted in the range 1 < Pe < 10 which is also found experimentally [51]. For Pe > 10 the slope turns positive corresponding to the regime where mechanical dispersion is prevalent but where molecular diffusion cannot be ignored. At Pe values greater than 1000 the slope approaches zero. This is where dispersion is purely mechanical and D is, theoretically, linear to U [56]. The observed decline at Pe > 103 occurs because the tailing of breakthrough curves leads to a poorer fit of numerical calculations against Eq. (9) and hence a deviation of D/Ua from the expected behavior. It should also be noted that points in Fig. 5 were calculated by varying both velocity and molecular diffusivity to verify the consistency of the dimensionless plot, i.e. that the cases corre-

1513

Fig. 5. (a) Dispersivity versus throat Peclet number Pe = Ua/Dm plotted in a dimensionless fashion as D/Ua. The shape of the curve reflects the transition from the diffusion- to convection-dominated regime of the fluid flow. (b) Dispersivity versus throat Peclet number Pe = Ua/Dm plotted as D/Dm. Two curves correspond to conformant (FE-1) and slightly inflated (FE-2) meshes.

sponding to different diffusivities but same Peclet number result in the same value of D/Ua. The curve in Fig. 5a emphasizes the transition between three different flow regimes and for that reason we only plot the results for the FE-2 mesh. Another way to plot the transport modeling results is in the log10(D/Dm) versus log10(Pe) scale, which allows easier quantitative comparison with other works dealing with the porous media dispersion. This plot is presented in Fig. 5b for both FE-1 and FE-2 meshes. It should be noted that in both cases we used the same average grain size value a for calculating the Peclet number which is justified by the fact that the relative difference of this parameter between FE-1 and FE-2 is close to one and will only result in a minor shift of the curve. Fig. 5b agrees well in terms of dispersivity magnitude with the results on longitudinal dispersion compiled in [3]. At low Peclet numbers diffusion drives the solute propagation and therefore the effective dispersivity is lower than molecular diffusivity (log10(D/Dm) is negative) as the pore space restricts the solute movement. The slope of the curve in this region gradually approaches zero as the very low Pe dispersivity depends only on the formation resistivity factor and porosity [57]. Again, at the mixed-flow regime convection–diffusion interplay is the most pronounced and hence the curve has a slope higher than unity. The difference between FE-1 and FE-2 manifests itself in both the dispersivity magnitudes and the power law coefficient d between D/Dm and Pe (a slope in the logarithmic scale). The results are d = 1.23 for FE-1 and d = 1.1 for FE-2 for the 10 < Pe < 200 region. This further illustrates the sensitivity of pore-scale modeling results to the proper construction of the numerical domain. The power law coefficient describes the influence of diffusion from the low-velocity boundary layers on the overall sample dispersiv-

1514

Y. Zaretskiy et al. / Advances in Water Resources 33 (2010) 1508–1516

ity. The widening of pore channels in FE-2 mesh reduced the rate of this diffusion and thus resulted in a lower value of d. Other studies on numerically reconstructed Fontainebleau sandstones report the power law coefficient of d = 1.34  1.56 [53]. For the sake of comparison, published experimental results on power law coefficient in Berea sandstones range from d = 1.13 [58] to d = 1.33 [59], which suggests that both values computed in this work are physically quite reasonable. Finally, for the FE-1 sample we observed that at around Pe > 200 the breakthrough curves demonstrate a strong non-Gaussian behavior with long tailing and that they can no longer be approximated with the Ogata–Banks solution. This is in agreement with the discussion on non-Gaussian dispersion behavior in [53]. We believe that it occurs earlier for FE-1 than FE-2 in terms of Peclet numbers because the FE-1 sample has narrower pore throats and hence contains more streamlines that do not reach the outflow boundary. In the case of small diffusivity when a solute particle starts moving along a dead-end streamline there is a low probability of it reaching the other side of a sample and therefore outflow concentration rises very gradually. One important application of the procedure described here is an analysis of relation between the diffusivity of a solute in bulk water and the diffusivity of the same solute in saturated porous media. This relation is crucial in the performance assessment of underground disposal systems for the radioactive waste. The average macroscopic diffusion coefficient is defined as Da = /d/s2  Dw [60], where / is diffusion-accessible porosity, d is constrictivity, s is tortuosity and Dw is solute diffusivity in bulk water. Porosity and constrictivity can be determined from the geometrical structure of the porous sample, while tortuosity can be obtained by calculating the stationary velocity profile and the corresponding streamlines [61]. The average diffusion coefficient can be computed by tracking the spatial evolution of the solute concentration within the sample when no advection occurs. Hence the validity and the range of applicability of the aforementioned relation can be studied.

4. Concluding remarks In this work we presented a direct pore-scale procedure for modeling the single-phase fluid flow and transport in reconstructed 3D porous medium. The key elements of this procedure are: (1) a FE solution of the stationary Stokes equation on unstructured meshes for an accurate representation of domain geometry; (2) an algebraic multigrid method for solving the resulting FE linear problem, which offers a linear number-of-operations scalability with respect to mesh refinement and hence significantly mitigates the computer resources limitation factor; (3) an optional domain decomposition for problem parallelization; (4) a combined FE–FV technique for solving the transient inert transport problem. Despite the fact that these components are not novel, when combined together they offer a powerful viable alternative to other methods for calculating single-phase flow in reconstructed pore geometries [33,62]. By applying this procedure to a reconstructed Fontainebleau sandstone sample we were able to demonstrate that the FE solution yields the same permeability value as the LB method while requiring a fraction of time necessary for LB solution to converge. We also highlighted the detrimental effect of reducing the amount of geometrical data by means of excessive domain boundary smoothing. Transport modeling demonstrated the correct transition between various flow regimes from diffusion- to convectiondominated and gave quantitatively similar results when upscaled sample dispersivity was compared to other published results. It also showed a slight decrease in the sample dispersivity as well as in the relative contribution of diffusion in the mixed-flow re-

gime in the case of a mesh with inflated pore throats. Finally, we observed that a mesh reconstruction procedure determines the onset of the breakthrough curves tailing with respect to the Peclet number. Acknowledgements We thank the following for the financial support for this work: the ExxonMobil Research Alliance ‘‘Fundamental Controls of Flow in Carbonates” (FC)2, the UK EPSRC Dorothy Hodgkin Programme, and the Edinburgh Collaborative of Subsurface Science and Engineering, a joint research institute of the Edinburgh Research Partnership in Engineering and Mathematics. We also thank Materialise company for providing us with the opportunity to work with MimicsÒ software free of charge. Finally, we are grateful to Jingsheng Ma for helping us to calculate the LB permeability value with his SHIFT code. This manuscript benefited greatly from the careful reviews by Rainer Helmig and five other anonymous reviewers. Appendix A. Stokes discretization Pressure and velocity components are approximated as

uhi ðxÞ

¼

X

NA ðxÞuiA ;

ph ðxÞ ¼

nodes

X

NA ðxÞpA ;

ðA:1Þ

nodes

where NA is a shape function associated with global node number A, and uiA and pA are the values of uhi ðxÞ and ph(x) at a node number A. In the general case velocity and pressure nodes, as well as the shape functions, can be different. In order to reduce the amount of computational work and use the same first order of approximation for both uhi ðxÞ and ph(x), we implement the Galerkin/least-squares (GLS) method [63]. This method exhibits stable and convergent approximation for any choice of the discrete velocity and pressure interpolations. Existence, uniqueness and error behavior of the resulting finite element solution has been demonstrated [64]. We substitute Eq. (A.1) into Eqs. (2) and (3), multiply the equations with weighting functions (identical to NA(x) in the standard Galerkin FE formulation) and integrate over the whole domain. This procedure yields the following system of linear equations

2

Kl

6 6 0 6 6 0 4 GTx

0

0

Kl 0

0 Kl

GTy

GTz

3 2 3 2b3 f cx u 6 x7 7 7 6 Gy 7 6 c f uy 7 6 by 7 7 76 6 7 ¼ 6 7: Gz 7 4 c bz 7 uz 5 6 5 f 4 5 b Kst p b 0 Gx

ðA:2Þ

The hat sign indicates a vector of nodal values, for example

  cx ¼ ux1 ux2    uxn T : u

ðA:3Þ

Every entry on the left-hand side matrix of (A.2) is a submatrix of the size of n  n, where n is the number of mesh points. Each element’s contribution to these submatrices is calculated as follows: e

Kl ¼ l Kest ¼ Gexi ¼

Z Xe

Z Xe

Z X

e

se

! b T @N b X @N dV; @xi @xi i ! b T @N b X @N dV; @x @x i i i

bT @N b dV: N @xi

b is a vector containing nodal shape Here Xe is an element volume, N b ¼ ½N 1 N2 N3 . . . Nn . Symbols u cxi and p b denote nodal functions, N unknowns of velocity components and pressure and fc xi are righthand side vectors containing the corresponding boundary conditions. Finally, submatrix Kst comes from the GLS stabilization technique

Y. Zaretskiy et al. / Advances in Water Resources 33 (2010) 1508–1516

1515

Fig. B.7. A simple 2D finite element mesh (dashed) with the corresponding and virtual node-centered finite volume mesh (solid). Mesh nodes are marked as filled black circles.

Fig. A.6. Velocity and pressure error in L2-norm versus characteristic element size h. Velocity field demonstrates the linear convergence while pressure converges sublinearly due to the implemented GLS stabilization technique.

where se is an element-wise stabilization parameter, which is best 2 chosen as se ¼ ahe =4l [65]. In this definition he is a measure of element size and a is a constant of the order of unity. The choice a = 2/3 appears to be optimal for linear elements. A.1. Benchmarking against analytical solution

Z

In order to verify the convergence of our FE scheme, we benchmarked it against an analytical solution for the stationary incompressible Stokes flow from [39]. The problem consists of determining the velocity field u = (ux, uy) and the pressure p in the 2D square domain X = [0, 1]  [0, 1] with a prescribed body force within X and a no-slip boundary condition on @ X. Fig. A.6 shows that the numerical error, measured in the L2 norm, decreases as the mesh is refined. As mentioned before, we use linear elements for both velocity and pressure and therefore expect that the error decreases as OðhÞ. Numerical results show that the error for the velocity field decreases at this rate while the pressure converges at half the rate. This is due to the pressure diffusion introduced in the stabilization term. The same difference in convergence rates for velocity and pressure would hold for a more conventional mixed FE method. Appendix B. ADE discretization The diffusion derivative, with the assumption of constant diffusivity and absence of volumetric sinks and sources, is given by

  @C ¼ Dm r2 C; @t diff

ðB:1Þ

which we solve using the FE method. The FE procedures necessary for solving the diffusion equation are well established [66] and result in the following linear system

 M b tþDt ¼ M C bt:  KD C Dt Dt



ðB:2Þ

The element contribution matrices are given by

KeD ¼ Dm Me ¼

Z Xe

Z Xe

! b T @N b X @N dV; @xi @xi i

bTN b dV: N

The advection derivative is given by

  @C ¼ r  ðuCÞ; @t adv

which we solve using the mass and shock preserving finite volume method. The finite volume mesh is constructed as a complement to the FE mesh around the corner nodes of each finite element [67,68]. The finite volume faces are constructed by connecting the finite element barycenters with the barycenters of its faces and edges, thus forming a cell around each node. A corresponding 2D FE–FV mesh representation is shown in Fig. B.7. We integrate Eq. (B.3) over a finite volume Vi, corresponding to a node i, and apply the Gauss–Ostrogradsky theorem

ðB:3Þ

Vi

@CðtÞ dV ¼  @t

I

n  u CðtÞ dA:

ðB:4Þ

Ai

Here Ai is the boundary area of the finite volume Vi and n is the outward pointing normal vector of a finite volume boundary. The area integral on the right-hand side of Eq. (B.4) can be represented as a sum over all FV faces, marked with index j (see Fig. B.7). By discretizing the time derivative on the left-hand side of Eq. (B.4) using a backward Euler scheme, we obtain the following formulation

C itþDt ¼ C ti 

Dt X ðnj  uj Þ C jtþDt Aj : Vi j

ðB:5Þ

A value of solute concentration at a jth face is taken to be equal to a value at an upstream node, i.e. a node for which (njuj) > 0. Value of velocity uj is calculated by interpolation of the nodal values onto the barycenter of an element which the face j belongs to. We also use the implicit formulation because the presence of very small elements in the meshes that represent porous samples leads to an excessively restricting CFL criterion. For meshes used in this work, it would require about 105 time steps in order to fully saturate them with an inert solute. With the implicit method it was possible to reduce the amount of time steps 100 times without sacrificing numerical accuracy. For all cases when the ADE was solved in this work, a maximum acceptable time step was empirically calculated such that the evolution of the concentration profile remained the same within a desired tolerance when the time step was halved. References [1] Meakin P, Tartakovsky AM. Modeling and simulation of pore-scale multiphase fluid flow and reactive transport in fractured and porous media. Rev Geophys 2009;47(3):RG3002. [2] Narsilio GA, Buzzi O, Fityus S, Yun TS, Smith DW. Upscaling of Navier–Stokes equations in porous media: theoretical, numerical and experimental approach. Comput Geotech 2009;36(7):1200–6. [3] Bijeljic B, Muggeridge AH, Blunt MJ. Pore-scale modeling of longitudinal dispersion. Water Resour Res 2004;40:W11501. [4] Acharya RC, Van Dijke MIJ, Sorbie KS, Van der Zee S, Leijnse A. Quantification of longitudinal dispersion by upscaling Brownian motion of tracer displacement in a 3D pore-scale network model. Adv Water Resour 2007;30(2):199–213. [5] Li L, Peters CA, Celia MA. Upscaling geochemical reaction rates using pore-scale network modeling. Adv Water Resour 2006;29(9):1351–70.

1516

Y. Zaretskiy et al. / Advances in Water Resources 33 (2010) 1508–1516

[6] Tartakovsky AM, Neuman SP. Effects of Peclet number on pore-scale mixing and channeling of a tracer and on directional advective porosity. Geophys Res Lett 2008;35(21):L21401. [7] Willingham TW, Werth CJ, Valocchi AJ. Evaluation of the effects of porous media structure on mixing-controlled reactions using pore-scale modeling and micromodel experiments. Environ Sci Technol 2008;42(9):3185–93. [8] Celia MA, Reeves PC, Ferrand LA. Recent advances in pore scale models for multiphase flow in porous media. Rev Geophys 1995;33:1049–58. [9] Aker E, Måløy KJ, Hansen A. Simulating temporal evolution of pressure in twophase flow in porous media. Phys Rev E 1998;58(2):2217–26. [10] Blunt MJ. Flow in porous mediapore-network models and multiphase flow. Curr Opin Colloid Interf Sci 2001;6(3):197–207. [11] Knackstedt MA, Sheppard AP, Sahimi M. Pore network modelling of two-phase flow in porous rock: the effect of correlated heterogeneity. Adv Water Resour 2001;24(3–4):257–77. [12] Nordbotten JM, Celia MA, Dahle HK, Hassanizadeh SM. Interpretation of macroscale variables in Darcy’s law. Water Resour Res 2007;43(8):W08430. [13] Spanne P, Thovert JF, Jacquin CJ, Lindquist WB, Jones KW, Adler PM. Synchrotron computed microtomography of porous media: topology and transports. Phys Rev Lett 1994;73(14):2001–4. [14] Dunsmuir JH, Ferguson SR, D’Amico KL, Stokes JP. X-ray microtomography: a new tool for the characterization of porous media. SPE Paper SPE-22860. In: SPE annual technical conference and exhibition; 1991. [15] Adler PM, Jacquin CG, Quiblier J. Flow in simulated porous media. Int J Multiphase Flow 1990;16(4):691–712. [16] Roberts AP. Statistical reconstruction of three-dimensional porous media from two-dimensional images. Phys Rev E 1997;56(3):3203–12. [17] Okabe H, Blunt MJ. Prediction of permeability for porous media reconstructed using multiple-point statistics. Phys Rev E 2004;70(6):66135. [18] Wu K, Nunan N, Crawford JW, Young IM, Ritz K. An efficient Markov chain model for the simulation of heterogeneous soil structure. Soil Sci Soc Am J 2004;68(2):346. [19] Wu K, Van Dijke MIJ, Couples GD, Jiang Z, Ma J, Sorbie KS, et al. 3D stochastic modelling of heterogeneous porous media–applications to reservoir rocks. Transport Porous Med 2006;65(3):443–67. [20] Øren PE, Bakke S. Process based reconstruction of sandstones and prediction of transport properties. Transport Porous Med 2002;46(2):311–43. [21] Chen S, Doolen GD. Lattice Boltzmann method for fluid flows. Ann Rev Fluid Mech 1998;30(1):329–64. [22] Peng G, Xi H, Duncan C, Chou SH. Lattice Boltzmann method on irregular meshes. Phys Rev E 1998;58(4):4124–7. [23] Pan C, Luo LS, Miller CT. An evaluation of lattice Boltzmann schemes for porous medium flow simulation. Comput Fluids 2006;35(8–9):898–909. on LB dependence on the relaxation parameter. [24] Hoogerbrugge PJ, Koelman J. Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. EPL (Europhys Lett) 1992;19:155. [25] Lucy LB. A numerical approach to the testing of the fission hypothesis. Astron J 1977;82:1013–24. [26] Koshizuka S, Tamako H, Oka Y. A particle method for incompressible viscous flow with fluid fragmentation. Comput Fluid Dynam J 1995;4(1):29–46. [27] Garmeh G, Johns R, Lake L. Pore-scale simulation of dispersion in porous media. SPE J 2009;14(4):559–67. [28] Cardenas MB. Three-dimensional vortices in single pores and their effects on transport. Geophys Res Lett 2008;35(18):L18402. [29] Cardenas MB. Direct simulation of pore level Fickian dispersion scale for transport through dense cubic packed spheres with vortices. Geochem Geophys Geosyst 2009;10(12):Q12014. [30] Fourie W, Said R, Young P, Barnes DL. The simulation of pore-scale fluid flow with real world geometries obtained from X-ray computed tomography. In: Proceedings of the COMSOL conference; 2007. [31] Malinouskaya I, Mourzenko VV, Thovert JF, Adler PM. Wave propagation through saturated porous media. Phys Rev E 2008;77(6):66302. [32] Gerbaux O, Buyens F, Mourzenko VV, Memponteil A, Vabre A, Thovert JF, et al. Transport properties of real metallic foams. J Colloid Interf Sci 2009. [33] Smolarkiewicz PK, Larrabee Winter C. Pores resolving simulation of Darcy flows. J Comput Phys 2010. [34] Akanji LT, Matthäi SK. Finite element-based characterization of pore-scale geometry and its impact on fluid flow. Transport Porous Med 2010;81(2): 241–59. [35] Matthäi SK, Geiger S, Roberts SG, Paluszny A, Belayneh M, Burri A, et al. Numerical simulation of multi-phase fluid flow in structurally complex reservoirs. Geol Soc London Special Publ 2007;292(1):405. [36] Stüben K. A review of algebraic multigrid. J Comput Appl Math 2001;128: 281–309. [37] Happel J, Brenner H. Low Reynolds number hydrodynamics: with special applications to particulate media. Springer; 1983.

[38] Heywood JG, Rannacher R, Turek S. Artificial boundaries and flux and pressure conditions for the incompressible Navier–Stokes equations. Int J Numer Methods Fluids 1996;22(5):325–52. [39] Donea J, Huerta A. Finite element methods for flow problems. John Wiley & Sons Inc.; 2003. [40] Stüben K, Clees T. SAMG users manual; 2005. [41] Wesseling P, Oosterlee CW. Geometric multigrid with applications to computational fluid dynamics. J Comput Appl Math 2001;128(1-2):311–34. [42] Brandt A, McCormick SF, Ruge JW. Algebraic multigrid (AMG) for automatic multigrid solutions with application to geodetic computations. Report. Fort Collins, Colorado: Institute for Computational Studies; 1982. [43] Stüben K. Albebraic multigrid (AMG): an introduction with applications. GMD Report 70; 1999. [44] Coumou D, Matthäi SK, Geiger S, Driesner T. A parallel FE–FV scheme to solve fluid flow in complex geologic media. Comput Geosci 2008;34(12):1697–707. [45] Karypis G, Kumar V. MeTiS – a software package for partitioning unstructured graphs, partitioning meshes, and computing fill-reducing orderings of sparse matrices – version 4.0. University of Minnesota; 1998. p. 102. [46] Cline HE, Lorensen WE. Marching cubes: A high resolution 3-D surface construction algorithm. Comput Graph 1987;21(4):163–9. [47] Rineau L, Tayeb S, Yvinec M. 3D mesh generation. In: CGAL user and reference manual, 3.6 edition, CGAL editorial board; 2010. . [48] Jiang Z. Quantitative characterization of the geometry and topology of pore space in 3D rock images. PhD thesis. Heriot-Watt University; 2008. [49] Krechel A, Stüben K. Parallel algebraic multigrid based on subdomain blocking. Parallel Comput 2001;27(8):1009–31. [50] Ma J, Wu K, Jiang Z, Couples GD. SHIFT: an implementation for lattice Boltzmann simulation in low-porosity porous media. Phys Rev E 2010;81(5):56702. [51] Fried JJ, Combarnous MA. Dispersion in porous media. Adv Hydrosci 1971;7:169–282. [52] Ogata A, Banks RB. A solution of the differential equation of longitudinal dispersion in porous media. Prof. Paper 411-A, US Geological Survey, US Government Printing Office, Washington, DC; 1962. [53] Salles J, Thovert JF, Delannay R, Prevors L, Auriault JL, Adler PM. Taylor dispersion in porous media. Determination of the dispersion tensor. Phys Fluids A: Fluid Dynam 1993;5:2348. [54] Haggerty R, Gorelick SM. Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. Water Resour Res 1995;31(10):2383–400. [55] Sorbie KS, Clifford PJ. Inclusion of molecular diffusion effects in the network modelling of hydrodynamic dispersion in porous media. Chem Eng Sci 1991;46(10):2525–42. [56] Brenner H. Dispersion resulting from flow through spatially periodic porous media. Philos Trans Roy Soc London Ser A, Math Phys Sci 1980;297(1430):81–133. [57] Sahimi M, Islam MR. Flow and transport in porous media and fractured rock. Wiley-VCH, Verlag GmbH; 1995. [58] Gist GA, Thompson AH, Katz AJ, Higgins RL. Hydrodynamic dispersion and pore geometry in consolidated rock. Phys Fluids A: Fluid Dynam 1990;2:1533. [59] Kinzel DL, Hill GA. Experimental study of dispersion in a consolidated sandstone. Can J Chem Eng 1989;67(1):39–44. [60] Skagius K, Neretnieks I. Porosities and diffusivities of some nonsorbing species in crystalline rocks. Water Resour Res 1986;22(3):389–98. [61] Matyka M, Khalili A, Koza Z. Tortuosity–porosity relation in porous media flow. Phys Rev E 2008;78(2):26306. [62] Ovaysi S, Piri M. Direct pore-level modeling of incompressible fluid flow in porous media. J Comput Phys 2010;229(19):7456–76. [63] Hughes TJR, Franca LP. A new finite element formulation for computational fluid dynamics. VII: The Stokes problem with various well posed boundary conditions: symmetric formulations that converge for all velocity, pressure spaces. Comput Methods Appl Mech Eng 1987;65(1):85–96. [64] Zhou TX, Feng MF. A least squares Petrov–Galerkin finite element method for the stationary Navier–Stokes equations. Math Comput 1993;60(202):531–43. [65] Codina R. On stabilized finite element methods for linear systems of convection–diffusion–reaction equations. Comput Methods Appl Mech Eng 2000;188(1-3):61–82. [66] Istok J. Groundwater modeling by the finite element method. American Geophysical Union; 1989. [67] Geiger S, Roberts S, Matthäi SK, Zoppou C, Burri A. Combining finite element and finite volume methods for efficient multiphase flow simulations in highly heterogeneous and structurally complex geologic media. Geofluids 2004;4(4):284–99. [68] Paluszny A, Matthai SK, Hohmeyer M. Hybrid finite element–finite volume discretization of complex geologic structures and a new simulation workflow demonstrated on fractured rocks. Geofluids 2007;7(2):1–23.