© L.M. Keller, published by IFP Energies nouvelles, 2021

1 Introduction

In the present work, model pore structures of Opalinus Clay were used to predict critical material properties related to gas transport. This topic is important, for example, because rocks such as Opalinus Clay are considered as potential host rocks for the disposal of radioactive waste (Andra, 2005; Nagra, 2002, 2004). Thereby, these transport properties are of interest in connection with the production of gas, for example through corrosion of the waste containers, and its subsequent transport through the rock material surrounding the radioactive waste. The produced gas can be transported from the place of formation through the rock by various transport processes (e.g. diffusion, advection, two-phase flow, etc.). Here, the focus is on the so-called two-phase transport where the produced gas displaces the water in the pores as a separate phase. This process is known as drainage. Two-phase flow is controlled by the microstructure of the porous medium, in particular by the characteristic pore throat radius and its distribution.

During drainage, the water in the pores is successively replaced by gas. This reduces the volume fraction of connected water bodies, which extend over the whole sample. As a consequence, the effective permeability of water is thereby continuously reduced and the effective gas permeability increases continuously as soon as a connected gas transport path is formed across the sample. Furthermore, there is a pressure difference between the phases, the capillary pressure that a gas must exceed in order to be transported into a water-saturated material. Transport of gas by two-phase flow depends on the following material properties of the porous medium: (i) the permeability, (ii) the porosity, (iii) the effective permeability and (iv) capillary pressure as function of saturation.

These properties are controlled on the pore scale but estimates of the suitability of a potential host rock require evaluations of gas transport on a larger length scale, where micropores cannot be considered as physical objects. Large-scale numerical simulations are often used to evaluate a potential host rock, which in turn require the aforementioned gas transport properties as input. These properties can be determined by laboratory experiments, but the available data related to Opalinus Clay are sparse. For example, water retention curves were determined to evaluate relationship between capillary pressure and saturation (Ferrari and Laloui, 2012; Muñoz et al., 2003; Nagra, 2013; Romero and Gómez, 2013; Zhang and Rothfuchs, 2007). Despite these laboratory data it remains unclear to what extend anisotropy affects the relationship between saturation and capillary pressure. In addition, data availability is poor in case of effective permeabilities and it is also unclear to what extent anisotropy affects effective permeability (Nagra, 2013). The aim of this study is the prediction of gas transport properties based on pore microstructures, which were reconstructed from 3D image data and in combination with numerical simulations. The predicted gas transport properties will be compared with experimental data. For those properties for which no experimental data are available (i.e. effective permeability during drainage, effect of anisotropy), the study makes predictions that hopefully will ignite further experimental work.

On the microscale, pores can be observed within the clay rich parts of Opalinus Clay. There, pores are formed due to geometric incompatibilities between the grain boundaries of small clay platelets. These pores are commonly referred to as intergranular pores and the radii associated with these pores range from about 2 nm to around 100 nm. In 3D these pores were imaged using Focused Ion Beam (FIB) nano tomography, which allowed reconstructing pore microstructures on the tens of micrometer scale (Keller et al., 2013). Unfortunately, FIB cannot resolve the whole pore space (Keller et al., 2011). With a typical voxel size of 5–10 nm about 30% of the pores can be imaged (Keller et al., 2011). In relation to the objectives of this study, this has the consequence that material properties can only be calculated incompletely. To overcome this problem we assume that the pore morphology in Opalinus Clay is approximately similar for the whole pore size spectrum, which spans over two orders of magnitude. This allowed the construction of a digital pore structure, which covers almost the complete pore size spectrum.

2 Derivation of two-phase flow properties from pore models

2.1 Image data

Only published image data were used in this study (Keller et al., 2013). Data from a sample of the shaley facies of Opalinus Clay (BDR sample) were used. The sample was taken from Opalinus Clay unit at the Mont Terri Rock Laboratory, Switzerland (Fig. 1, Tab. 1). For sample BDR also pore size data, which determined by using the nitrogen gas–adsorption technique in combination with the BET theory (BET), are available (Keller et al., 2013, Tab. 1).

Fig. 1 Reconstructed microstructures based on FIB. See Table 1 for voxel resolution and porosity. (a) The 3-D reconstruction of the analyzed volume of OPA from the shaley facies based on SE (secondary electron) images. (b) The 3-D reconstruction of the analyzed volume of OPA from the sandy facies based on SE (secondary electron) images.

All used methods (imaging methods, image processing and characterization, etc.) are described in detail in our previous publications. Here, only a short summary is given. For electron microscopy analysis clay samples need to be dried prior to analysis. Conventional air-drying causes artifacts such as the frequently observed desiccation cracks, which can lead to misleading conclusions in terms of transport paths in clay rocks. In addition, freeze-drying of moist materials may cause preparation artifacts such as ice formation or surface roughness during mechanical polishing. Special methods such as high-pressure freezing and subsequent freeze-drying were used in order to avoid these artifacts (Bachmann and Mayer, 1987; Keller et al., 2011). The porosity that can be resolved by FIB corresponds to about 20–30% of the total pore space (Keller et al., 2011, Tab. 1). With a typical voxel size of 5–10 nm only pores >5–10 nm can be imaged. These comparatively larger pores control the gas entry pressure and the capillary pressure at the beginning of the drainage of water from the pores. Hence, the capillary pressure curve determined on the base of FIB image data is largely incomplete because a large part of the present pore sizes are not considered.

2.2 Building digital pore structures

The calculation of the capillary pressure curve and related permeabilities requires as input a pore structure representing a large part of the pore size spectrum. In order to construct such a digital pore structure we used a partial volume of the 3D image data as a starting point. This volume had a size of 300 × 300 × 300 voxels at a voxel size of 10 nm and a porosity of 2.72 vol.%. BET pore size data show that the investigated sample from the shaley facies had a porosity of about 12.0 vol.%. Therefore, the porosity of the starting pore structure was increased in order to account for the resolution limitation of FIB (see above). This has been achieved by the following procedure:

• Argue that the pore structure is self-similar.

• If the condition that the pore structure is self-similar would be satisfied, this would justify the superposition of a finer pore structure with coarser pore structures, where the pore structures are topologically identical (see also Wang et al., 2014).

In order to provide evidence that the clay pore structure is self-similar we follow the approach of Katz and Thompson (1985). These workers correctly predicted porosity from fractal dimensions and used it as an argument that the pore structure is therefore a fractal. Yu and Li (2001) analytically derived an expression between fractal dimension and porosity as,

(1)where d is the Euclidian dimension (d = 3 in 3D), D _f is the fractal dimension and λ is the pore size scale (from smallest (λ _min) to the largest diameter (λ _max)) and the statistical self-similarity exists for a set of fractal pores in this range. In addition, Kou et al. (2009) point out that the relation λ _min ≪ λ _max must be satisfied for a media to be considered as fractal. Regarding the pore structure of sample BDR we have from BET measurements λ _min = 1.4 nm and λ _max = 60.8 nm, which yields λ _min/λ _max = 0.02, which according to (Yu and Cheng, 2002) satisfies the previous relation.

In this study the fractal dimension was calculated by using the box counting method. Thereby, the volume that was analyzed by FIB (Fig. 1) is covered with an orthogonal cubic lattice with increasing lattice constant r. For each box size r the number n of boxes containing any part of the pore structure is counted. To obtain a uniformly distributed data density in a logarithmic plot, it was ensured that the box size grows exponentially, resulting in equidistant data points in a logarithmic plot (Dathe et al., 2001). Regarding fractal objects the relation,

(2)is represented by a straight line and D _f is the slope of that line and c is the intercept. D _f was determined by a linear fit and by considering all data points, which yielded D _f  = 2.41 (Fig. 2).

Fig. 2 A graph of log n (number of boxes) against log r (box sizes). The data were obtained by the box counting method (see text).

Using the determined value of D _f in equation (1) a porosity of 12.4 vol.%. was predicted. This value is about 1 vol.% higher when compared to the value that was determined by BET (Tab. 1). This prediction is considered good and is therefore used as further (in addition to the above relation) evidence that the pore structure of sample BDR is approximately self-similar.

Now that evidence has been presented that pore space in OPA is self-similar, we used the same pore structure at different levels of details or resolutions. The procedure for the construction of a pore structure in Opalinus Clay, which contains a large part of the pore-size spectrum is outlined in Figure 3. We used an original and self-similar (see above) as well as FIB derived pore structure with the size of 300 × 300 × 300 voxels with a voxel size of 10 nm as a starting point.

Fig. 3 Illustration of the process for creating the pore structure in the Opalinus Clay.

Firstly, this image stack was converted into one with a size of 600 × 600 × 600 voxels and a voxel size of 5 nm. Secondly, the original image stack was converted into a 300 × 300 × 300 voxel image stack with a voxel size of 5 nm. From the latter, eight pieces were combined into a cubic image stack. To ensure continuity across the image boundaries, the cubes were arranged so that all boundaries are mirror planes. This compiled image stack was combined with the former image stack (green pores in Fig. 3). This can also be described as a superposition of a finer pore structure with coarser pore structures, where the pore structures are topologically identical. This step is supported by the evidence that the pore structure is self-similar. This procedure was continued up to a voxel size of 1.25 nm and an image stack size of 2400 × 2400 × 2400 voxels. The final pore model is large with regard to the size of the image stack and is therefore only of limited use for the calculation of petrophysical properties. With respect to Pore Size Distribution (PSD) and capillary pressure curve, the calculations were performed for four different sub-volumes with the size of 2400 × 2400 × 640 voxels. Two of each of these sub-volumes were used to calculate the properties parallel and perpendicular to the bedding plane, while the squares of the cuboids are either perpendicular or parallel to the bedding plane. Figure 4a shows the calculated Continuous Pore Size Distributions (CPSD) (Münch and Holzer, 2008), which were compared to the ones that were determined by BET. The CPSD indicates the pore volume fraction that can be filled with spheres of a certain radius. For the calculation it is irrelevant whether the pore space is connected or whether the pore bodies are isolated, therefore the CPSD indicates nothing about the interconnectedness of the pores. It can be assumed that the nitrogen gas penetrates almost the entire pore space of the dried sample during the absorption analysis. Furthermore, the absorption process, as far as we know, is not affected by the ink-bottle effect (see below). Therefore, we compare the BET pore size distribution with the CPSD to see if the model structure in terms of pore size distribution is consistent with the reality. It turned out that the model PSD and the BET target PSD (green lines in Fig. 4a) are similar in terms of the pore-size spectrum and porosity (Fig. 4a).

Fig. 4 (a) Different types of Pore Size Distributions (PSD) of samples and model pore structures as indicated. (b) The MIP-PSD was converted to a capillary pressure curve using the relationship between capillary pressure and pore radii. The corresponding pore radii are also shown.

The capillary pressure behavior related to drainage or imbibition cannot be inferred from the CPSD and BET pore size distribution because these processes are related to a pressure-driven intrusion process that has its beginning at an entry surface into the material. Such a process depends on the pore geometry, which follows the entrance surface. In addition, this process is affected by the ink-bottle effect that occurs after a material (e.g. water, air, mercury) is pressed into the pore space. The pore size distribution related to an intrusion process is the so-called MIP-PSD. Thereby, in a laboratory experiment mercury is gradually pressed into the material with increasing pressure. Because the narrowest pore passages control the applied pressure, the resulting MIP-PSD corresponds rather to a bottleneck size distribution than to the real PSD (see Münch and Holzer, 2008). Here, the MIP process was simulated for the model pore structures (Münch and Holzer, 2008) to calculate their MIP-PSD (Fig. 4). The orientation of the intrusion surface was either parallel or perpendicular to the bedding plane. It must be noted that the MIP intrusion process does not affect the whole pore space but only the part connected to the entrance surface. The same applies to the drainage process (see below). Furthermore, the MIP-PSD underestimates the volume fraction of larger pores because the volume of larger pores located downstream of a narrow pore bottleneck is added to the radius of the bottleneck. In Figure 4 it can be seen that this effect is more pronounced when the direction of intrusion is perpendicular to bedding compared to intrusion parallel to bedding. Consequently, the bottlenecks connecting larger pore bodies in the direction perpendicular to the bedding have smaller radii than those parallel to the bedding. Based on the MIP-PSD the corresponding capillary pressure can be estimated. Based on the MIP-PSD, the capillary pressure P _c can be estimated using the relationship between capillary pressure and pore radius r given as P _c  = 2γcosθ/r, where γ is the interfacial (surface) tension between water and air (72.8 mN/m) and θ is the contact angle. As for the material constants, there are uncertainties regarding the contact angle. Based on the reported values of capillary pressures in OPA (see below), it can be assumed that this clay-rich rock is hydrophilic. In such a case, the contact angle is small (10°–30°) (Borysenko et al., 2009). Thus, we have used 30° for the contact angle in the above equation. A consistent picture emerges that the capillary pressure to be overcome increases more for an intrusion perpendicular to bedding after the entry than for an intrusion parallel to the bedding (Fig. 4b). This is related to the constrictivity anisotropy described above. The largest pores in the entry surfaces vary in the partial volumes, resulting in some variation in the entry pressure (2–5 MPa).

2.3 Calculated two-phase flow properties

The GeoDict software was used to calculate numbers of these properties (Wiegmann et al., 2010). This software uses 3D pore microstructures in voxel representation as input and calculates a variety of flow characteristics by solving flow partial differential equations. GeoDict determines the capillary pressure curve (Fig. 5) on the basis of the pore morphology method (Hipert and Miller, 2001). For the calculation of the relative permeabilities of water and air, the required phase distribution (Fig. 6) was taken from the previously simulated drainage or imbition process related to the capillary pressure curve calculation. After the saturation steps, flow direction and boundary condition have been set, the flow partial differential equations are solved by the software for values of flow properties (Fig. 6) and permeability by using an iterative solver (Wiegmann et al., 2010).

Fig. 5 Capillary pressure curves related to four digital pore structures as indicated (simulations). These curves are compared with laboratory data (WRC, MIP). The latter were subdivided according to depth (see text).

Fig. 6 (a) Stages during drainage of water from the pore structure related to Figures 8a and 8c showing the distribution of water and air related to different water saturations. (b) Visualization of air-flow properties showing the distribution of air flow velocities at different water saturation stages.

The larger pores of the size spectrum within the gas entrance plane control gas entry pressure. For this pore size spectrum there is match between the PSD determined from BET measurements and the PSD related to digital pore structures. Gas entry occurs for gas pressures >~2.4 MPa (see also Sect. 3). Figure 5 shows the capillary pressure curves related to the model pore structures.

In Figure 5, the simulated curves were compared with laboratory data. The relationship between capillary pressure and water saturation can be derived from the Water Retention Curves (WRC) measured in the laboratory (Ferrari and Laloui, 2012; Ferrari et al., 2014; Muñoz et al., 2003; Nagra, 2013; Romero and Gómez, 2013; Zhang and Rothfuchs, 2007). The capillary pressure curve was also evaluated based on data acquired by Mercury Intrusion Porosimetry (MIP) (Romero and Gómez, 2013). The laboratory data were divided into data points referring to samples taken at shallow depths (~300 m) and data points referring to samples taken at greater depths (~800–900 m). The simulations agree reasonably with the laboratory data for flow parallel to bedding and for samples taken from shallower depths. The laboratory data determined for rocks from greater depth show slightly higher capillary pressures. The simulations predict that especially at the beginning of drainage the capillary pressures are higher for a flow perpendicular to the bedding than for a flow parallel to the bedding. This must have to do with the fact that the larger pores are connected by narrow bottlenecks in direction perpendicular to bedding. Therefore, for a given gas pressure the residual saturation depends on the direction in which the gas pressure is applied.

Due to the long calculation time, effective permeabilities related to the drainage and imbibition cycles were calculated only for two partial volumes of the size of 600 × 600 × 600 voxels (Figs. 6–8). Figure 7 shows the corresponding PSD. The pore structure model related to Figures 8b and 8d has the larger pores, which as a consequence leads to a lower air entry pressure at the beginning of the intrusion process. Thereafter, the structural model related to Figures 8a and 8c has about the same capillary pressure because the comparatively smaller pores in the entrance area are connected to the rest of the pore space with similarly large pore constrictions. This shows that not only the bulk constriction distribution determined from the MIP-PSD, where the mercury enters the sample from all sides, but also the directional constriction connectivity plays a role.

Fig. 7 Pore size distributions related to the two partial volumes that were used to calculated the relative permeability shown in Figure 8.

Fig. 8 Capillary pressure curves and related effective permeabilities calculated for two different pore structures: (a) and (c) drainage with subsequent imbibition for the pore structure with lower permeability; (b) and (d) drainage with subsequent imbibition for the pore structure with somewhat higher permeability.

Drainage with subsequent imbibition was calculated for flow parallel to the bedding (Fig. 8). The influence of the anisotropy of the rock on the capillary pressures and the associated permeabilities was studied by calculating drainage with flow parallel and perpendicular to the bedding (Fig. 9).

Fig. 9 (a) Capillary pressure curves and (b) related effective permeabilities for drainage parallel and perpendicular to bedding.

Regarding the used model pore structures, residual saturation of around 0.35 along the drainage curve at P _c  ~ 120 MPa was predicted (Figs. 5 and 8). At this point on the drainage curve water permeability drops to zero and no more water can be displaced. A continuous path of air through the sample starts to form at P _c  ~ 21–25 MPa (gas emergence pressure) in case of the smaller pore models (see also Sect. 3). The related water saturation at gas emergence pressure ranges between 0.8 and 0.93. According to Figures 8c and 8d the formation of continuous gas transport paths in the shaley facies of Opalinus Clay requires the displacement of around 10–20% of the pore water (Figs. 8c and 8d). Afterwards, air permeability increases from zero along the drainage curve to the final air permeability at an air saturation of around 0.65 (Fig. 8c). After that and along the imbibition air permeability curve, air permeability drops to zero at a water saturation in the range of 0.53–0.60. Hence, the volume fraction of trapped air after imbibition is around 0.40–0.47 (Figs. 8c and 8d). The volume fraction of free air during subsequent imbibition is around 0.15–0.25 (Figs. 8c and 8d).

Initial water permeability parallel to the bedding is about a factor of 10 higher than water permeability perpendicular to bedding (Fig. 9b). The same applies to final air permeability at the point on the drainage curve that corresponds to residual water saturation, where no further water is displaced. In Figure 9, a continuous path of air through the samples starts to form at water saturation of around 0.80 regardless whether flow occurs parallel or perpendicular to bedding. However, the corresponding gas emergence pressure is considerably different for the two flow directions. Parallel to bedding gas emergence occurs at around 25 MPa whereas perpendicular to bedding gas emergence occurs at 75 MPa (see also Sect. 3). In summary, water and air permeability behavior during drainage are strongly affected by the anisotropy of the rock.

3 Discussion

The pore size distribution of the constructed model pore structure related to Opalinus Clay agrees with pore size distributions measured in the laboratory. In addition, the calculated capillary pressure curves agree with those determined in the laboratory. This shows that the constructed pore model has properties that not only match the real pore size distribution, but also the crucial bottleneck size distribution. The approach to pore model construction, although quite simple, is apparently effective and the pore models provide plausible predictions. The predicted capillary pressure curves are not necessarily useful when laboratory data is available but the prediction that the relationship between capillary pressure and saturation in Opalinus Clay is directionally dependent has not been considered for Opalinus Clay until now (Nagra, 2013). This is completely different for other geomaterials. For example, for sandstones, the effects of anisotropy on multiphase transport properties have been studied and found to be not negligible (Bakhshian et al., 2020). Since the pore space of Opalinus Clay has a comparatively strong anisotropy (Keller et al., 2011), it is not surprising that the gas transport properties also depend on the anisotropy. In this work the permeabilities related to drainage and imbibition cycles were given as effective values and not as relative as usual. This handling has been proposed by Bear and Breaster (1987) who have shown that the concept of relative permeabilities is only applicable to isotropic materials. Like these authors, we see no advantage in the use of the concept of relative permeability in case of Opalinus Clay.

The simulations predict that gas will enter Opalinus Clay at a pressure of >~2.4 MPa (Tab. 2). These values are consistent with laboratory data related to samples taken from greater depths (Marschall et al., 2005). Furthermore, the predictions agree reasonably with the linear relationship between intrinsic permeability and gas entry pressure established by Davies (1991). For example, the microstructure on which Figures 8a and 8c are based has a permeability of 2.2e-20 m² and an entry pressure of about 6 MPa, which agrees with the laboratory data of Davies (1991) (Tab. 2).

The gas pressure required to form a continuous path for gas transport depends on the model size. The model length in flow direction is the same for all models, but the models have different cross-sectional areas. The probability of the presence of a single connected path across the sample that is formed by pores with larger radii increases with increasing sample size. Therefore, the large pore models (Fig. 5) show an emergence pressure of 10–13 MPa for gas flow parallel to the bedding, while in the small models (Fig. 8) it is a factor of 2 higher (Tab. 2). The large models are therefore considered more realistic. It should also be noted that the emergence pressure of the large models for gas flow perpendicular to the bedding is in the range of 30 MPa (Fig. 5, Tab. 2). The water that has to be displaced in the course of the formation of a continuous gas path varies between 5 and 20 vol.% (Tab. 2). Parallel to bedding, the gas emergence pressure is in the range of 13 MPa and is associated with a pressure limiting bottleneck radius of about 8 nm, while perpendicular to the bedding; the gas emergence pressure is about 30 MPa with an associated pressure limiting bottleneck radius of about 3 nm. Taking into account the stress state in Opalinus Clay, Marschall et al. (2005) concluded that the gas pressure in the Opalinus Clay cannot exceed about 10 MPa without causing mechanical response of the rock. This finding in combination with the results of this study suggests that a gas transport parallel to the bedding might be possible without mechanical response, while a transport perpendicular to the bedding is hardly possible without mechanical response. It should be further noted that during the imbibition following the drainage it is not possible to remove all the gas from the pore space because a fraction of it is unconnected and immobile (Fig. 8).

The material properties described in this study can be used for simulations of gas transport on coarser length scales for example on the mesoscale. On the mesoscale (i.e. mm–cm scale), heterogeneities are related to spatially different contents of non-clayey minerals. Carbonates and quartz are major non-clayey constituents, which can be regarded as non-porous and thus, the porous clay matrix is heterogeneously distributed on various length scales. On the mesoscale the geometry of the clay matrix as well as the spatial distribution of porosity and permeability in the clay matrix control fluid flow. On the mesoscale micropores cannot be considered as physical objects and flow is thus modeled by assuming continuity of fluid phases and Darcy’s law. As shown here, the required material properties for this upscaling such as porosities, capillary pressure and permeabilities can be determined using pore models in combination with computer simulations. Of course, the calculated properties have to be verified experimentally, which is still pending for some of the properties in this study.

4 Conclusion

The present work is based on a combination of microstructural analysis and computer simulations and shows that in case of Opalinus Clay the development of the effective permeability during drainage and the corresponding development of the capillary pressure depend on the transport direction. This material behavior is controlled on the pore scale and in particular by the geometric anisotropy of the pore space. Regarding the suitability of Opalinus Clay for safe containment of nuclear waste, the anisotropic behavior of the effective permeability and capillary pressure during drainage and associated two-phase flow has been neglected until now. Here we show that this neglect is not appropriate because the results of this study show that the two-phase transport in Opalinus Clay and the related material properties such as capillary pressure and effective permeability are anisotropic during drainage.

Acknowledgments

This work was funded by the Swiss National Cooperative for the Disposal of Radioactive Waste (NAGRA) and by the European Joint Program on Radioactive Waste Management (EURAD).

References

All Tables

All Figures

	Fig. 1 Reconstructed microstructures based on FIB. See Table 1 for voxel resolution and porosity. (a) The 3-D reconstruction of the analyzed volume of OPA from the shaley facies based on SE (secondary electron) images. (b) The 3-D reconstruction of the analyzed volume of OPA from the sandy facies based on SE (secondary electron) images.
In the text

	Fig. 2 A graph of log n (number of boxes) against log r (box sizes). The data were obtained by the box counting method (see text).
In the text

	Fig. 3 Illustration of the process for creating the pore structure in the Opalinus Clay.
In the text

	Fig. 4 (a) Different types of Pore Size Distributions (PSD) of samples and model pore structures as indicated. (b) The MIP-PSD was converted to a capillary pressure curve using the relationship between capillary pressure and pore radii. The corresponding pore radii are also shown.
In the text

	Fig. 5 Capillary pressure curves related to four digital pore structures as indicated (simulations). These curves are compared with laboratory data (WRC, MIP). The latter were subdivided according to depth (see text).
In the text

	Fig. 6 (a) Stages during drainage of water from the pore structure related to Figures 8a and 8c showing the distribution of water and air related to different water saturations. (b) Visualization of air-flow properties showing the distribution of air flow velocities at different water saturation stages.
In the text

	Fig. 7 Pore size distributions related to the two partial volumes that were used to calculated the relative permeability shown in Figure 8.
In the text

	Fig. 8 Capillary pressure curves and related effective permeabilities calculated for two different pore structures: (a) and (c) drainage with subsequent imbibition for the pore structure with lower permeability; (b) and (d) drainage with subsequent imbibition for the pore structure with somewhat higher permeability.
In the text

	Fig. 9 (a) Capillary pressure curves and (b) related effective permeabilities for drainage parallel and perpendicular to bedding.
In the text