Pore-Morphology-Based Simulation of Drainage in Porous Media Featuring a Locally Variable Contact Angle

Since the first publications by Hazlett (Transp Porous Med, 20:21–35, 1995) and Hilpert and Miller (Adv Water Res, 24:243–255, 2001), the pore-morphology-based method has been widely applied to determine the capillary pressure–saturation curves of porous media. The main advantage of the method is the simulation of a primary drainage process for large binary images using moderate computational time and memory compared to other two-phase flow simulations. Until now, the pore morphology model was restricted to totally wetting materials or those with a constant contact angle. Here, we introduce a similarly computationally efficient extension of the model that now enables the calculation of capillary pressure–saturation curves for porous media, where the contact angle varies locally within, due to a composite structure.


Introduction
Following the first work based on pore morphology or full morphology (FM) models, the FM approach has been used in several studies to determine the capillary pressure-saturation curves for different porous media. This includes soil samples (e.g., Silin et al. 2011) as well as materials for technical applications like the gas diffusion layer (GDL) and micro-porous layer (MPL) of the polymer electrolyte fuel cell (PEFC), as shown in Schulz et al. (2007) and Wargo et al. (2013). The key advantage of the model consists of the low memory requirements and the extremely low calculation times, especially when compared to other multiphase flow simulations such as lattice-Boltzmann methods and the level set approach recently presented by Jettestuen et al. (2013).
The FM model delivers results comparable to pore network models and full multiphase flow simulations, as numerous publications have recently shown (Vogel et al. 2005;Mukherjee et al. 2011). While the computation time and memory requirements are comparable between pore network and FM models, pore network models require a transformation of the 3-D microstructure into an idealized network of pores. In contrast, the FM approach can be applied directly to any given binary microstructure regardless of the porosity and shape of the pore space. However, it should be noted that the FM approach only models the quasi-static state of a drainage process.
Early FM simulations assumed a totally wetting fluid, which corresponds to a contact angle of θ = 0 • . However, in practice, media with non-vanishing contact angles are used frequently, particularly for technical applications. Additionally, porous materials often exist as composites in which each constituent solid has a unique wettability. A number of experiments show the influence of wetting conditions on the capillary pressure-saturation relation (Moseley and Dhir 1996;O'Carroll et al. 2005). Pore network models have also been developed to accommodate structures with mixed wettability (Valvatne and Blunt 2004;Ryazanov et al. 2014). Therefore, it is desirable to extend the FM model, which can analyze a digital structure directly on a pixel-by-pixel basis, to handle porous composites with locally variable surface wettability.
To date, partially wetting fluids in the FM model have only been taken into account through a rescaling of the capillary pressure with a factor of cos θ Schulz et al. (2007). This paper outlines how to incorporate locally variable contact angles in the FM model and is structured as follows. First, the background of the established FM model following Hilpert and Miller (2001) and Schulz et al. (2007) is described. Subsequently, an extension of the FM model for partially wetting fluids using a constant contact angle over the entire porous medium is introduced. The extension of the FM model is then demonstrated for the simulation of porous media, where the wetting conditions differ locally within the structure. The influence of locally variable contact angles on the capillary pressure-saturation relation is exemplified.

Approach
The FM model, which has been established by Hilpert and Miller (2001), simulates the quasistatic steps of a drainage process in totally wetting porous media. Initially, the porous medium is filled entirely with the wetting fluid. The wetting fluid is connected to a wetting phase (WP) reservoir on one side of the medium, while the opposite side faces a non-wetting phase (NWP) reservoir. Boundary conditions for the remaining sides can be chosen freely, e.g., consisting of closed walls or periodic boundaries. The choice of these boundary conditions impacts the determined capillary pressure-saturation curves only marginally given that the size of the sample is large enough (Schulz 2003).
The FM approach (in 2-D) uses a digital disk with a diameter D as a probe in order to determine the pore space that is accessible by the NWP. The diameter D of the disk corresponds to a capillary pressure of where γ is the interfacial tension. The disk moves from the NWP reservoir into the pore space without intersecting the solid phase. The pore space that is accessible to the NWP at the corresponding capillary pressure p c from Eq. 1 is thus determined. The FM model neglects some physical effects such as trapped or irreducible wetting phases. Nevertheless, it is very close to the experimental setup of mercury intrusion porosimetry (Becker et al. 2008;Silin et al. 2011) and has been shown to produce results comparable to pore network modeling and more computationally intense flow simulations (Vogel et al. 2005;Mukherjee et al. 2011). Figure 1 shows a snapshot of a typical FM simulation in 2-D, for which the porous medium consists of non-overlapping disks (gray) with the same diameter. From left to right, the radius of the digital probing disk is (a) D = 18 pixels, (b) D = 14 pixels, and (c) D = 10 pixels. Using Eq. 1, the corresponding capillary pressure for an air-water system with a surface tension of γ = 0.0723 N/m and resolution of = 0.5 mm/pixel is calculated. This yields capillary pressures of (a) p c = 16.07 Pa , (b) p c = 20.66 Pa and (c) p c = 28.92 Pa.

Implementation of the FM Model
For practical reasons, the implementation of the FM model is decomposed into several calculation steps. The approach is based on digital image processing methods and incorporates erosion and dilation procedures. The steps of the general FM model following Schulz et al. (2007) are listed below, and specific details on the morphology operations are explained in Sect. 2.2.
1. The pore space is eroded by digital disks with increasing radius starting with the smallest radius. 2. From the erosion result, all parts which are not connected to the NWP reservoir are removed. 3. The eroded set of remaining connected pores is considered accessible to the NWP and is dilated using the same digital disk used in step 1. The saturation is then recorded as the volume fraction of the pore space filled by NWP. 4. The entire process is repeated for the next pressure step using a larger disk.

Mathematical Expression
As shown by Hilpert and Miller (2001), the steps in the FM model can be expressed in a formal way by introducing the Minkowski addition (⊕) and Minkowski subtraction ( ). In general, the Minkowski addition of a binary image A is defined as where B is called a structuring element (Ohser and Mücklich 2000). In the FM model (in 2-D), B = S(D) is a digital disk with even diameter D. Since S(D) is symmetric, the Minkowski addition is equivalent to a morphological dilation of the image A (Ohser and Mücklich 2000). Additionally, the Minkowski subtraction of image A by structuring element B is given by Again, with the choice of B = S(D), i.e., a digital disk with diameter D, Eq. 3 is the morphological erosion of the image A. 1 Using the framework of morphological operations, the NWP saturation at capillary pressure p c (Eq. 1) is determined by In Eq. 4, P S(D) describes the erosion of the pore space P using a disk of diameter D.
Additionally, C [P S(D)] is defined as the part of P S(D) that is connected to the NWP reservoir. The result of this operation is dilated with the disk of diameter D, expressed as ⊕ S(D). Eq. 4 thus corresponds to the steps outlined in Sect. 2.1. The saturation of the NWP is eventually calculated as the volume fraction, where Vol[P] is the volume of the pore space. The morphological operations with disks are closely related to the concept of parallel bodies as shown by Mecke (2000). Therefore, an efficient implementation of the FM model relies on the Eucledian distance transform of a digital image. An efficient algorithm for the calculation of the Eucledian distance map has been presented by Saito and Toriwaki (1994).

Previous Work
In previous publications (Schulz et al. 2007;Becker et al. 2008;Wargo et al. 2013), we already proposed an extension of the FM model to porous media which are partially wetting with a constant contact angle. In these studies, the contact angle differs from 0 • , and the capillary pressure is given by the corresponding Young-Laplace equation where θ is the contact angle. Using Eq. 5, Schulz et al. (2007) and Becker et al. (2008) were able to align FM model results with measured data for fuel cell GDL samples by rescaling the capillary pressure with cos θ during post-processing of the results.

Further Model Extension
This paper proposes further modification of the FM model to treat partially wetting porous media directly in the simulation, by instead scaling the diameter of erosion by | cos θ |. As outlined below, this important update to the FM model allows for the direct analysis of a partially wetting structure, including cases where the contact angle is constant throughout the entire medium, as well as cases, where the contact angle varies locally within the pore structure. For demonstration, the modified model is first applied to cases, where the contact angle is larger than 0 • but constant. This model extension is based on the introduction of a second diameter D , which is defined by By employing D in the erosion procedure, the corresponding NWP saturation is now calculated by Note the slight difference between Eqs. 4 and 7, in which S(D) becomes S(D ). Accordingly, the modified FM model determines the pore space that is accessible to the NWP if the pore diameter is larger than D (instead of D). This is illustrated in Fig. 2 for a 2-D capillary, where the position and shape of the interface for totally and partially wetting conditions are compared. In a totally wetting medium, the shape of the interface in the pore neck is a semicircle (Fig. 2a). However, at larger contact angles the interface is only a segment of the semicircle, as illustrated in Fig. 2b for a value of θ = 20 • .
It should be noted that the model extension described in Eqs. 6 and 7 does not influence the computational efficiency of the FM approach. The memory consumption increases only to the extent necessary for the coding of locally different wetting conditions, as described in Sect. 3.4.

Drawbacks of the Extended FM Model
As can be seen in Fig. 2b, a straight forward dilation of the eroded pore space with the disk of diameter D would also fill some of the solid with non-wetting fluid. Therefore, the dilation is not applied to the solid phase. In Eq. 7, this is expressed by the modified Minkowski addition ⊕ P , indicating its restriction to the pore space P. In practice, this can lead to artifacts if the 'thickness' of the solid phase is smaller than (D − D )/2, in which NWP could be erroneously dilated into pore space beyond the solid boundary. These artifacts could be mitigated through an additional connectivity check after the dilation step, but represent a drawback which is getting more severe if the contact angle tends toward 90 • .
Additionally, a contact angle close to 90 • corresponds to a large value of the fraction D/D as given in Eq. 6. As a consequence, there will be a large gap between the position of the pore throat and the resulting interface between wetting and non-wetting phases (see Fig. 2b). This can result in an over-prediction of the NWP saturation and shifting of the WP|NWP interface. Thus, in the following examples, the contact angle is varied between 0 • and 60 • to minimize these effects.
Furthermore, it should be noted that although Fig. 2b shows that the extended FM model can correctly predict a partially wetting meniscus, the straight capillary walls of the test case in Fig. 2 will not exist in most practical applications. For more realistic cases (e.g., soil, fuel cell GDL), the pore walls will be jagged/curved, and the shape of menisci will inevitably fail to match the given contact angle visually. Therefore, it should be stressed that the primary strength of the extended FM approach is its ability to calculate capillary pressure-saturation curves for large datasets with locally varying contact angle in a rapid, computationally efficient manner.

Simulation of Locally Different Contact Angles
Perhaps the most important advantage of the presented extension is that the FM model can now handle locally different contact angles. In this case, instead of a binary image, the model input is a grayscale image, where each grayscale value corresponds to a distinct material phase i with contact angle θ i . The main difference in the calculation steps (outlined in Sect. 2.1) is the choice of the diameter used for the erosion of the pore space. Now that there are multiple contact angles within the structure, the diameter of the structuring element may change size when moving from pixel to pixel during the erosion procedure due to the locally different contact angles. Figure 3 compares the phase distribution results at the same capillary pressure for a medium of solid disks with constant contact angle (60 • ) versus disks in which the contact angle varies locally on the surface (0 • and 60 • ). For the latter case (Fig. 3, bottom), random pixels on the surface of the disks are assigned a contact angle of 0 • . All other pixels are assigned a contact angle of 60 • . As shown in Fig. 3c, the random pixels result in tiny surface protrusions of WP after the erosion procedure due to the increased wettability of the associated surface pixels. The influence of the random surface pixels is captured by the extended FM model and yields significantly different saturation results when comparing Fig. 3b versus Fig. 3d.
It should be noted that for input structures with multiple contact angles, it may be necessary to consider only the 'surface' pixels of the solid (i.e., pixels adjacent to the pore space) during the erosion procedure. This was not necessary in the case of Fig. 3c, since the erosion diameter associated with the randomly selected surface pixels (D = 28) is larger than the erosion diameter for the remaining solid (D = 14). However, if the associated contact angles were reversed, the influence of the random surface pixels would be eclipsed by the adjacent internal disk pixels. For such a case, only the surface pixels should be considered during the erosion. The surface pixels can be identified easily by eroding the solid features of the input image using a disk radius of 1 pixel and then subtracting the resulting image from the original.

Application to a Test Case with a Constant Contact Angle
As described above, the modified FM model can calculate the capillary-pressure saturation curves of a porous medium with different wetting conditions. While the approach can easily be applied in 3-D using spherical structuring elements, only 2-D test cases are presented herein for clearer demonstration and illustration purposes. Thus, we have chosen a digital porous medium with dimensions of 2,000 × 2,000 pixels. The structure consists of 1,235 nonoverlapping solid disks each with a diameter of 36 pixels, which corresponds to a porosity of 70.0 % as shown in Fig. 4.

Capillary Pressure-Saturation Curves of the Test Case
The extended FM model is applied to the test case shown in Fig. 4 to calculate the capillary pressure-saturation ( p c − s n ) curve of a primary drainage process, assuming a constant contact angle over the entire solid phase. Four simulations were performed using contact angles ranging from totally wetting (θ = 0 • ) to a maximum value of θ = 60 • . As in the example shown in Fig. 1, a resolution of = 0.5 mm/pixel is chosen and the capillary pressure is determined for an air-water system (surface tension γ = 0.0723 N/m). This yields capillary pressures in the range of 5-25 Pa as shown in Fig. 5. Interestingly, the shape of the p c − s n curves is almost constant and independent of the contact angle, and only the The points of the p c −s n results in Fig. 5 correspond to FM model simulations with digital disks of even diameter (D ∈ (2 · n), n = 2, 3, 4, . . .). This limitation is necessary to ensure that smaller disks are strict subsets of larger ones as discussed by Hilpert and Miller (2001). The limitation to even diameters leads to fewer steps in the drainage simulation but can be compensated by increasing the resolution of the image.

Interpretation for Partially Wetting Porous Media with a Constant Contact Angle
As shown for the test case, the results obtained by the extended FM model appear similar to a rescaling of the capillary pressure with the cosine of the contact angle, as described in Section 3.1. For the investigated contact angles of up to 60 • , only a slight difference is shown in the NWP saturation for the scaled capillary pressure results (Fig. 5, right). This result confirms that the rescaling of the capillary pressure with cos θ as proposed by Schulz et al. (2007) is sufficient in the case of a constant contact angle. The small increase in discrepancy as the NWP saturation increases from right to left occurs due to the limited resolution when working with structuring elements of increasingly smaller size, since the scaled diameter (D ) must Fig. 6 The porous medium from Fig. 4 with locally different wetting conditions. The gray disks are totally wetting (θ = 0 • ), and the black disks are partially wetting (θ = 60 • ) be rounded to a nearest integer prior to performing the erosion procedure on the pixelated structure. This rounding error is amplified as the size of the morphological operating diameter decreases.
Other publications demonstrate a similar influence of the contact angle on the primary drainage curves. O'Carroll et al. (2005) have investigated Ottawa sand size fractions with different wetting conditions, and the reported primary drainage curves show a similar trend towards smaller capillary pressures as contact angle is increased. Jettestuen et al. (2013) used a level set approach to simulate the drainage in regular sphere packages. In their work, the simulated domain is rather small with only four spheres resolved by 66 × 66 × 62 pixels. Nevertheless, Jettestuen et al. (2013) observed the same trend of decreased entry pressure as contact angle increases. They report that for varying contact angles, the main differences occur at very low saturation. At these low saturation values, the form of the interfaces between the non-wetting and wetting phase is dominated by pendular rings. As already mentioned by Hilpert and Miller (2001), FM results for very low saturations have to be carefully interpreted when pendular rings are dominant. In consequence, a qualitative comparison of the results in Fig. 5 to the data from Jettestuen et al. (2013) is not possible.

Application of the Extended FM Model to Porous Media with Locally Variable Wetting Conditions
The main advantage of the described model extension is its application to porous media where the contact angle is not constant. To help demonstrate the influence of a locally different contact angle, the 2-D porous medium from Fig. 4 has been modified using disks with two separate wetting conditions, shown in Fig. 6. In the first configuration (Fig. 6, left), all disks with an x-coordinate to the left of the vertical center (x <1,000) are assumed to be totally wetting (θ = 0 • ), and the remaining disks on the right side are assigned a contact angle of θ = 60 • . In the second configuration (Fig. 6, right), approximately 50 % of the disks are randomly assigned a contact angle of 60 • while the remainder are totally wetting. As described in Sect. 3.2, the extended FM model allows for a locally different calculation of the erosion. The erosion of the pore space is carried out via a disk of diameter D if the solid is totally wetting. If the solid is only partially wetting, a disk of scaled diameter D < D (Eq. 6) is used. A separate D is determined for each contact angle present in the dataset. Fig. 7 Snapshot of NWP distribution in a porous medium with locally different contact angles. The solid on the left side in light gray is totally wetting. The contact angle for the dark gray solid in the right half is 60 • . The NWP in black enters the pore space from the top. The WP is white After the connectivity check, all remaining NWP is dilated with the unscaled diameter D. Finally, the capillary pressure is calculated using Eq. 1, in this case with a resolution of = 0.5 mm/pixel and γ = 0.0723 N/m. An example of the phase distribution resulting from the extended FM model for the left configuration of Fig. 6 is shown in Fig. 7. At the same capillary pressure, notably different saturations are shown for the pore space of the totally wetting half (light gray disks on the left side) versus the partially wetting half where θ = 60 • (dark gray disks, right). The NWP saturation distribution in Fig. 7 was determined using digital disks of diameter D = 46 pixels and D = D · cos 60 • = 23 pixels in the morphological operations. The right half of Fig. 7 shows how the pore space around the partially wetting solid is more accessible to the NWP than the left side which is totally wetting.

Capillary Pressure-Saturation Curves for Porous Media with Locally Variable Wetting Conditions
The distribution of the disks with different wetting conditions has a strong impact on the capillary pressure-saturation relation as shown in Fig. 8. The configuration of the vertically separated totally wetting and partially wetting disks (Fig 8, 'vertical') yields an entry pressure of about 6 Pa. This is the same pressure found for the porous medium with a constant contact angle of 60 • (see Fig. 5, left). At a saturation of about 50 %, the partially wetting porous solid is entirely filled with the non-wetting fluid. The totally wetting solid fills with the NWP only for capillary pressures larger than approximately 10 Pa, causing a step to appear in the capillary pressure-saturation curve around a saturation of 0.5 for the 'vertical' configuration. The p c − s n curve is, therefore, similar to a parallel arrangement of two media with equivalent porosity but different pore size (i.e., Miller-similar medium). The significant differences between the 'vertical' and the 'random' configuration illustrate the influence of locally different contact angles on the wetting conditions. Interestingly, when the disks with different wetting conditions are placed randomly, the capillary pressure-saturation curve is similar to the one calculated on the basis of a constant contact angle of approximately θ = 40 • (see Fig. 8, 'random'). The value of 40 • is close to an average of the cosines of the contact angles θ = arccos cos 0 • + cos 60 • 2 = arccos(0.75) ≈ 41.41 • , which is believed to better capture the average effect of the random mixed case shown on the right in Fig. 6.  Fig. 8 Capillary pressure-saturation curves for the test media shown in Fig. 6. The contact angle differs locally between 0 • and 60 • . The 'vertical' configuration corresponds to the image on the left side in Fig. 6 and 'random' to the image on the right, respectively. For comparison, the curve for a constant contact angle of 40 • is also shown

Summary
The capillary pressure-saturation relation is an essential property for the description of multiphase flow in porous media. This property is strongly influenced by the wettability of the solid, i.e., the contact angle, which has been shown previously in several experiments and simulations. The full morphology (FM) model has been established as a computationally efficient tool to determine the p c − s n curves of 2-D and 3-D binary images of porous media. The FM approach simulates the primary drainage of a wetting/non-wetting system and until now, was limited to totally wetting porous media and those with a constant contact angle. This study presents an extension of the FM model which takes into account locally different contact angles θ i . The extension is implemented directly through the morphological erosion procedure, which is carried out in the extended FM model using a disk(s) of diameter D (scaled by | cos θ i | for each solid phase i) instead of D, where D < D. This modification has no major influence on the computational efficiency and memory consumption of the FM model.
For validation purposes, the extended FM model is first applied to partially wetting structures with different uniform contact angles, and the resulting p c − s n curves are compared. The results confirm that the different uniform wetting conditions correspond to a scaling of the capillary pressure by a factor of 1/ cos θ , which validates the pressure-scaling approach previously proposed by Schulz et al. (2007) for FM-based models. The extended FM approach was then demonstrated on structures in which the contact angle varies locally, which is the primary advantage of the model extension. For a structure containing two different contact angles (e.g., 0 • and 60 • ), the simulated capillary pressure-saturation curves are found to depend strongly on the specific distribution of wettability. The presented examples demonstrate that the extended FM model offers a rapid, computationally efficient approach for determining the capillary pressure-saturation relations of porous materials with locally variable surface wettability.
Although the application of the extended FM model to porous media with contiguously variable contact angles (i.e., in which contact angle varies along continuous pore walls) is not exemplified specifically by the structures tested in this paper, the findings of this study demonstrate the potential of this approach to handle virtually any number of contact angles existing in a single dataset. Such a spatial distribution of contact angles (i.e., constituent mate-rials) in a porous composite can be measured directly using recent advances in tomographic imaging techniques, generated based on knowledge of the structure, etc.
Additionally, it is worth noting that, while the examples presented in this work were limited to 2-D cases to allow for clearer demonstration and illustration of the approach, there is no significant difference in the application to 3-D images. In the case of 3-D, the structuring element that is used in the erosion and dilation of the image is a digital sphere of diameter D or D (Eqs. 6 and 7). Also, as in 2-D applications, the Euclidian distance map may be used to calculate the erosion and dilation of the image efficiently in 3-D. When based on the algorithm presented by Saito and Toriwaki (1994), this computation is O(n 4 ) for a 3-D binary image with n = n x = n y = n z . However, algorithms of higher efficiency in the order of O(n · log n) exist which are based on solving an eikonal equation (Surazhsky et al. 2005).