Local Capillary Pressure Estimation Based on Curvature of the Fluid Interface – Validation with Two-Phase Direct Numerical Simulations

With the advancement of high-resolution three-dimensional X-ray imaging, it is now possible to directly calculate the curvature of the interface of two phases extracted from segmented CT images during two-phase flow experiments to derive capillary pressure. However, there is an inherent difficulty of this image-based curvature measurement: the use of voxelized image data for the calculation of curvature can cause significant errors. To address this, we first perform two-phase direct numerical simulations to obtain the oil and water phase distribution, the exact location of the interface, and local fluid pressure. We then investigate a method to compute curvature on the oil/water interface. The interface is defined in two ways. In one case the simulated interface which has a sub-resolution smoothness is used, while the other is a smoothed interface which is extracted from synthetic segmented data based on the simulated phase distribution. Computed mean curvature on these surfaces are compared with that obtained from the fluid pressure computed directly in the simulation. We discuss the accuracy of image-based curvature measurements for the calculation of capillary pressure and propose the best way to extract an accurate curvature measurement, quantifying the likely uncertainties.


Introduction
Capillary pressure, Pc, is a pressure discontinuity across the interface between oil and water, defined as Pc = Po−Pw, where Po and Pw are the pressures of oil and water phase, respectively. Traditionally, capillary pressure for oil/water systems has been measured in a laboratory using the porous plate method in which the pressure of each phase is measured using two external pressure transducers. Based on the Young-Laplace equation, capillary pressure locally is defined as: where σ is the interfacial tension between two phases and κm is the mean curvature of the interface.
With the advancement of high-resolution threedimensional X-ray imaging, it is now possible to directly measure the curvature of the interface extracted from segmented CT images during two-phase flow experiments to derive capillary pressure. Armstrong et al. [1] demonstrated this approach using synchrotron-based tomographic datasets of oil/water drainage and imbibition cycles on a bead pack structure [2]. They compared the capillary pressure obtained from curvature measurements with that obtained from pressure transducers. Fairly good agreement was obtained for imbibition, whereas the curvature measurement showed a systematically lower value than that obtained from the transducers for drainage cycles. Later, using the same dataset, Li et al. [3] presented that their proposed curvature measurement method improved the agreement with the transducer based capillary pressure. Using a similar curvature measurement method, Herring et al. [4] estimated the capillary pressure for a range of curvature between 0 and 0.225 voxel -1 based on their air/water drainage and imbibition experiments on a Bentheimer sandstone. However, there is an inherent difficulty of this imagebased curvature measurement: the use of voxelized image data for the calculation of curvature can cause significant errors, resulting in a wide range of measured values, with some negative curvature values, which are not expected in a water-wet system. Hence, it is not clear how the distribution of measured curvature values represents the true range of local capillary pressure.
We investigate the accuracy of curvature measurement on the basis of pore-by-pore comparison using direct numerical simulations of two-phase flow. The color-gradient lattice Boltzmann method is employed to generate an oil and water phase distribution in pore structures of a bead pack and Bentheimer sandstone. From the simulated phase distribution, synthetic segmented data is generated, then curvature computation on the interface extracted from this segmented data is performed by employing several smoothing methods. These curvature values are then compared with those obtained from the simulated local fluid pressure. We discuss the accuracy of image-based curvature measurements for the calculation of capillary pressure and suggest the best method with associated errors for capillary pressure estimation.

The color-gradient lattice Boltzmann method
The color-gradient LB model proposed by Halliday et al. [5] was used. Our LB model was constructed with a 3D19Q lattice model which consists of a set of 19 discrete lattice velocity vectors, ei, in three-dimensional space. We defined particle distributions of two immiscible fluids, labeled red and blue, as fi r and fi b , respectively. The fluid density, ρr and ρb, and velocity, u, at position x and time t are obtained by: Using the color function, the interfacial tension between two fluids was computed as a spatially varying body force, F, based on the continuum surface force (CSF) model [7] given by: where σ is the interfacial tension and κ is the curvature of the interface. This spatially varying body force F was incorporated into the lattice Boltzmann equation through the body force term ϕi. For the MRT model, this was performed by transforming the forcing term proposed by Guo et al. [8] using the scheme presented in Yu and Fan [9]. After the application of the interfacial tension (F) to the particle distributions, the recoloring algorithm proposed by Latva-Kokko and Rothman [10] applied to these distributions to ensure phase segregation and maintain the interface. This results in a slightly diffusive interface whose thickness is about 2 to 3 lattice units. Further details of our LB model can be found in Akai et al. [11,12]. The only difference is that we used the MRT collision operator, while the single-relaxation-time (SRT) collision operator [13] was used in Akai et al. [11]. At solid-fluid boundary lattice nodes, a full-way bounce back boundary condition was implemented to achieve a non-slip boundary condition. In addition, the wettability of solid surface was modeled by specifying contact angles using the wetting boundary condition presented in Akai et al. [11]. This wetting boundary condition accurately assigns contact angles for 3D arbitrary geometries with smaller spurious currents compared to the widely-used fictitious density boundary condition [11].

Curvature computation on voxelized images
There are mainly two approaches to compute curvature from voxelized data such as micro-CT images. One approach calculates curvatures from the gradient of 3D float data (e.g., raw or processed grayscale data), while the other approach estimates curvature through the fitting of a quadratic form locally to a surface extracted from voxelized, segmented data [14]. We used the latter approach as this has been presented in previous studies [14,1,15,4,3,16].
We started with three-phase segmented label data (oil, water and solid) obtained from raw grayscale CT images. Using the marching cubes algorithm [17], the oil/water interface was extracted from the label data. Since this surface had a staircase shape, it had to be smoothed before computing curvature. In this study, we compare three smoothing methods: Constrained Gaussian smoothing (CGS), Laplacian smoothing (LPS) and boundary preserving Gaussian smoothing (BPGS). CGS is applied when the surface is extracted using the marching cubes algorithm. A Gasussian kernel filter with different kernel sizes is applied to label data, then isosurface is extracted. In this process, the constraint to preserve the location of an original label is imposed.
LPS [18] and BPGS [19] are applied after the extraction of the surface with the marching cubes algorithm. The extracted surface with a staircase shape is modeled as a triangulated surface. Then, the vertices of the triangle elements are moved with these smoothing methods. In LPS, the position of a vertex is moved to the average position of its neighboring vertices. BPGS also moves the position of a vertex based on the position of its neighboring vertices. A scale factor which determines the degree of the movement in one iteration is defined. Two consecutive smoothing steps with a positive and negative scaling factor are performed in one iteration. This smoothing produces a surface which preserves the original boundary without shrinkage [19]. The degree of the smoothing in LPS and BPGS is controlled by the number of iterations. In this study, CGS and BPGS were performed with commercial image processing software, Avizo, while LPS was performed with Paraview.
After the generation of a smoothed triangulated surface, the elemental triangles were fitted by a quadratic form: .  (8) Then, the principal curvature values and directions of principal curvature were obtained from the eigenvalues and eigenvectors of the fitted quadratic form in Eq. 8 [1]. Avizo was used to perform this computation. Users can choose the number of neighboring triangles to be used for the fitting at the center of a triangle. We used a fixed value of 4 neighbors in the following analyses.

Results
In section 3.1, we provide validation of our direct numerical simulations for different grid resolutions. Next, in section 3.2, using a single oil droplet test case, we investigate the accuracy of the curvature computation method described in section 2.2. Here, a smoothed interface which was extracted from synthetic segmented data based on the simulated phase distribution is used to compute curvature, then the resultant values are compared with those obtained from the simulated fluid pressure. Finally, in section 3.3, the same approach is applied to the simulated phase distributions in two realistically complex pore structures of a bead pack and Bentheimer sandstone.

Validation of the two-phase lattice Boltzmann model
To validate our two-phase lattice Boltzmann model and investigate its accuracy as a function of grid resolutions, the oil/water interface in a corner of triangular pore space was simulated. A 2D pore structure with an isosceles triangle shape as shown in Fig. 1 was used. Here, the length of Lx and Ly were set to 70 μm. R is the radius of curvature of the interface; θ is the contact angle; β is the half angle of the corner of the triangle, which is given by: tanβ=Lx/2Ly. This pore structure was modeled with 4 grid sizes of Δ=1.0, 2.0, 3.5 and 5.0 μm.
The identical density and viscosity of the water and oil phases were set to 1,000 kg/m 3 and 1 mPa, respectively. The interfacial tension and contact angle were set to 18 mN/m and 45°, respectively. Initially, the lower part of the pore structure was filled with oil to a specified oil saturation, while the other part at the top corner was filled with water. Then, simulations were performed for 50,000 time steps until they reached equilibrium.
In this pore geometry, the radius of curvature, R, can be analytically derived based on the geometrical relationship, which is given by: where Sw is the water saturation and the other parameters are shown in Fig 1. This analytically derived radius of curvature was used to compare with the simulated radius of curvature obtained using the input contact angle of θ=45°. Based on the simulated fluid configurations, the exact location of the interface corresponding to ρ N =0 (see Eq. 6) was extracted as shown in Fig. 2. Then, the radius of the curvature of this interface was obtained by fitting a circle to the interface, since in 2D in equilibrium the analytical shape of the interface is a part of a circle. Table  1 shows the comparison between the radius of curvature obtained from the analytical solution and that obtained with a circle fit to the simulated interface for different grid sizes. The relative error in the radius of curvature is less than 3% for these grid resolutions. This suggests that good agreement in capillary pressure is also obtained since the capillary pressure is directly linked to the radius of curvature through Eq. 1.

Curvature of a single oil droplet in triangular pore space
To investigate the curvature computation algorithm described in section 2.2, a simple 3D test case was performed. A cylindrical pore structure with the length of 171.5 μm with an isosceles triangular cross-section as in the previous test case was used. This pore structure was modeled with a grid size of 3.5 μm. The same fluid properties and contact angle as in the previous section were used. Initially, oil was placed in the central region of the pore space as shown in Fig. 3a. Then, simulations were performed for 50,000 time steps until they formed a single oil droplet after equilibrium (Fig. 3b).
Two types of the oil/water interface were prepared based on the simulation results: a "simulated interface" and a "smoothed interface". The simulated interface was obtained by extracting the contour line of the color function ρ N =0. This surface originally had a subresolution smoothness. The other surface, the smoothed interface, was obtained from synthetic voxelized label data. The simulated distribution of the color function was segmented into label data with the threshold of ρ N =0, then the oil/water interface was extracted using the marching cubes algorithm. This surface has a staircase shape due to the shape of a voxelized grid system. Therefore, it needs to be smoothed before computing curvature. As discussed in section 2.2, three smoothing methods were applied: Constrained Gaussian smoothing (CGS), Laplacian smoothing (LPS) and boundary preserving Gaussian smoothing (BPGS). The parameters used for the smoothing are summarized in Table 2.
Based on the simulated fluid pressure, the capillary pressure was obtained using: where Po avg and Pw avg are the average fluid pressure in the oil and water phases, respectively; No and Nw are the number of grid blocks having ρ N >0.9 and ρ N <−0.9, respectively; P is the simulated fluid pressure in each grid block. Here, to exclude the interface region, we computed an average pressure for each phase for |ρ N |>0.9. Then, Pc sim was converted to the simulated mean curvature using Eq. 1. We refer to this mean curvature as the fluid pressure derived mean curvature, κm P . Fig. 4 shows the oil/water interface after smoothing, colored by computed mean curvature. For all of the smoothed surfaces, we see a variation in the computed mean curvature, although the correct mean curvature should give a uniform value as shown in the curvature computed on the simulated interface, Fig. 4(j), because the oil droplet is in capillary equilibrium. Furthermore, we see errors around the edge of the interface (close to the three-phase contact lines) for all of the smoothed interfaces, Fig. 4(a)~(i). Therefore, we decided to discard the data points whose distance from solid surface is fewer than 3 voxels. This 3 voxels distance cutoff was used for all of the following analyses. Table 3 summarizes the average and standard deviation of the computed mean curvature and the relative difference to the fluid pressure derived mean curvature, κm P =0.133 voxel -1 for the smoothing level 2 as shown in Fig. 4(b), (e), (h). We note that Li et al. [3] proposed that in addition to the distance cutoff, taking an average weighted by the distance from the solid surface further improves the accuracy of curvature measurement using synthetic test cases of the interface in a cylindrical structure with a circular-cross section; however, taking a distance-weighted average did not improve the estimation of curvature in our cases as we can conclude from the overestimated curvature values observed in the central region of the oil droplet in Fig. 4(a) ~ (i).
In fact, the optimum smoothing method and its level of smoothing is dependent on the shape and the size of an object. Therefore, here, we only provide the qualitative features of the three smoothing methods. CGS preserves the shape of an object, but significant smoothing cannot be applied even with increasing the kernel size, because the resultant surface is constrained by the original voxel data. As a result, CGS tends to give a wider range of variation. LPS can apply significant smoothing by   increasing the number of iterations; however, this could change the shape of the interface as shown in a top part of the interface in Fig. 4(f) where we see the bending of the interface towards the opposite direction. This is caused by the propagation of errors in the three-phase contact line through many iterations. However, combining with the distance cutoff, this erroneous part can be effectively removed while preserving a well smoothed surface in the middle of the interface. BPGS gives the similar results to LPS, while better keeping the shape of an object compared to LPS. For the simulated interface which was directly obtained from the simulated color function without any smoothing, the computed curvature showed good agreement with the fluid pressure derived curvature with a smaller standard deviation (Table 3). Moreover, since this surface did not show the errors close to the threephase contact line, the average and standard deviation without the distance cutoff also gave similar values, i.e., the average of 0.132 voxel -1 and standard deviation of 3.99×10 -3 voxel -1 .

Local capillary pressure estimation for complex porous media
Two realistically complex pore structures were used: a synthetic bead pack pore structure and micro-CT images of Bentheimer sandstone. For these structures, simulation domains were composed of 256×256×256 voxels with a grid size of 3.56 µm as shown in Fig 5. For the analysis of the simulation results, the pore structures were divided into pore regions using the separate pore algorithm in Avizo, resulting in 402 and 272 pore regions for the bead pack and Bentheimer sandstone, respectively. The mean pore radius which accounts for 50% of the pore volume was 37 μm and 30 μm for the beadpack and Bentheimer sandstone, respectively.
The identical density and viscosity of water and oil phase were set to 1,000 kg/m 3 and 1 mPa, respectively. The interfacial tension and contact angle were set to 25 mN/m and 45°, respectively. All 6 faces of the cubic simulation domain were covered with solid voxels, i.e., a no-slip boundary condition was applied to all 6 faces. Initially, 50% of oil and 50% of water were randomly placed in pore voxels, then the simulations were performed with no external force for 250,000 time steps (corresponding to 0.1 seconds) until they reached equilibrium conditions. Fig. 6 shows the simulated phase distribution at the equilibrium condition. Similar to the analysis presented in the previous section, the simulated interface and the smoothed interface were prepared from the simulation results. Curvature computation was performed on these surfaces. Figs. 7 and 8 show the histogram of the computed curvature for the bead pack and Bentheimer sandstone, respectively. The fluid pressure derived mean Here, the surface is colored by the value of mean curvature. The correct surface should present uniform mean curvature as seen in the simulated interface since the droplet is in capillary equilibrium for which the capillary pressure is uniform. Table 3. Voxel based interfaces with the smoothing level 2 (see Table 2) and a distance cutoff of 3 voxels. Here, the relative difference was computed to the fluid pressure derived mean curvature of κm P =0.133 voxel -1 .   Table  4 for the bead pack and Table 5 for the Bentheimer sandstone.
Here, we discuss the following two key observations. First, the computed curvatures obtained with the smoothed interface show the wide range of the distribution when the cutoff is not applied, while the curvatures obtained with the simulated interface show a much narrower range. For the histograms without the cutoff, the difference between the distributions obtained with the smoothed interfaces and the simulated interface suggests that there are many erroneous values in curvature computation on the smoothed surfaces. The histograms obtained with CGS appeared to have the most similar distribution to that obtained with the simulated interface. However, as shown in Table 4 and 5, their average values are much higher than that obtained from the fluid pressure and simulated interface because of the long tails of their distribution toward values higher than 0.3, which are not shown in Figs. 7 and 8. Second, after the application of the distance cutoff, all the histograms obtained with the smoothed interface become similar to that obtained with the simulated interface. However, in both the smoothed and the simulated interface, the data points of high curvature values have been lost. This is because the distance cutoff removes the data points of not only the edges of the interface but also the entire parts of the interface in small pores, which tend to give a high curvature value. Although a high curvature value is not captured when the distance cutoff is applied, the smoothed interface seems to provide a good estimate of local mean curvature values of the interface. This will be   (a and b). The histograms computed on the simulated interface without and with the distance cutoff (c and d). The fluid pressure derived mean curvature obtained with Eq. 10 is also shown by the vertical line. further discussed in the following analysis. Among the three smoothing methods, when the cutoff is applied, LPS gave the closest average value of the mean curvature to that obtained from the simulated interface with the smallest standard deviation.
To further investigate the accuracy of curvature computation on the interface, the comparison was made on a pore-by-pore basis. First, we computed local capillary pressure applying Eq. 10 for each pore region. This local capillary pressure was converted to the mean curvature for each pore region using Eq. 1. Next, we also obtained the mean curvature for each pore region by taking the average of the computed curvature values on the interfaces in that pore region. Figs. 9 and 10 show the comparison between the dimensionless mean curvature obtained from the simulated fluid pressure and that obtained from the computed curvature on the interface for the bead pack and Bentheimer sandstone, respectively. In both, the curvature computed on the simulated interfaces gave consistent values with those obtained from the fluid pressure. This means that when the interface is reasonably smooth, the curvature computation gives an accurate estimation of local capillary pressure. For the curvature computed on the smoothed interface, CGS shows highly overestimated curvature values for some pore regions, while LPS and BPGS estimate curvature within a range of about ±20% difference from the curvatures obtained with the fluid pressure for the mean curvature smaller than

Conclusions
Direct numerical simulations were conducted to obtain oil and water phase distribution, the exact location of the oil/water interface and local fluid pressure. From the simulation results, two types of the oil/water interface were generated. One is a simulated interface which has a sub-resolution smoothness, while the other is a smoothed interface. The smoothed interfaces were generated by applying three types of smoothing methods to a staircase shape surface extracted from synthetic segmented data based on the simulated phase distribution. Then, curvature computation was performed on these surfaces. The resultant mean curvature value from these surfaces was compared with that obtained from the simulated fluid pressure. This was performed to find the best way to extract accurate curvature measurement, and to quantify the likely uncertainties.
First, we tested our approach using a single oil droplet in triangular pore space. Computed curvature values on the simulated interface showed a narrow distribution whose average value was consistent with that obtained from the fluid pressure. This means that when the surface is sufficiently smooth, the fitting of a quadratic form equation to the surface properly computes local curvatures. On the other hand, computed curvature values on the smoothed interface showed a wide distribution with erroneous values close to the three-phase contact lines. Therefore, we discarded values whose distance from solid surface was fewer than 3 voxels. After the application of the distance cutoff, the average of computed mean curvature values became closer to that obtained from the fluid pressure. However, their standard deviation was on the order of 10 -2 voxel -1 against the mean curvature value of 0.133 voxel -1 , which was 10 times higher than that obtained from the simulated interface. These deviated values were caused by error, since an oil droplet in capillary equilibrium should give a uniform value of mean curvature.
Next, realistically complex pore structures of a bead pack and Bentheimer sandstone were used to simulate oil droplets in capillary equilibrium. When the distance cutoff was not applied, curvature computed on the simulated interface gave an average value of mean curvature consistent with that obtained from the fluid pressure for both the bead pack and Bentheimer sandstone. However, when the distance cutoff was applied, the average value became lower than that from the fluid pressure. This is because the distance cutoff removes the data points of not only the edges of the interface but also the entire parts of the interface in small pores, which tend to give a high curvature value. For the smoothed interface, the distributions obtained without the cutoff were quite different from that obtained from the simulated interface. Many negative values of computed mean curvatures suggest significant errors in these cases.
After the application of the distance cutoff, the distribution became closer to that obtained with the simulated interface.
Finally, we compared computed mean curvature on a pore-by-pore basis. As opposed to the test case performed on the single oil droplet, oil droplets with different mean curvature values corresponding to different pore sizes can be obtained in these simulations. Therefore, we need to investigate whether the range of computed mean curvature is caused by error or it reflects a variation in local capillary pressure. Good agreement between mean curvature obtained from the fluid pressure and that computed from the simulated surface for each pore region suggested that the range of distribution observed in the simulated surface properly captures variation in local capillary pressure. For the smoothed interface with the cutoff, LPS and BPGS estimated curvature within a range of about ±20% difference from the curvatures obtained with the fluid pressure for mean curvatures smaller than about 0.15 voxel -1 , while CGS showed highly overestimated curvature values for some pore regions.
In conclusion, the application of the distance cutoff is necessary to remove the errors close to the three-phase contact line. Among the three tested smoothing methods, LPS appeared to be the best method since it gave the closest average of mean curvature values to that obtained from the simulated interface with the smallest standard deviation. In addition, it gave good estimates of the local mean curvature for each pore with a ±20% error up to 0.15 voxel -1 at this resolution.
Currently, we are using the similar approach to quantify the accuracy of curvature computation for drainage and imbibition events on the bead pack and Bentheimer sandstone. This will be discussed in a future publication.