1 Introduction

Multiphase flow in porous media is an important scientific phenomenon that occurs in many natural and manufactured settings (Blunt 2017; Sahimi 2011; Pruess and Garcia 2002). Two-phase flow through porous media occurs during hydrocarbon recovery in reservoirs and carbon dioxide storage, for instance (Blunt 2017; Pruess and Garcia 2002; Dullien 2012). Traditionally experiments on two-phase flow through porous media have been conducted at the core scale (several cm) (Culligan et al. 2006; Bourbiaux and Kalaydjian 1990). Capillary pressure (\(P_{c}\)) and relative permeability (\(k_{r}\)) are determined to quantify multiphase flow properties (Blunt 2017; Sahimi 2011; Akbarabadi and Piri 2013; Al Mansoori et al. 2010; Pentland et al. 2010).

Advances in high-resolution three-dimensional x-ray imaging of porous media have made it possible to extract the geometry of the pore structure and to observe the fluids within the pore space (Hazlett 1995; Ketcham and Carlson 2001; Prodanović et al. 2007; Silin et al. 2011). Moreover, the use of synchrotron radiation allows images to be acquired every few 10s of seconds, allowing the micron-scale (or pore-scale) pattern of invasion to be visualised (Singh et al. 2017; Wildenschild and Sheppard 2013; Cnudde and Boone 2013; Bultreys et al. 2019; Berg et al. 2013; Andrew et al. 2015).

X-ray imaging can capture the structure of the pore space with micron resolution for input into numerical models (Blunt 2001; Blunt et al. 2013). Different numerical methods to explore pore-scale physics have been proposed, which can be classified into two groups; direct numerical simulation (DNS) and pore network modelling. In DNS methods, the simulation grid can be obtained directly from an image. DNS includes lattice-Boltzmann methods (Boek and Venturoli 2010; Ramstad et al. 2012; Akai et al. 2019; McClure et al. 2018), finite volume simulation (Ferrari and Lunati 2013; Raeini et al. 2014), and the level set approach (Prodanović and Bryant 2006). In contrast, pore network modelling provides an upscaled representation of the pore-scale physics (Blunt 2017). The simulation domain, of interconnected pores and throats, has to be extracted from an image of the pore space (Blunt 2017; Raeini et al. 2017).

One challenge of any modelling approach is to match the pore-by-pore displacement sequence (Bultreys et al. 2019). If a model is tuned or calibrated to reproduce experimentally measured continuum scale flow properties, for example, relative permeability and capillary pressure-saturation functions, the match is likely to be nonunique which means that its predictive power outside the range of the experiments used for comparison is limited (Aghaei and Piri 2015; Bondino et al. 2013; Masalmeh et al. 2015; Bultreys et al. 2019).

The first studies which compared fluid distributions with experiment were confined to quasi-two-dimensional micro-fluidics studies, as opposed to displacement in three-dimensional rock samples (Joekar-Niasar and Hassanizadeh 2012; Zhao et al. 2019). Qualitative comparisons between experimental fluid distributions and model predictions have been shown (Armstrong et al. 2015; Aghaei and Piri 2015; Armstrong et al. 2016). More recently though, comparisons of imaging experiments with pore-scale models have been the subject of several studies. Bultreys et al. (2018) mapped the measured fluid distribution of steady-state micro-CT experiments to the corresponding pore networks and performed comparisons of model and experiment in water-wet Bentheimer sandstone. In another study by Øren et al. (2019) pore-by-pore comparison between the distribution of fluid measured in a micro-CT experiment of CO\(_2\) storage in a water-wet sample showed good agreement using a sophisticated model to estimate capillary entry pressures in pores. Later, Raeini et al. (2019) quantified the mismatch between measured fluid distributions and the results of their model. They showed that even in replicate experiments on the same sample a volume-weighted mismatch of 10–20% in pore occupancy is seen, which sets a benchmark for assessing the accuracy of model predictions (Andrew et al. 2013, 2014). These studies performed comparisons either at the end of a displacement when the fluids were at rest, or during steady-state flow, and hence essentially predicted fluid configurations at fixed points, rather than the correct time sequence of pore-scale filling. In another study, Bultreys et al. (2019) performed a pore-by-pore analysis to examine the performance of a network model against results from time-resolved micro-CT data sets. The experimental data in their study were two publicly available synchrotron studies of waterflooding in strongly water-wet Ketton limestone (Singh et al. 2017) and Gildehauser sandstone (Rücker et al. 2015). One major uncertainty in predictive models, namely the local distribution of wettability—the contact angle on a pore-by-pore basis—was not addressed. Water-wet datasets were considered where it was known that the contact angles everywhere were less than 90\(^{\circ }\); however, most reservoir rocks have undergone a wettability alteration on contact with crude oil with a range of local contact angles in the pore space (Kovscek et al. 1993; AlRatrout et al. 2018; Buckley 2001).

In a previous study (Foroughi et al. 2020), we proposed a framework to adjust local contact angles to match the fluid configurations imaged in steady-state displacement experiments. The initial estimate of contact angle was based on the average angle which obeyed energy balance for quasi-static displacement with no viscous dissipation (Blunt et al. 2019). We applied our framework to two steady-state experiments on water-wet and mixed-wet Bentheimer samples conducted by Lin et al. (2018, 2019). We showed that the proposed algorithm could decrease the mismatch error in pore occupancy between the model and experiment to below that seen between repeat experiments. The relative permeability and capillary pressure predicted by the model then matched the experimental results. It was also shown that there was a spatial correlation in the contact angle distribution with a correlation length of approximately four pores, with a tendency for the larger pores to be more oil-wet. Direct simulation using the lattice Boltzmann method has also been used to predict the pore-scale distribution of fluid, compared to experimental measurements. To obtain a good match a distribution of contact angle was needed (Akai et al. 2019; Zhao et al. 2018; Bakhshian and Hosseini 2019).

In this paper we extend previous work in two ways. Firstly, we study carbonates (Scanziani et al. 2020; Alhosani et al. 2020a), including a sample from a producing reservoir, which has a more complex pore structure than the sandstones studied hitherto. These samples are more oil-wet than the mixed-wet sandstone studied previously (Foroughi et al. 2020). Secondly, we apply our calibration algorithm to time-resolved synchrotron datasets which allows a more precise validation of the filling sequence. By doing so, we further constrain the model to increase its predictive power.

In our previous work (Foroughi et al. 2020), we considered the steady-state experiments in which both relative permeability and capillary pressure were measured but only 7 or 8 fluid configurations were imaged corresponding to equilibrium for a specific water fractional flow. In this paper, we analyse synchrotron datasets with 65 and 78 images where the displacement sequence is accurately captured. This presents a much more rigorous test of our calibration algorithm. Furthermore, we will track the interface between the fluid phases and compare with the model predictions to quantify the accuracy of the predicted displacement sequence.

The paper is organised as follows. In Sect. 2, pore network modelling and the micro-CT multiphase flow experiments are described. Then, the workflow to perform pore-by-pore comparisons of micro-CT experimental measurements against the pore network model is presented. This methodology is implemented on two time-resolved x-ray micro-CT experiments in Sects. 3.1 and 3.2 together with comparisons between predicted and measured capillary pressures. In Sect. 3.3, we examine the correlation in contact angle in our calibrated models.

2 Materials and Methods

2.1 Dynamic Experiments

The measurements analysed in this study are from two unsteady-state (constant injection flow rate) multiphase flow experiments. Both cases studied waterflooding, and the results were published recently (Scanziani et al. 2020; Alhosani et al. 2020a). Both samples were placed in contact with crude oil for an extended period of time at high pressure and temperature. This mimics conditions in subsurface hydrocarbon reservoirs: surface-active components of the crude oil adhere to the solid surface changing its wettability. We call these samples oil-wet because—as we show later—the flow behaviour can only be modelled by assuming that almost all the contact angles lie above \(90^{\circ }\). In both experiments, the experimental fluids used to represent the oil and water phases were doped n-decane (15% wt. 1-iododecane) and doped deionised water (20% wt. potassium iodide), respectively. The experimental fluids were doped to increase their x-ray attenuation values, to enhance the image contrast between the phases.

Both experiments were conducted on small cylindrical samples. The first sample was a Ketton limestone rock: the captured image had a diameter of 3.67 mm and a length of 3.34 mm. The other sample was a carbonate reservoir rock, extracted from a large producing hydrocarbon reservoir in the Middle East. The diameter and length of the corresponding captured image were 3.85 mm and 3.62 mm, respectively. The Ketton sample had a helium porosity of 28% of which 14% was from macropores (defined as pores that could be explicitly seen at the resolution of the image). The measured permeability was \(2.8 \times 10^{-12}\) \(\mathrm{m}^2\). The reservoir carbonate had a measured helium porosity of 26%, with macroporosity accounting for 16%. The permeability was \({2.7\times 10^{-13}}\) \(\mathrm{m}^2\) (Alhosani et al. 2020b).

For both cases, the experiments were initiated by flushing oil at a flow rate of 0.1 \(\mathrm{ml/min}\) to replace the crude oil used to alter the wettability. Therefore, for both samples, the initial state was oil (doped n-decane) and a low initial water saturation present in micro-porosity, smaller macropores, and in the corners of the pore space. Temperature and pressure conditions were 50 \(^{\circ }C\) and 8 MPa. Then, the dynamic experiment was started by injecting brine at a very low flow rate. For the Ketton sample, the brine was injected at 0.5 \(\upmu \mathrm{L/min}\) while for the carbonate reservoir sample this value was 0.15 \(\upmu \mathrm{L/min}\). For both experiments, the densities of brine and oil were 1145.0 \(\mathrm{kg\,m}^{-3}\) and 715.2 \(\mathrm{kg\,m}^{-3}\) respectively; the viscosity of brine and oil were measured to be 0.547 \(\mathrm{mPa\,s}\) and 1.088 \(\mathrm{mPa\,s}\), respectively. The interfacial tension between oil and brine was measured as 52.1 \(\mathrm{mN\,m}^{-1}\). The capillary number (\({Ca =\upmu q/\sigma }\), where \(\sigma\) is the oil-brine interfacial tension, \(\upmu\) is the viscosity of brine, and q is the Darcy velocity of brine) was \(3.3\times 10^{-9}\) in the Ketton case and \(2.1\times 10^{-9}\) for the reservoir sample.

As can be seen, we have low capillary numbers for both cases which shows that flow occurred under capillary-dominated conditions. Furthermore, this allowed images to capture the pore-scale displacement with a temporal resolution of 74 s and 70 s for Ketton and the reservoir samples, respectively (this is the time required to acquire a complete tomogram). During water injection, 65 scans were acquired for Ketton and 76 scans for the reservoir rock. For Ketton each scan contains \(1050 \times 1050 \times 953\) voxels with a side length of 3.5 \(\mu m\). For the reservoir sample each scan contains \(1100 \times 1100 \times 1035\) voxels also with a side length of 3.5 \(\mu m\). Details of the segmentation of the images acquired into rock, oil and brine phases can be found in Scanziani et al. (2020) and Alhosani et al. (2020a).

2.2 Pore Network Extraction and Modelling

Although direct numerical methods based on solving the Navier–Stokes equation numerically can be used to simulate two-phase flow at the pore scale, they are limited in terms of the domain size and the simulation time (Bultreys et al. 2019; McClure et al. 2018; Alpak et al. 2018). In this work, we use a pore network model that allows the simulation of two-phase flow on a large number of pores in only a few minutes (Blunt et al. 1992; Joekar-Niasar et al. 2010). The pore network model is an upscaled representation of the pore space and fluid displacement. In this conceptualised model, the pore space is divided into a network of large void space, pores, which are connected through restriction called throats (Dong and Blunt 2009; Raeini et al. 2017).

To determine the network for each of the datasets, the segmented dry scan image of the pore space is run through a maximal ball pore network extraction code (Dong and Blunt 2009; Raeini et al. 2017). To determine the pore and throat volume, it is necessary to partition pore volume between pores and throats. The border between a pore and its connected throat is at distance \(l_i\) from the pore centre which defined as follows:

$$\begin{aligned} l_i=l_t^p (1-\alpha r_t/r_p ), \end{aligned}$$
(1)

where \(l_t^p\) is the Euclidean distance between the pore and throat, and \(r_p\) and \(r_t\) are the pore and throat inscribed radii. \(\alpha\) is the pore-throat segmentation coefficient on previous work \(\alpha =0.6\) is optimum (Dong and Blunt 2009). After determining the border of pore and throat, each voxel of the pore space is assigned to an element: that is a pore or a throat. Therefore, we can determine pore and throat volumes, \(V_i\).

For both datasets reported in Sect. 2.1, based on the available high-quality segmented images from dry scans, pore networks were extracted. For the Ketton sample, the extracted network consists of 12,319 throats and 5905 pores, while for the carbonate reservoir rock sample, the extracted network consists of 137,730 throats and 101,202 pores. A cross section of the samples and the extracted networks are displayed in Fig. 1. The pore and throat radii distributions are also shown in Fig. 1. As can be seen, the reservoir rock sample is composed of smaller pores and throats compared to the Ketton rock sample, which is why there are more pores and throats in the image of the reservoir sample.

Fig. 1
figure 1

Two-dimensional cross sections of the three-dimensional dry scans of Ketton (a) and the reservoir sample (b). c The extracted network for Ketton, and d for the reservoir sample where throats are shown as lines and pores as balls: the colours provide an indication of throat radius. e Volume-weighted pore and f throat radius distributions for both Ketton (green solid lines) and the carbonate reservoir sample (red dashed lines)

Full size image

After extracting the pore network, flow simulation can be performed. As discussed in Sect. 2.1, in both experiments flow occurs under predominantly capillary-controlled conditions. Therefore, the quasi-static pore network model developed by Valvatne and Blunt (2004) is used here to simulate waterflooding. Filling occurs in order of the entry capillary pressure of the elements, which are calculated semianalytically based on the geometry of the element and the contact angle. Different types of displacement can occur, including snap-off, piston-like advance, and cooperative pore body filling.

Water is allowed to enter both end faces of the network model, through which oil can also escape to account for back-flow in the porous medium outside the field of view. When a phase is not connected to either side of network through the centres of elements or through the fluid layers in the corners of elements, it is trapped. For the initial conditions, the modification proposed by Bultreys et al. (2019) is used in which the fluid occupancy at the start of waterflooding in the experimental data is used to assign the initial conditions of the model. This was done by mapping the experimental fluid occupancy to the model.

The extracted network here is based on the resolved porosity of the rocks from dry scan images. As stated in Sect. 2.1, the measured helium porosities of Ketton and the reservoir rock are 28% and 26%, respectively. In contrast, the resolved porosities of samples are 14.3% and 16%, respectively. We assume that unresolved micro-porosity remains water-saturated in the experiments and does not contribute significantly to the flow, or impact the displacement sequence during waterflooding. The networks are well-connected and fewer than 1% of the pores are isolated. The predicted permeability of the Ketton and reservoir networks is \(2.20 \times 10^{-12}\) \(\mathrm{m}^2\) and \(2.99 \times 10^{-13}\) \(\mathrm{m}^2\) respectively, which are similar to the measured values of \(2.8 \times 10^{-12}\) \(\mathrm{m}^2\) and \(2.7 \times 10^{-13}\) \(\mathrm{m}^2\). For these samples, it appears that indeed unresolved porosity does not significantly affect the flow properties. However, in other cases, such as Estaillades limestone, the unresolved pore space accounts for more than half the total porosity and does have a significant effect on flow and connectivity (Lin et al. 2021). Incorporating micro-porosity into our pore network model will be a topic of future work, based on multiscale models (Jiang et al. 2013; Mehmani et al. 2013; Mehmani and Prodanović 2014; Bultreys et al. 2015; Blunt 2017; Bijeljic et al. 2018).

2.3 Pore-by-Pore Mismatch Analysis

The calculation of the mismatch in occupancy involves the following steps: (1) mapping voxelised micro-CT images of the pore space to the extracted network; (2) determining the mismatch between model and experiment at the network level based on pore occupancy; and (3) quantifying the mismatch. These steps will be discussed in more detail below.

As discussed in previous studies (Bultreys et al. 2018; Raeini et al. 2019; Øren et al. 2019; Bultreys et al. 2019; Foroughi et al. 2020), the extracted network is used as an image analysis tool: for the inscribed sphere centred on each element i, the phase that occupies the majority of voxels in that sphere determines the filling state of that element, \(\phi _i\), where for the oil phase \(\phi _i=0\) and for water \(\phi _i=1\). After mapping, we could determine the average water occupancy, \(\psi\), for each image as follows:

$$\begin{aligned} \psi =\frac{\sum _{i}V_i\phi _i}{\sum _{i}V_i}, \end{aligned}$$
(2)

where the sum is over all elements (pores and throats) and \(V_i\) denotes the volume of element, i. Consider the experimentally measured fluid distribution corresponding to image j which has an average water occupancy \(\psi ^{ex,j}\). To calculate the pore-by-pore discrepancy of the model compared to the experiment, we export the model fluid distribution at an average water occupancy of \(\psi ^{m,j}\) (volume-weighted fraction of water-filled elements), where the difference between these two average occupancies is as small as possible—that is we take the time step in the simulation with the occupancy closest to the experiment. Then the mean absolute discrepancy is defined as follows:

$$\begin{aligned} M_{AD}^j=\frac{\sum \limits _{i} V_i|{\phi _{i}^{m,j}-\phi _{i}^{ex,j}} |}{\sum \limits _{i} V_i}. \end{aligned}$$
(3)

For each element, the absolute difference between the filling state in the model, \(\phi _{i}^{m,j}\), and experiment, \(\phi _{i}^{ex,j}\) when there is a mismatch is 1. The discrepancy is weighted by the volume of the element. A schematic representation is shown in Fig. 2 to illustrate the estimation of the mismatch.

Fig. 2
figure 2

A schematic representation of the mismatch in pore occupancy between model and experiment

Full size image

Throughout the displacement, the following integral measure to capture the discrepancy can be defined:

$$\begin{aligned} A_{M}(\psi _f) = \sum \limits _{j=1}^{n_{s}} \frac{(M_{AD}^j+M_{AD}^{j-1})}{2}(\psi _j-\psi _{j-1}), \end{aligned}$$
(4)

where \(\psi _f\) is average water occupancy at the end of invasion obtained from the last scan, and \(n_s\) is number of scans after the initial condition (\(\psi _f \equiv \psi _{n_s}\)). Note that for initial scan (j=0), the experimental occupancy is read into the model and therefore \(M_{AD}^0=0\).

To test the predictive capability of the model, we compare to a case where elements are filled randomly — a model with no predictive power. It can be shown easily that for a large system, the mean absolute discrepancy of the random model is:

$$\begin{aligned} M_{AD}^R(\psi )=2\psi (1-\psi ). \end{aligned}$$
(5)

As before, the integral of \(M_{AD}^{R}\) gives a measure of discrepancy throughout the invasion:

$$\begin{aligned} A_{R}(\psi ) = \int _{\psi _{i}}^{\psi } M_{AD}^R(\psi ')d\psi ' =\psi ^2-\frac{2}{3}\psi ^3-\psi _i^2+\frac{2}{3}\ \psi _i^3 . \end{aligned}$$
(6)

Therefore, similar to the previous study which was on steady-state experiments (Foroughi et al. 2020), we can define an additional characteristic to quantify discrepancy in dynamic data as follows:

$$\begin{aligned} E=\frac{A_M(\psi _f)}{A_R(\psi _f)}. \end{aligned}$$
(7)

For the case that we have completely correct predictions \(E=0\), while \(E=1\) corresponds to the condition in which model prediction is no better than random assignment.

2.4 Inverse Problem and Method of Solution to Characterise Wettability

The entry capillary pressure is a function of geometry, wettability and flow history. The geometry of the network is extracted from a high-quality dry image and we consider it sufficient to describe the exact geometry of the pore space. We import the initial filling state of the model from the results of the experiment to capture the previous displacement history. However, wettability is an uncertain parameter that requires characterisation.

Blunt et al. (2019) have derived a thermodynamically consistent contact angle which is calculated from differences between successive pore-scale images ignoring viscous dissipation such that all the energy injected into the system is converted to surface energy. The equation is derived when phase 1 (water in this case) displaces phase 2 (oil) in a porous medium. By assuming conservation of Helmholtz free energy between two local equilibrium states, this equation is as follows:

$$\begin{aligned} {\varDelta }a_{1s}\mathrm{cos}(\theta _t)=\kappa \phi \varDelta S_1+{\varDelta }a_{12}, \end{aligned}$$
(8)

where a is the interfacial area per unit volume, \(\phi\) is the porosity, S is the saturation and \(\kappa\) is the total curvature of the fluid–fluid interface. The subscript s denotes the solid, and \(\varDelta\) denotes changes in saturation and area. In this paper we ignore unresolved micro-porosity, assuming that it remains water-saturated throughout the experiment. S is the saturation measured in the resolved macroporosity. With the advent of high-resolution time-resolved three-dimensional X-ray imaging, all the terms in this expression can be measured directly and used to determine the thermodynamic contact angle, \(\theta _t\). It is the angle to use in a model of displacement which, like our network model, assumes capillary-dominated flow with no viscous dissipation. Although this contact angle describes the average wettability state of the system reliably, we know from geometric contact angle studies that contact angle varies spatially (Alhammadi et al. 2017; AlRatrout et al. 2017, 2018).

We have ignored viscous dissipation in this analysis. Its impact has been analysed extensively using lattice Boltzmann simulation, which does account for this effect, by Akai et al. (2020). They showed that ignoring viscous dissipation made a difference of typically less than \(10^{\circ }\) in the predicted contact angle for waterflooding in mixed-wet and oil-wet porous media. As we show later, the thermodynamic contact angle provides an initial guess for our optimisation. Later during the process of optimisation, we adjust the value of contact angle for each pore individually to capture the observed filling sequence. Using the thermodynamic contact angle as an initial guess allows us to reproduce the observed pore occupancy and accurately predicts capillary pressure. We will show that using different initial guesses of contact angle leads to poorer predictions, particularly of capillary pressure.

We define an inverse problem. There is a set of observed fluid distributions (mapped micro-CT images) which we assume is determined by the local distribution of contact angle, which we wish to find. We can calculate the pore-by-pore discrepancy (\(M_{AD}\), Eq. 3) defined in the previous section as our objective function. To minimise this objective function, the contact angle assigned to elements of the network is updated using the following relationship:

$$\begin{aligned} {\theta }_i^{k+1} ={\theta }_i^{k}+(\phi _i^{ex,j+1}-\phi _i^{m,j+1})\mathrm{d}\theta . \end{aligned}$$
(9)

where \(\theta _i\) is the contact angle of element i, d\(\theta\) is an updating step, which we choose to be 0.5\(^{\circ }\). Eq. 9 is used to solve a nonlinear optimisation problem. Readers are referred to our previous study for more details of the algorithm (Foroughi et al. 2020). This algorithm could be classified as an array of gradient descent optimisations in which the direction of the gradient is obtained heuristically by assuming a monotonic relationship between capillary entry pressure (and hence filling state) and contact angle.

The thermodynamic contact angle is used as the initial guess applied to all the network elements. The optimisation is performed step-wise. The contact angles are adjusted to match image j using Eq. 9 until \(M_{AD}\) is less than a specified error criterion. In this work, we choose 0.075 which is the lowest discrepancy observed in repeat experiments (Raeini et al. 2019; Andrew et al. 2013, 2014; Foroughi et al. 2020). Once this threshold is reached, we continue to the next image \(j+1\) and again adjust the contact angles. The maximum number of iterations, \(k^{\rm max},\) is equal to 500, while on average each step takes less than 100 iterations to converge.

3 Results and Discussion

3.1 Ketton

The thermodynamic contact angle computed by Scanziani et al. (2020) was 126°. In the following, firstly we analyse the results of pore-by-pore modelling based on the proposed optimisation algorithm, and then we will validate the predicted capillary pressure. Finally, we predict the relative permeability.

3.1.1 Pore-by-Pore Analysis

To explore the sensitivities of the model, we consider five contact angles as initial guesses in our optimisation. These values are 100°, 111°, 126°, 135° and 150°; 126° is the thermodynamic contact angle. The discrepancy, \(M_{AD}\), Eq. 3 is shown in Fig. 3. Note that when we try to match the pore occupancy of a given time, it has an effect on the earlier pore and throat filling states to some extent. As a result, through matching \(M_{AD}\) for the later time series, the mismatch \(M_{AD}\) earlier may degrade and exceed the error criterion. This effect is stronger when the initial guess is farther from thermodynamic contact angle.

Fig. 3
figure 3

\(M_{AD}\) , Eq. 3, as a function of occupancy comparing the discrepancy between experiment and the model for 5 initial estimates of contact angle (Ketton sample). \(\theta _t\) is the thermodynamic contact angle

Full size image

During the displacement, we can also study how filling proceeds. Since the majority of the pore space initially occupied by oil becomes oil-wet, displacement proceeds approximately as an invasion-percolation like process with the water filling the widest available throat followed by the adjoining pore (Scanziani et al. 2020; Alhosani et al. 2020a; Wilkinson and Willemsen 1983). A throat is considered available if one of its neighbouring pores is occupied by water. In Fig. 4, the available throats throughout the displacement and the ones filled are shown as a function of radius. There is a tendency for the larger throats to be filled preferentially. This is captured accurately by the model. Note, however, that a very small number of smaller throats are invaded; as discussed later, these tend to have slightly lower contact angles, making them more favourable for invasion.

Scanziani et al. (2020) observed that although there was a tendency for water to invade throats with a larger radius, some smaller elements were also filled. In this work, we consider the mismatch based on a volume-weighted discrepancy which tends to disguise the filling of smaller elements; however, in both the experiment and the model we do see some filling of elements of all size. This behaviour can only be accurately captured by allowing a range of contact angle so that both pore geometry and wettability control the filling sequence.

Fig. 4
figure 4

Invaded throats during waterflooding compared to all the available throats for both experiment and model (Ketton sample)

Full size image

Another way to show the accuracy of the model is to plot the distribution of filled elements as a function of inscribed radius. Figure 5 shows the pore and throats filled with water at the end of waterflooding for the optimised case using the thermodynamic contact angle as the initial guess. The model correctly captures the filling of elements of all size, with some discrepancy for the pores and throats of intermediate radii. Overall there is a tendency to fill the larger elements: these are oil-wet with contact angles above 90°.

Fig. 5
figure 5

Pore and throat sizes that are filled with water at the end of waterflooding comparing the model and experiment (Ketton sample)

Full size image

Figure 6 shows the optimised contact angle distribution using the thermodynamic angle, 126° as the initial estimate. To match the data, a distribution is required about this average with some smaller values for elements that fill early in the displacement sequence. However, almost all the values are greater that \(90^{\circ }\) indicating oil-wet conditions overall. Our optimised distribution for pores has a mean contact angle of \(123.9^{\circ }\) and a standard deviation of \(8.4^{\circ }\). Similarly for throats, the mean and standard deviation of contact angle are \(123.4^{\circ }\) and \(7.6^{\circ },\) respectively.

Fig. 6
figure 6

Optimised contact angle for pores and throats for the Ketton sample. The thermodynamic contact angle of 126\(^{\circ }\) was used as the initial guess for both pores and throats

Full size image

3.1.2 Capillary Pressure and Relative Permeability by the Calibrated Model

In the previous section it was shown that our methodology leads to a distribution of contact angle which leads to accurate predictions of pore-by-pore occupancy. However, we need to show that this results in a good prediction of macroscopic properties.

The capillary pressure for the experiment is calculated from the oil-water interfacial curvature and the Young–Laplace equation as reported by Scanziani et al. (2020). Figure 7 shows the predicted and measured capillary pressures. If the thermodynamic contact angle is used, the predictions match the experiment to within the uncertainty in the data. Furthermore, we can make predictions over the whole saturation range. Using other initial estimates of contact angle, while they can give good predictions of pore occupancy, see Fig. 3, with the correct rank ordering of contact angles, they are either all too high or too low, giving erroneous estimates of capillary pressure. Note that the current experiments are unsteady-state experiments, where only a narrow saturation range is explored: there is little displacement after water breakthrough.

Figure 7 also shows the predicted relative permeabilities. Here, since these are not measured, we do not have experimental data for comparison. The predictions show oil-wet type behaviour (Blunt 2017) with a high maximum water relative permeability, a wide saturation range over which the oil relative permeability is low, since oil is confined to layers and the smaller elements, but with a very low residual oil saturation, again thanks to the maintenance of oil connectivity through wetting layers. The water relative permeability shows some finite size effects, since we are modelling a small experimental sample. In the experiments, there is an inevitable trade-off between overall system size and image resolution; we could reduce size effects through the study of a larger sample.

Fig. 7
figure 7

a Capillary pressure based on our optimisation algorithm for five initial estimates of contact angle. The error bars represent the uncertainties in the experimental measurements, b Prediction of relative permeability using the thermodynamic contact angle of 126\(^{\circ }\) as the initial estimate c logarithmic-scale representation of relative permeability

Full size image

3.2 Carbonate Reservoir Rock

The second fast synchrotron data set studied here is the experiment conducted by Alhosani et al. (2020a). This experiment was performed on a carbonate reservoir rock sample under similar conditions to the previous Ketton experiment. This is a more challenging dataset to segment and analyse. In this paper the extraction and smoothing of interfaces was repeated using commercial image analysis software: this resulted in an average thermodynamic contact angle of \(111^{\circ }\). The results of image analysis for 49 consecutive time steps are presented in Table 1 using Eq. 8 between successive images to calculate the contact angle. These results are based on the images captured before water breakthrough. According to these results, the mean value of the thermodynamic contact angle is 111.0\(^{\circ }\) and its standard deviation is 6.2\(^{\circ }\).

Table 1 Data from the oil-wet reservoir carbonate experiment
Full size table

3.2.1 Pore-by-Pore Analysis

Our optimisation algorithm was implemented to minimise \(M_{AD}\) as before. We used the same initial estimates of contact angle as for Ketton: 100\(^{\circ }\), 111\(^{\circ }\), 126\(^{\circ }\), 135\(^{\circ }\) and 150\(^{\circ }\). Note that as discussed in the previous section, 111\(^{\circ }\) is the correct thermodynamic contact angle based on our calculations. The results shown in Fig. 8 indicate that the higher contact angles of 135\(^{\circ }\) and 150\(^{\circ }\) as initial guesses led to poor results for pore-by-pore discrepancy, which means that the optimisation routine could not converge during the displacement. For the initial guess of 111\(^{\circ }\) which is thermodynamic contact angle, the value of E is 0.15 and \(M_{AD}\) at the end of waterflooding is 0.07.

Fig. 8
figure 8

\(M_{AD}\), Eq. 3, as a function of occupancy comparing the discrepancy between experiment and the model for optimised distributions of contact angle based on the three initial estimates shown (carbonate reservoir rock sample)

Full size image

In terms of the filling sequence, Fig. 9 shows that the model accurately captures the size distribution of throats available for filling and those that are actually filled. As seen for Ketton, and consistent with a rock that is, overall, oil-wet, there is a tendency for the larger elements to be preferentially filled with water. Although not evident in Fig. 9, the filling is more strictly in order of size with a stronger tendency for the larger pores to be filled preferentially. Here we do observe an invasion percolation displacement (Alhosani et al. 2020a). This is not because the rock is more oil-wet—indeed the thermodynamic contact angle is lower than for Ketton—but that the pore and throat size distributions, see Fig. 1, mean that size, rather than wettability, has a greater impact on the filling sequence in this heterogeneous reservoir rock, as compared to a more homogeneous quarry limestone (Ketton).

Fig. 9
figure 9

Invaded throats during waterflooding compared to all the available throats for both experiment and model (carbonate reservoir sample)

Full size image

Figure 10 shows the pores and throats filled at the end of waterflooding as a function of radius, comparing the experiment and model with the thermodynamic contact angle as the initial estimate. As with Ketton, the agreement is excellent with discrepancies confined to a few elements of intermediate size.

Fig. 10
figure 10

Pore and throat sizes that are filled with water at the end of waterflooding for the optimised contact angles comparing the model and experiment (carbonate reservoir rock)

Full size image

Figure 11 shows the distribution of contact angles: here we see a distribution around the thermodynamic value. Unlike in Ketton, most values are close to the average, although the pores, in general, do have slightly larger contact angles than the throats. Almost all values are greater than \(90^{\circ }\) indicating oil-wet behaviour. Our optimised model predicts a mean of \(110.3^{\circ }\) and a standard deviation of \(7.3^{\circ }\) for pores, while for throats mean and standard deviation are \(108.9^{\circ }\) and \(5.7^{\circ },\) respectively.

Fig. 11
figure 11

Optimised contact angle for pores and throats for the carbonate reservoir rock sample. The thermodynamic contact angle of 111\(^{\circ }\) was used as the initial guess for both pores and throats

Full size image

3.2.2 Capillary Pressure and Relative Permeability by the Calibrated Model

Figure 12 shows the measured capillary pressure compared to the model predictions for different initial estimates of contact angle. Here, the uncertainty in the measurements is significant. However, only through using the correct average estimate of wettability can the predictions match the data to within experimental error. Also, shown are the predicted relative permeabilities which show characteristic oil-wet behaviour, as expected, with a rapidly rising water relative permeability and a low oil relative permeability at high water saturation, albeit with a low residual saturation. Similar to the Ketton results in Fig. 7c, in Fig. 12c, there is also evidence of finite size effects, with a very low water relative permeability before breakthrough.

Fig. 12
figure 12

a Capillary pressure based on our optimisation algorithm for different initial estimates of contact angle indicated. The error bars represent the uncertainties in the experimental measurements, b Prediction of relative permeability based on optimised distribution using thermodynamic contact angle of 111\(^{\circ }\) as the initial estimate, c logarithmic-scale representation of relative permeability

Full size image

3.3 Correlation of Contact Angle

In previous work, with a rock having more mixed-wet characteristics, including contact angles both above and below \(90^{\circ }\), it was shown that the spatial correlation of contact angle had a significant impact on the prediction of relative permeability (Foroughi et al. 2020). If the same distribution of contact angles was randomly assigned, the predicted properties were qualitatively different. We test this with these more oil-wet datasets. We retain the optimised contact angle distributions, Figs. 6 and 11, but assign them spatially randomly to pores and throats. Ten random realisations are compared to the results shown previously in Figs. 7 and 12.

We find that in our datasets the impact on predicted capillary pressure and relative permeability is modest, with only a slight downwards shift in the oil relative permeability at low water saturation for Ketton and a reduction in water relative permeability for the reservoir sample. This indicates that simply using a spatially randomly distributed contact angle around the thermodynamic values enables good predictions of capillary pressure and possibly relative permeability to be made (Figs. 13, 14).

Fig. 13
figure 13

a Capillary pressure and b relative permeability for the Ketton sample obtained using the optimised distribution of contact angle, see Fig. 6, compared to ten randomly redistributed realisations using the same distribution

Full size image
Fig. 14
figure 14

a Capillary pressure and b relative permeability for the carbonate reservoir sample obtained using the optimised distribution of contact angle, see Fig. 11, compared to ten randomly redistributed realisations using the same distribution

Full size image

In our previous study (Foroughi et al. 2020), we examined two types of correlation in the optimised distribution for wettability. Firstly, we determined the Pearson correlation coefficient between the contact angle and the inscribed radius of the elements. Secondly, the spatial correlation of the optimal contact angle distribution was evaluated. The Pearson correlation coefficient between two arbitrary variables x and y denoted by \(\rho _{xy}\) is defined as follows:

$$\begin{aligned} {\rho _{xy}=\frac{\sum \limits _{i}{(x_i-\bar{x})(y_i-\bar{y})}}{\sqrt{\sum \limits _{i}{(x_i-\bar{x})^2}}\sqrt{\sum \limits _{i}{(y_i-\bar{y})^2}}}}. \end{aligned}$$
(10)

Here, x is the contact angle for element i while y is the inscribed radius of the element (pore or throat). Note that the sum is over all pores or all throats, so we define correlation coefficient for pores and throats separately. \(\bar{x}\) and \(\bar{y}\) are average contact angle and average radius, respectively.

For Ketton, the Pearson correlation coefficients for pores are 0.09 and 0.10 for throats. For the reservoir sample, the corresponding values for pores and throats are 0.22 and 0.46, respectively. These results show that there exists a weak correlation between wettability and pore and throats size, especially for reservoir sample. For the Ketton sample, this geometric correlation is almost negligible. The weak correlation in Ketton explains why, during the displacement some smaller and possibly less strongly oil-wet pores are filled, in contrast to the strict filling in order of throat size observed in the reservoir sample (Scanziani et al. 2020; Alhosani et al. 2020a). In contrast, for the mixed-wet sandstone sample investigated previously the correlations were 0.70 and 0.68 for pores and throats respectively, demonstrating that there was a strong tendency for larger elements to be more oil-wet (Foroughi et al. 2020).

Secondly, we will test whether the wettability is spatially correlated with clumps of more water-wet and more oil-wet pores. The spatial correlation \(\xi (d)\) of contact angle from the centre of two elements i and j in the network model where d is the Euclidean distance between these two locations

$$\begin{aligned} {\xi (d)=\frac{\sum \limits _{j}\sum \limits _{i}{l_{ij}(\theta _i-\theta _j)^2}}{2\sigma ^2{\sum \limits _{j}\sum \limits _{i}{l_{ij}}}}}, \end{aligned}$$
(11)

where the sums over i and j are over all pores and all throats. \(l_{ij}(d)\) is defined as follows:

$$\begin{aligned} {l_{ij}(d)= \left\{ \begin{array}{ll} 1 &{} d-{\epsilon }< d_{i \rightarrow j} < d+{\epsilon } \\ 0 &{} \mathrm{otherwise} \end{array} \right. } \end{aligned}$$
(12)

\(\sigma ^2\) is the variance of the contact angle, \(d_{i \rightarrow j}\) is the Euclidean distance between centres of the element i and the element j. \(\epsilon\) is the statistical length, in this study a length of 2 \({\upmu } m\) was used for \(\epsilon\). Note that the value of \(\xi (d)\) for infinite range or high values of r converges to 1. \(\xi (d)\) = 0 indicates a perfect correlation between two spatially separated values, whereas \(\xi (d)\) = 1 indicates no correlation, and \(\xi (d)>\)1 represents an anticorrelation.

The results shown in Fig. 15 show no clear evidence of a correlation length beyond 100 \(\mu m\); again this contrasts with results on a mixed-wet sandstone where the correlation length was at least 100 \(\mu m\), spanning, on average, four pore lengths (Foroughi et al. 2020).

The overall conclusion of this analysis is that for the carbonate systems studied here, which have undergone a significant wettability alteration to render most surfaces contacted by oil oil-wet, there is no substantial correlation of contact angle in space, or a relationship between contact angle and pore or throat radius. This means that in modelling studies it is possible to simply assign contact at random, using a distribution centred on the thermodynamically measured average. Displacement is principally an invasion percolation process, with the water advancing through the largest available throats and then filling the neighbouring pores (Alhosani et al. 2020a; Wilkinson and Willemsen 1983). The displacement sequence, pattern of filling and the macroscopic properties, such as relative permeability, are controlled by the pore space geometry, rather than the exact arrangement of wettability.

In contrast, in a mixed-wet sandstone, where the contact angles are both smaller and larger than \(90^{\circ }\), and with a relatively homogeneous pore size distribution, the contact angles have a much more critical role on the filling sequence. Pores and throats of all size are filled at the same time controlled by local contact angle, with a preference to fill more water-wet regions (Foroughi et al. 2020; Lin et al. 2019). Capturing, correctly, exactly the spatial and size-dependent nature of the wettability was crucial to match the displacement sequence and make accurate predictions of relative permeability and capillary pressure (Foroughi et al. 2020; Lin et al. 2019).

In terms of macroscopic properties, the predicted and experimentally measured relative permeabilities for the previously studied sandstone displayed mixed-wet characteristics (Foroughi et al. 2020), with a low water relative permeability until a high water saturation was reached, and a crossover saturation, when the oil and water relative permeabilities were equal, of greater than 0.5, indicating favourable waterflood recovery. In contrast, the relative permeabilities for the oil-wet carbonate samples, Figs. 7 and 12, are classically oil-wet with crossover saturations of less than 0.5 with an implication of poor waterflood recovery (Blunt 2017).

Fig. 15
figure 15

Spatial correlation, \(\xi (d)\) of Ketton and reservoir carbonate samples. \(\xi (d)\) = 0 indicates a perfect correlation between two spatially separated values, whereas \(\xi (d)\) = 1 indicates no correlation, and \(\xi (d)>\)1 represents an anticorrelation

Full size image

4 Conclusions and Future Work

We used pore-network modelling to predict the displacement sequence of two synchrotron imaging datasets in the literature. In both cases, a sequence of images was acquired for waterflooding an oil-wet carbonate.

We applied an energy balance, assuming no viscous dissipation to estimate an average, thermodynamic, contact angle. This angle was initially applied to all the network elements, the pores and throats. We then applied a simple optimisation algorithm to adjust the angles on a pore-by-pore basis to minimise the mismatch between the observed and predicted filling sequence.

We found that for both datasets, the displacement was invasion percolation like with a preference to fill the largest available throats in the reservoir sample. For Ketton again the larger elements were filled preferentially but with some filling of smaller elements as well. We could match the filling accurately using a distribution of contact angle approximately centred on the initial average value with a standard deviation of approximately \(8^{\circ }\) or less.

With the optimised distribution of contact angle we predicted capillary pressure to within the uncertainty of the experimental measurements, based on estimates of interfacial curvature. We were able to predict both capillary pressure and relative permeability across the full saturation range. The relative permeabilities showed oil-wet characteristics with a low residual oil saturation, but a rapidly rising water relative permeability and a saturation crossover where the two relative permeabilities were equal of less than 0.5.

For the two oil-wet samples considered, we found that there was only a weak or negligible correlation between contact angle and pore size, and no evidence of a significant spatial correlation in contact angle. Using the matched distribution of contact angle distributed randomly in space had little impact on the predicted capillary pressures and relative permeabilities. This result contrasted with previous work on a mixed-wet sandstone with a relatively homogeneous pore size distribution (Foroughi et al. 2020). Here, the larger pores tended to be more oil-wet with a spatial correlation in contact angle of approximately four pore lengths. The distribution of contact angle had a critical impact on displacement sequence and the averaged flow properties—using randomly assigned angles—led to very poor predictions of relative permeability and capillary pressure. The significance of spatial correlation and the correlation between pore size and wettability therefore depends on the sample studied. Our results suggest that for the oil-wet carbonates studied these effects are less important than for a mixed-wet sandstone.

In further work, we could apply this methodology to wider range of flow conditions and rock samples. The work could provide guidelines for the necessary and sufficient conditions to make accurate predictions of capillary pressure and relative permeability. Furthermore, a more accurate assessment of wettability, allowing contact angles to be estimated on a pore-by-pore basis may remove the need to use the optimisation procedure to match the displacement sequence.