Development of a numerical workflow based on μ-CT imaging for the determination of capillary pressure – saturation-specific interfacial area relationship in 2-phase flow pore-scale porous-media systems : a case study on Heletz sandstone

In this case study, we present the implementation of a finite element method (FEM)-based numerical porescale model that is able to track and quantify the propagating fluid–fluid interfacial area on highly complex microcomputed tomography (μ-CT)-obtained geometries. Special focus is drawn to the relationship between reservoir-specific capillary pressure (pc), wetting phase saturation (Sw) and interfacial area (awn). The basis of this approach is highresolutionμ-CT images representing the geometrical characteristics of a georeservoir sample. The successfully validated 2-phase flow model is based on the Navier–Stokes equations, including the surface tension force, in order to consider capillary effects for the computation of flow and the phase-field method for the emulation of a sharp fluid–fluid interface. In combination with specialized software packages, a complex high-resolution modelling domain can be obtained. A numerical workflow based on representative elementary volume (REV)-scale pore-size distributions is introduced. This workflow aims at the successive modification of model and model set-up for simulating, such as a type of 2-phase problem on asymmetricμ-CT-based model domains. The geometrical complexity is gradually increased, starting from idealized pore geometries until complex μ-CT-based pore network domains, whereas all domains represent geostatistics of the REV-scale core sample pore-size distribution. Finally, the model can be applied to a complex μ-CT-based model domain and the pc–Sw–awn relationship can be computed.


Introduction
Understanding the evolution of the fluid-fluid interfaces in 2-phase flow in porous-media systems is relevant for a series of engineering and technological applications (e.g.carbon capture and storage (CCS), nuclear waste repository, oil recovery; Hassanizadeh and Gray, 1990;Joekar-Niasar et al., 2008;Niessner and Hassanizadeh, 2008;Reeves and Celia, 1996).The quantity of the (typically unknown parameter) interfacial area between two phases restricts kinetic interphase mass and energy transfer (Ahrenholz et al., 2011).In numerical modelling with classical macro-scale flow models, these processes are not properly taken into account (Muccino et al, 1998;Ahrenholz et al., 2011).Here, averaged quantities neglecting the interfacial area are often used (Tatomir et al., 2015), which are described to lead to negligence of kinetics of mass transfer, assuming local equilibrium (Niessner and Hassanizadeh, 2009).In order to resolve this, the capillary pressure(p c )-wetting phase saturation(S w )interfacial area(a wn ) relationship using the concept of spe-A.Peche et al.: Development of a numerical workflow based on µ-CT imaging cific interfacial area per volume described by Hassanizadeh and Gray (1990) can be used.
Based on an extended form of Darcy's law (considering fluid-fluid friction force and interfacial forces), the mathematical fundamentals of pore-scale 2-phase flow regarding the p c -S w -a wn relationship, are given in Hassanizadeh and Gray (1980, 1990, 1993) and Niessner and Hassanizadeh (2008).The derivation of this porous-mediaspecific p c -S w -a wn relationship can be realized using physical experiments e.g. in the form of micromodels as described amongst others in Karadimitriou et al. (2013Karadimitriou et al. ( , 2014)), Karadimitriou and Hassanizadeh (2012) and Lenormand and Touboul (1988).Corresponding numerical models are given, amongst others, in Ahrenholz et al. (2011) and Porter et al. (2009).These specific approaches succeed in deriving the p c -S w -a wn relationship, but are based on Lattice-Boltzmann simulations, which are considered to be elegant but are limited by computational resources (White et al., 2006).Generally, the derivation of this parametric relationship on the pore scale requires a high computational effort as well as a complex numerical model.Arzanfudi and Al-Khoury (2015) present a finite element method (FEM)-based model that is able to simulate CO 2 leakage through a wellbore.In that one-dimensional approach, the moving CO 2 -water interface is successfully traced using the level-set method.Akhlaghi Amiri and Hamouda (2013) evaluate two approaches for dynamic interface tracking: the level-set method and the phasefield method for two-dimensional FEM-based modelling of 2-phase flow with viscosity contrast in a dual-permeability porous medium.Both approaches are able to quantify the specific interfacial area between the fluids.They conclude that the phase-field method is more successful in complicated porous media, gives more realistic results regarding pressure gradients and fluid profiles and they observed less computation time compared with the level-set method.
Several experimental and numerical approaches were developed to measure the interfacial area between the two fluid phases.Current research in the field of CCS focuses on the design and development of field investigation techniques allowing short response times in on-site plume monitoring and detection.Detection and quantification of the CO 2 -brine fluid-fluid interface could be realized using kinetic interface sensitive (KIS) tracers as described in Schaffer et al. (2013).Several mathematical models for immiscible 2-phase flow and KIS tracer transport in porous media were developed in Tatomir et al. (2013Tatomir et al. ( , 2014) ) and Tatomir et al. (2015), whereas the latter also considers compositional effects.In Tatomir et al. (2015) a 2-phase four -component flow and transport model is realized with a kinetic mass transfer of tracers between the two fluids, taking the dissolution of CO 2 and brine into account.The p c -S w -a wn relationship is seen as an extension to the standard Brooks-Corey model (Brooks and Corey, 1964).
The modelling approaches by Tatomir et al. (2013Tatomir et al. ( , 2014) and Tatomir et al. (2015) require the a priori knowledge of the fluid-fluid-solid system-specific p c -S w -a wn constitutive relationship (Tatomir et al., 2013).
The present article presents the corresponding model developed for deriving this specific p c -S w -a wn constitutive relationship.In this study, a finite element method (FEM)based pore-scale model is used in order to determine the dynamic evolution of the CO 2 -water interface in geometries with a gradually increasing level of complexity.Within this case study, we present the workflow that enables the derivation of p c -S w -a wn relationships of sandstone core-samplebased model domains.We describe a feasible and relatively simple numerical model that enables the dynamic fluid-fluid interface tracking on micro-computed tomography (µ-CT)obtained geometries.Initial and boundary conditions as well as spatial and temporal discretization are derived using a novel approach, taking simplified idealized geometries based on REV-scale pore geometry statistics into account.
Besides using the resulting p c -S w -a wn relationship as input for the model given by Tatomir et al. (2015), this work contributes to designing the KIS tracer laboratory experiments described in Tatomir et al. (2015).Furthermore, the derived relationship aims at improving the characterization of the Heletz sandstone reservoir (Edlmann et al., 2016;Niemi et al., 2016;Tatomir et al., 2016).

Heletz georeservoir geology and mineralogy
The sandstone sample on which this case study is based originates from Heletz, Israel, which is a scientifically motivated deep saline CO 2 storage pilot site (Niemi et al., 2016).This CO 2 injection site extends to depths of 1300 to 1500 m over an area of 22 km 2 .The reservoir itself consists of sandstone layers overlaid by an impermeable cap rock consisting of shale and marl.A detailed description of the site geology is given in Tatomir et al. (2016) and Niemi et al. (2016).With regard to the sample mineralogy, a thin section analysis was carried out to characterize the solid matrix, which was found to consist of Quartz grains interconnected by a Calcite and Clay cementation matrix.X-ray diffraction confirmed these results as described in Tatomir et al. (2016).Sorting and rounding as in Tucker (2011) classified the sample to be very well sorted and rounded with low to high sphericity respectively.The pore space can analogously be defined as homogenously distributed, creating a network of pores of mainly similar sizes.Further information is given in the poresize distribution analysis described later.A photograph of the sample as well as a thin section and scanning electron microscopy (SEM) images are given in Fig. 1.

Mathematical model
All modelling in this study was carried out using COM-SOL Multiphysics (Version 4.3b  (NSEs) as governing equations.Defining a Newtonian fluid of constant density and laminar flow enables us to solve the governing equations in the form of momentum balance (Eq. 1) and continuity (Eq.2) as follows (Cengel and Cimbala, 2013): where ∇ is the nabla operator and t, p, ρ, µ, u, g, T and F st denote time, pressure, density, dynamic viscosity, velocity, gravity, matrix transpose superscript and a surface-tension force component respectively.The latter was added in order to include pore-scale-dependent capillary effects.For simplification, gravity was neglected.Sharp fluid-fluid interface tracking was realized by implementing the phase-field method as described amongst others in Kim (2012).Here, interfaces are defined as regions of rapid change in an auxiliary parameter P or phase field (Popatenko, 2007).Phase change occurs at −1 < P < 1.
The change of phase component concentration over time can be described as the product of diffusion coefficient and chemical potential known as the Cahn-Hilliard equation (Cahn and Hilliard, 1958).In the present model, this equation is used in the form of Eq. (3) using the terms mobility γ , predefined interface thickness ε, mixing energy density λ and a phase field help variable comprising the chemical potential.The latter is calculated according to Eq. (4).
The typical interface thickness parameter is ε = h c /2, whereas h c is the characteristic spatial mesh element size at interface position.Calculating the surface tension coefficient σ = 2 √ 2 3 λ ε in the form of energy per area (Woodruff, 1973) derives the mixing energy density λ = 3σ ε 2 √ 2 . The mobility is solved mathematically using γ = χ ε 2 with the manually chosen mobility tuning parameter χ (here: χ = 1).Volume fractions of fluid 1 (f1) and fluid 2 (f2) can then be calculated according to V f1 = 1−P 2 and V f2 = 1+P 2 .Prior to computation, ρ and µ of both fluids have to be defined.In the phase transition zone, ρ mix and µ mix are calculated according to Eq.( 5) and Eq.( 6). (5) 4 Model validation The model described in this study was validated by reproduction of three analytical solutions.First validation refers to single-phase flow and second and third validation at 2-phase flow and fluid-fluid interface tracking.
The NSE model is validated by simulating plane Poiseuille flow.This type of problem describes pressure-driven steady, incompressible laminar viscous flow through a twodimensional channel (Rogers, 1992).Within the channel, a parabolic velocity profile forms and the maximum velocity (u max ) is in the channel centre.This can be calculated using Eq. ( 7) derived from the NSEs (Pillai and Ramakrishnan, 2006): Here, h, mu and ∂p ∂x denote the channel height, dynamic viscosity and change of pressure over channel length respectively.More detailed descriptions about plane Poiseuille flow are given amongst others in Chughtai et al. (2009) and Shahriari et al. (2013).Comparing the results of the analytical solution with the corresponding u max obtained from modelling (based on a channel geometry of length and height of approx.x = 146 and h = 48 µm, whereas the latter represents the mean pore diameter of the complex µ-CT based geometry) gave a deviation of u max = 0.05 %.Thus, the computation of plane Poiseuille flow can be seen as precise.Here it has to be stated that the channel length was extended on both sides by a correction length x corr = 1/2 × x in order to neglect boundary in-and outflow effects.These extensions were neglected for the validation analysis.A visualization of the model set-up as well as model results in the form of isolines of velocities is given in Fig. 2.
For the validation of the phase-field method regarding the emulation of the sharp fluid-fluid interface, the Rayleigh-Taylor instability has been numerically designed and modelled.In principle, the Rayleigh-Taylor instability is described as the gravity-driven unstable displacement of a light fluid f2 with a heavy fluid f1 (Lord Rayleigh and Strutt, 1883;Taylor, 1950;Mikaelian, 2014).Initial separation of both fluids is given by a cosine-shaped interface of the amplitude η 0 .Over time, f1 breaks into f2 causing a characteristic interface amplitude change ∂η ∂t .With regard to the characteristic  et Al., 1978;Kordilla, 2014), the interface amplitude change can be described mathematically as follows (Read, 1984;Youngs, 1984): A, g and t are Atwood number (dimensionless density ratio: ), gravitational acceleration and time respectively, while α q is a growth rate coefficient that should range within 0.01 and 0.08 (Glimm et al., 2001).Modelling was carried out within a geometry, based on the auxiliary parameter d p = 0.048 m describing a rectangular shape of width w = d p and height h = 2d p .The initial cosine-shaped interface was defined as a function y ) and f2 as supercritical CO 2 (ρ CO 2 = 700 kg m −3 ).Boundary conditions were set to no slip and modelling was solved transiently over a time interval t = 0.75 s using a direct MUMPS (MUltifrontal Massively Paral-  lel sparse direct Solver).Schematic model set-up and model outcome are visualized in Fig. 3.
Regarding the results, α q = 0.068 was found to be a good fit.Comparing the modelling result with the corresponding analytical solution gives a R 2 = 0.99.The associated plot is given in Fig. 4. In conclusion, the model is aligned with the analytical solution to a very high degree.
Further model validation with regard to NSEs and phasefield method in the form of a two-dimensional capillary rise within a tube was carried out.Here, fluid rises within a capillary, contracting air due to imbalanced intermolecular adhesive forces between liquid and solid as well as intermolecular cohesive forces within the liquid.The corresponding analytical solution derived for the fluid column height h is given in the Young-Laplace equation (Young, 1805;Laplace, 1805;Norde, 2003) as follows: where γ , θ, ρ, g and r denote surface tension, contact angle, fluid density, gravitational acceleration and capillary radius respectively.Modelling was carried out on a geometry consisting of a rectangular base with a height and width of 2 × 10 −3 m 4 × 10 −3 m and an overlaying capillary extending to a height of 40 × 10 −3 m and a width of 2 × 10 −3 m.Pressure-driven replacement of air of ρ air = 1.2 kg m −3 by H 2 O (ρ H 2 O as above) was simulated, where a capillary pressure of p c = 145 Pa was defined on one lateral side of the rectangular base (inlet) below the capillary and a pressure of zero at the capillary top.Remaining boundaries were set to no slip.The model was solved transiently over a total time of 2 s using the iterative generalized minimal residual method (GMRES).Model set-up and capillary rise at pressure equilibrium (t = 2 s) are illustrated in Fig. 5.
Comparing the model results in the form of height of water column in the capillary with the corresponding solution of Eq. ( 9) gave a deviation of h = 4.67 %, which can be seen as sufficiently accurate.The mathematical model used within this study can be seen as successfully validated.

Numerical workflow
The numerical workflow steps (e.g. generation of the µ-CT based digital sample approximation, statistical analysis of the pore space properties, the 2-phase flow model set-up) are described in the following section.
In order to test and modify the FEM-based model described in Sect. 3 on the sample-specific pore geometries (e.g.transitions from big pores to small bottlenecks), statistical parameters of the µ-CT-based digital sample approximation had to be determined.A porosity-based REV-analysis was carried out in order to derive the domain size for a poresize distribution (PSD) analysis.Subsequent implementation of a PSD-analysis gave statistical information about average, maximum and minimum pore diameters.This statistical information was used to design idealized model domains and to extract subgeometries from the digital sample approximation.As a next step, the model was implemented stepwise, first on simple geometries (idealized single pore, idealized pore network) and later on complex geometries (µ-CT-based pore, µ-CT-based pore network).This enabled us to test the model performance with regard to very steep pressure gradients (and critical flow) building up on small bottlenecks.If the equation system could not be solved for a given geometry, modifications in the model were done.These modifications include choice of temporal discretization, initial and boundary conditions and, in extreme cases, simplifications in the mathematical model (e.g.neglecting the contact angle by changing apparent slip to no-slip boundary conditions).This was done until a model and model set-up for the successful implementation on complex µ-CT-based geometries were derived.Figure 6  The three-dimensional imaging of the sandstone has been performed with a nanotom S 180 µ-CT (180 kV, 500 mA) device, located at the Leibniz Institute for Applied Geophysics (Department 5, Petrophysics and Borehole Geophysics) in Hanover.The nanotom is a compact CT system for pore-scale imaging purposes, i.e. for high-resolution imaging within the micrometre (typically 1-2 µm) to submicrometre range (about 700 nm for very small samples), featuring high image sharpness due to a significantly reduced penumbra effect (Brunke et al., 2008).The image data were processed with the AVIZO Fire software suite (www.fei.com/software/avizo3d/),including scan artefact reduction and image filtering.Due to the low image noise and the fact that only single-phase segmentation of the pore system was needed, segmentation was performed by the fast and robust "automatic threshold selection method" as described by Otsu (1979).Afterwards, the pore space was transferred into a mixed hexagonal-tetrahedral grid by using the extended volume marching cubes algorithm as implemented within the software ScanIP and as described by Young (2008).This step included an optimization of the segmented grid, which significantly reduced the number of grid cells (almost by a factor of 10) and hence computational time and power, without "destroying" topological information of the grid-surface.The final mesh has been exported as an ASCII-grid ( * .mphtxt),which has been imported to COM-SOL for modelling.

REV-analysis
The interfacial area between the two fluid phases for the macro-scale models is represented as volume-averaged.Therefore, the representative elementary volume (REV) with regard to effective interconnected porosity was identified and extracted from digital sample approximations.Generally, a REV-analysis identifies the minimum porous-media volume with similar mechanical and hydraulic properties of the total porous-media volume (Singhal and Gupta, 2010).Further information about the REV can be found in Bear (1988).A detailed description of the workflow and results specifically for Heletz sandstone cores is given in Tatomir et al. (2016).
Results indicate a REV range of a subvolume edge length bandwidth ranging from approx.120 to 280 voxel.Within this range, fluctuations of porosity are minimal.Subvolumes bigger than 280 voxel edge length show an increase in porosity and can be identified as domain where macroscopic effects are dominant.Subvolumes smaller than 120 voxel can be further distinguished in a transition phase ranging from 50 to 120 voxel, where fluctuations in porosity are relatively low and in domains smaller than 50 voxel, where microscopic effects are dominant and fluctuations are high.

Pore-size distribution analysis
Based on the µ-CT images, a series of operations combining pore identification, binarization, separation and statistical analysis regarding geometric extent of features was implemented on an arbitrarily extracted REV-scale subvolume of the digital sample approximation using the software Avizo Fire (Version 8.1).Pore identification and characterization was based on their specific greyscale value.Further separation was done according to the watershed transformation technique as described amongst others in Sheppard et al. (2005) and Soille (2003).Finally, all relevant data were determined by applying the special "label analysis" operation.The results of the PSD are given in the form of a histogram (Fig. 7) and Table 3, where the latter summarizes most important statistical information, such as minimum, maximum and mean pore diameters.A detailed description of the workflow is given in Tatomir et al. (2016).
It can be observed that the majority of pores are represented by diameters up to 90 µm.In fact, approximately 95 % of pores are below a threshold of the doubled mean pore diameter.This value was considered instead of the maximum pore diameter when designing the idealized geometries due to optimized geometry meshing.It has to be stated that remaining bigger pores may be few; however these pores play a more significant role for fluid flow representing preferential

Model domains and initial model set-up
Model domains were called idealized pore, idealized pore network, µ-CT-based pore and µ-CT-based pore network.A detailed description of each domain is given in the following section.
The idealized pore measures the mean pore diameter (d pore ).A geometric bottleneck with a d throat = 1/2 × d pore was positioned in the middle of the domain in order to test pore throat in-and outflow effects.Domain in-and outlet were scaled to d in,out = 1/4 × d pore .
The idealized pore network domain led from a pore with a diameter d pore simultaneously into two pores, each represented by the minimum pore diameter (d min ) and the doubled mean pore diameter (2 × d pore ).The latter is approximated as a maximum pore diameter for the 95th percentile of pore diameters in the digital sample approximation.Pressuregradient-driven flow was conducted z axially for both idealized geometries.Flow in-and outlet were positioned at z in = 0 µm and z out = z max and the initial interface was lo- Density ρ (kg m −3 ) 973.9 470 Dynamic viscosity µ (Pa s −1 ) 3.59e − 4 0.55e − 4 cated at z = 50 and z = 10 µm for the idealized pore and the idealized pore network.
The µ-CT-based pore geometry was explicitly chosen due to its pore diameter which approximates d pore with d ≈ 40 µm.Furthermore, it was chosen due to the presence of a dead-end pore in order to observe trapping effects.The µ-CT-based pore network geometry was chosen due to its bandwidth of pore diameters 25 µm d 220 µm approximating the corresponding band width of the REV-scale µ-CT geometry.Again, the presence of dead-end pores was considered.Pressure-driven flow was conducted y and z axially for the µ-CT-based pore geometry and the µ-CT-based pore network geometry respectively.Domain in-and outlet were set analogously to the idealized geometries and the initial interfaces were positioned at z = 0.1 × z max and z = 0.25 × z max .Regarding the domain pore sizes, both µ-CT-based geometries directly relate to geometric and later flow characteristics that are present in the REV-scale geometry.All geometries at initial condition are shown in Fig. 8.
The fluid properties for supercritical CO 2 and water are given in Table 1.
The temperature is set to 333.15 K, which is the measured temperature in the Heletz reservoir (Niemi et al., 2016).Analogous to typical on-site injection, a one-axial pressure drop was chosen to induce fluid flow, whereas the overall pressure gradient was defined as capillary pressure (p inlet − p outlet = p c ). Domain in-and outlet were defined as the constant pressure head (Dirichlet) boundary condition (p = const) and viscous stress was neglected [µ(∇u + (∇u) T ) = 0].The model was implemented on a range of capillary pressures 1 MPa < p c < 15 MPa with a p = 2 MPa, whereas the maximum p c approximates a field injection shut-in pressure.In order to take wetting and non-wetting-fluid effects into account, a contact angle θ = 30 • = π/6 rad was introduced.This value for a supercritical CO 2 -H 2 O-SiO 2 system at p c = 10 MPa originates from a study conducted by Espinoza and Santamarina (2010).Here it has to be stated that effects due to a change in θ are of importance when quantifying the interfacial area (a wn ) on a pore scale.A change in contact angle results in a change of interfacial curvature (Jiang and Zeng, 2013).An increase in θ leads to a bigger curvature, which leads to an increase in a wn and vice versa.The computation of θ could be realized by defining an apparent slip boundary condition on all remaining boundaries.This enforces zero velocity at the wall.The product of velocity u and interface normal (using the previously introduced Such an apparent slip introduces a larger velocity difference at the boundary while the fluid speed seems (but in fact does not) to extrapolate to zero below the surface (Granick et al., 2004;Forster and Smith, 2010).In this model, an apparent slip is enabled when adding the boundary condition "wetted wall".This boundary condition adds a frictional boundary force F fr .
where β is slip length and µ viscosity.The following boundary integral consists of the surface tension coefficient σ and the delta function The test function test (u) enables us to neglect (Eq.12) when applying a no-slip boundary condition.Subsequently, this set-up enables θ to be defined.

Final model set-up
For the model based on idealized geometries, an appropriate time-step size was found to be t = 5 × 10 −8 s.Total time intervals were set to t max = 1 × 10 −4 and t max = 1 × 10 −5 s for the idealized single pore and idealized pore network.
Implementing the initial model and model set-up on the µ-CT-based geometries was only possible after modifications.The wetted wall boundary condition had to be changed to a no-slip boundary condition in order to avoid imbibition effects.The contact angle was neglected and thus not considered for later results.The total time interval had to be decreased when applying bigger capillary pressures.In case of the pore benchmark geometry, the initial time interval of t max = 5 × 10 −5 s was reduced to t max = 1 × 10 −5 s for each p c ≥ 5 × 10 6 Pa.In the case of the µ-CT-based pore network geometry, the time interval choice was more variable.It is given in Table 2. Analogous to above, time-step size was kept at t = 5 × 10 −8 s.In order to minimize computation time, the absolute error tolerance of the solver was changed from 0.001 to 0.01.While the µ-CT-based pore geometry modelling could be successfully executed on these settings, model performance on the µ-CT-based pore network geometry had to be stabilized by introducing a Dirichlet-type constant outlet pressure boundary condition.This setting enabled us to solve the equation system and avoid otherwise appearing singularities.Choice of this pressure was found to depend on the inlet capillary pressure.Complete displacement of H 2 O was not possible in the µ-CT-based pore network geometry model, due to trapping effects.Resulting minimum wetting phase saturations S w,min for each p c are given in Table 2.

Mathematical derivation of the p c -S w -a wn Relationship
Since the pressure drop within the domain was kept constant for each simulation, the calculation of p c was realized by subtraction of outlet from inlet pressure.
The wetting phase saturation S w was determined by dividing the volume integral of wetting phase filled domain cells B w by the total domain volume (volume integration V 1dV ).This is mathematically expressed as follows: In order to quantify the specific interfacial area, an isosurface data set was generated.This data set represented the 0.5 concentration of each of the non-wetting and wetting phases.The corresponding auxiliary phase field parameter P is zero.All cells represented by this isosurface are denoted with B P =0 .A subsequent surface integration of this isodata set normalized to the total domain volume gave a wn as follows:

Dynamic evolution of CO 2 -water interface
An exemplary time series of phase and interface propagation within the µ-CT-based pore network geometry is given in Fig. 10.The same is given for the µ-CT-based single pore geometry as Fig. 9. Complete H 2 O replacement (S w = 0) could not be achieved in any of the given domains, which is partly due to trapping and bubble-forming effects.A more detailed description about such effects is given in Saeedi (2012).
Residual water in dead-end pores and bubbles can be observed in the last time step visualized in both Figs. 9 and 10.
Figure 11 visualizes the change of S w over time for µ-CT-based pore network geometry over all applied capillary pressures.The figure is given in log-log axis scale and the initial wetting phase saturation for all scenarios (S w = 0.93) at time step zero is neglected.
It can be clearly observed that an increased capillary pressure results in a significant decrease of wetting phase saturation within the domain.However, at pressures p c > 10.70 Pa, Given p c -S w -a wn relationships show that the maximum observed interfacial area correlates to an increase in capillary pressure.Peaks in a wn distribute for all p c approximately at a range 0.7 > S w > 0.35.It can be stated that both relationships show the typical bell-shaped form, which is widely postulated in the literature (Joekar-Niasar et al., 2010;Karadimitriou and Hassanizadeh, 2012;Karadimitriou et al., 2013).

Discussion
A case study on Heletz sandstone was carried out to present and test a workflow for determining the constitutive capillary pressure-saturation interfacial area relationship based on µ-CT obtained geometries.Based on the digital sample approximation specific geostatistics, modelling domains were designed and extracted.Modelling on these domains enabled to test the model and model set-up on geometry-specific steep pressure gradients (at small bottlenecks) causing critical flow and problems solving the equation system.To avoid the latter, the presented workflow can be seen as a systematic procedure to successively modify the implemented 2-phase flow NSE model towards a successful model implementation for such a type of modelling problem.
The design of the idealized and extracted µ-CT-based geometries can be improved by taking the whole PSD (including pores representing e.g.quantiles) into account and not only average, maximum and minimum pore sizes.For the modelling problem within this case study, a two immiscible fluid phase system consisting of water and CO 2 is considered.For practical field applications, prior to injection the Heletz aquifer is filled with brine, which has a higher density than pure water and a different viscosity.Thus, the model does not represent field conditions.However, envisioned laboratory experiments (e.g.Tatomir et al., 2015) require the determination of the p c -S w -a wn relationships.Therefore, the main objective of this case study was to develop a feasible workflow for obtaining p c -S w -a wn relationships from µ-CTbased geometries of the Heletz sandstone.
The scanned sandstone sub-core samples do not represent in situ reservoir conditions, such as the pore space compaction occurring at reservoir depth (below 1600 m), which implies reduced rock permeability and porosity.The digital sample approximations used for later modelling represent ex situ rock characteristics.Comparing permeability results from unpressurized and pressurized core samples, Tatomir et al. (2016) shows an approximate overestimation of permeability by a factor of 1.5 (on the same investigated core samples).However, when compared with the permeability results of Hingerl et al. (2016), who used a 9 MPa pressurization, the overestimation reaches a factor of 4. Furthermore, the comparison of capillary pressure-saturation relationships obtained from µ-CT modelling and from laboratory results of mercury porosimetry showed important differences (Tatomir et al., 2016).
Furthermore, the accuracy of mesh generation via µ-CT scan must be seen as a function of scanning resolution.Micro-scale features below scanning resolution are lost and subsequent effects of these features on flow characteristics cannot be taken into account.From the mineralogical investigations of the Heletz sandstone the presence of a clay coating on quartz grains was observed.Dry clay minerals are known to swell when in contact with a fluid-like water.This process undergoes with a further change in pore space characteristics, which was neglected for simplification.
The composed numerical model could be successfully validated.NSEs were simplified by defining flow incompressibility.Pressure-driven flow at steep pressure gradients in highly irregular geometries lead to the formation of an extremely heterogeneous flow velocity field.Narrow geometric features such as bottlenecks result in extremely steep pressure gradients and high-flow velocities.Velocities were often found to exceed the fluid-specific speed of sound c.Calculations based on supercritical CO 2 in fully CO 2 -saturated locations showed a maximum exceedance of v ≈ 8×c (MACH ≈ 8).Since incompressibility can be assumed when v < c (Wesseling, 1991), this threshold is considerably exceeded and critical flow must be implied.Hence, neglecting compressibility leads to model inaccuracies.The same applies for the model simplification defining flow to be laminar.At narrow geometric features, calculated Reynolds numbers exceed the threshold for turbulent flow by an approximate factor of 6.7 (Re max ≈ 20 190).Thus, neglecting turbulent flow leads to further model inaccuracies.A further source for inaccuracies can be seen in defining fluid-fluid immiscibility.As described by Tassaing et al. (2004), conditions of water dissolved in CO 2 under modelling conditions would range up to several mol per litre stating miscibility.Modifications on the mathematical model and model set-up carried out during the workflow, such as the negligence of the contact angle were done in order to be able to solve the equation system.The introduction of a judiciously chosen background pressure enabled to avoid extremely steep pressure gradients at narrow geometric features.At this point, modifications and simplifications must be seen as inevitable under the given circumstances.Simplifications including the negligence of compressibility, miscibility and turbulent flow were done in order to avoid a considerable increase in model complexity.Avoiding these simplifications would have led to excessive computational effort and long computation times.Regarding this and the fact that this study aimed at presenting a feasible, easily applicable tool for the tracking and quantification of the dynamic fluid-fluid interface for a very challenging modelling problem with high (geometric) complexity, such simplifications must be seen as a necessity.

Summary and conclusions
In this study we have developed a numerical workflow for the determination of capillary pressure-saturation fluid-fluid interfacial area on highly complex µ-CT obtained geometries.The presented workflow consists of a successive number of model runs on geometries with increasing degree of complexity.
The p c -S w -a wn relationship is obtained from modelling idealized geometries with average, maximum and minimum pore size and estimated pore throats resulting from PSD. Pore-throats are important as they control the entrance of the non-wetting phase in the fully wetting phase saturated domain.This is defined at macro-scale (REV-scale) as the entry pressure.
We have not conducted any comparison of with the real REV geometry because the geometry was too complex to resolve with the 2-phase NSE model.However, by preserving the average properties of the medium and the average porethroat diameters, we attempted to obtain the same kind of accuracy.Future work will concentrate in simulating larger domains on bigger clusters and determine the accuracy of the technique.Also, no comparison with a pore network model was conducted, but it is known that the idealizations imposed by such pore network models also result in incorrect estimation of the pc-S (e.g.Tatomir et al., 2016;Leu et al., 2014;Wildenschild et al., 2005;Wildenschild and Sheppard, 2013).
These results contribute to improving the characterization of the Heletz reservoir where CO 2 injection operations are planned within the next year.

Figure 1 .
Figure 1.Top left: photograph of a well core sample.Top right: thin section image of the sample confirming the well sorted and rounded character.Labels mark the main mineral components quartz, calcite, clay and the pore space.Bottom left: SEM image showing a plane view on a mineral conglomerate of quartz grains interconnected by clayey and carbonatic cement.Bottom right: SEM image highlighting the partially significant amount of cementation (modified from Tatomir et al., 2016).

Figure 2 .
Figure 2. Model set-up for the model validation using plane Poiseuille flow.Coloured lines indicate isolines of velocity, whereas blue and red colours define lower and higher velocities respectively.The subgeometry used for the analysis is marked by the black box.
∂η ∂t , Mikaelian (2014) defines three stages of the Rayleigh-Taylor instability in the form of linear, nonlinear and turbulent regimes.At the transition of nonlinear to turbulent regime (later time steps at η(wavelength λ) = λ 2π , Anuchina

Figure 3 .
Figure 3. Schematic model set-up and phase propagation over time of the model validation using Rayleigh-Taylor instabilities.All Rayleigh-Taylor instability stages in the form of linear, nonlinear and turbulent regime are visualized.Blue and red colours represent heavier (ρ heavy here: H 2 O) and lighter fluid (ρ light here: CO 2 ).For comparability, the cosine-shaped initial interface is given continuously as a black line, while the propagating interface is illustrated in grey.

Figure 4 .
Figure 4. Interface amplitude change η of Rayleigh-Taylor instability over time.The η COMSOL represents the COMSOL modelling result as a grey line while η Read represents the analytical solution for the late nonlinear to early turbulent RT regime as a dashed line.The bold dashed line visualizes the representative part of the analytical solution.

Figure 5 .
Figure 5. Model set-up for the model validation via Young-Laplace capillary rise at pressure equilibrium.Blue and red colours define water and air respectively.Height of water column is h at a capillary radius r and an inlet capillary pressure p c .Orientation axes values are given in metres.

Figure 6 .
Figure 6.Schematic overview of the workflow for deriving the final 2-phase flow model and model set-up for complex µ-CT-based model domains.
illustrates the workflow schematically.6 Generation of the µ-CT-based digital sample approximation, geostatistical analyses, model domain description and initial model set-up 6.1 Domain construction and mesh generation via µ-CT imaging

Figure 7 .
Figure 7. Result of the pore-size distribution for an REV-scale subvolume arbitrarily placed in the sample approximation as number of pores vs. pore diameter, whereas pore diameters were unified in 30 µm intervals.Due to representation, this figure excludes the pore representing the maximum pore diameter d Pore, max , which exceeds the pore diameter range significantly.

Figure 8 .
Figure 8. Heletz sandstone digital sample approximation based on idealized and µ-CT-based model domains used within this study.All domains are visualized at initial condition.Blue colours represent water, red colours CO 2 .The orientation grid is given in µm scale for the idealized pore network domain and for all other domains in millimetres.

Figure 9 .
Figure 9. Image series showing the propagation of supercritical CO 2 displacing H 2 O (above) as well as the corresponding propagation of the fluid-fluid interface (below) at p c = 15e6 Pa.The series is given for selected time steps within the µ-CT-based single pore geometry.Blue, red and green colours represent H 2 O, CO 2 phase and the fluid-fluid interface respectively.Grid values are given in millimetres.

Figure 10 .
Figure 10.Image series showing the propagation of supercritical CO 2 displacing H 2 O (above) as well as the corresponding propagation of the fluid-fluid interface (below) at p c = 15e6 Pa.The series is given for selected time steps within the µ-CT-based pore network geometry.Blue, red and green colours represent H 2 O, CO 2 phase and the fluid-fluid interface respectively.Grid values are given in millimetres.

Figure 11 .
Figure 11.Plot of S w over t for each p c .

Figure 12 .
Figure 12. 3-D bar plots of the p c -S w -a wn relationship obtained from the µ-CT-based pore geometry (above) and the µ-CT-based pore network geometry (below).Note the different scaling in a wn .

Table 1 .
Fluid properties representing supercritical CO 2 and H 2 O in the models.

Table 2 .
Outlet pressure p out , maximum model time interval t max and minimum-reached wetting phase saturation S w, min for each corresponding capillary pressure p c for the µ-CT-based pore network geometry model.

Table 3 .
Results of the pore-size distribution highlighting values for number of pores as well as mean, minimum and maximum pore diameter.