Boyun Guo, Philip Wortman
Energy Institute of Louisiana, University of Louisiana at Lafayette, Lafayette, USA.
DOI: 10.4236/ojapps.2026.169196   PDF    HTML   XML   6 Downloads   63 Views   Citations

Abstract

Shale permeability is difficult to measure due to shale’s nature of ultra-low transmissibility. It is highly desirable to develop an efficient method for measuring shale permeability. This study proposes a new method for measuring shale permeability based on optical-imaging spontaneous imbibition. The Krüss Drop Shape Analyzer initially developed for measuring the configuration of sessile drop was used for measuring the dynamic spontaneous imbibition in shale cores. The measured dynamic sessile drop contact angle, sessile drop height, and sessile drop wet diameter were used to calculate the dynamic volume and dynamic surface area of the sessile drop. The dynamic surface area data were used for calculating the dynamic fluid evaporation volume. After correction for fluid evaporation, the change in the dynamic volume of the sessile drop was translated into the time dependent depth of spontaneous imbibition based on material balance. A closed form analytical equation was derived in this study to describe the time-dependent depth of spontaneous imbibition. The trend of the spontaneous imbibition profile predicted by this equation is very consistent with the trend of the measured spontaneous imbibition profile with an R-value between 0.9426 and 0.9885 for four US shales. Matching the spontaneous imbibition profile given by the analytical equation to the measured spontaneous imbibition profile allowed for determination of shale permeability. The time required for determining shale permeability depends on the time of imbibition measurement, typically a few hours, mostly for humidity control. This is considered as an efficient method for measuring shale permeability.

Keywords

Unconventional Reservoirs, Shale Formations, Permeability Measurement, Spontaneous Imbibition, Optical Imaging

Share and Cite:

1. Introduction

Permeability of porous media is a proportionality factor in Darcy’s law relating superficial velocity to pressure gradient and fluid viscosity. Permeability is widely used for charactering fluid transmissibility of oil and gas reservoirs, which is vitally important for efficient exploration and enhanced oil and gas production from the reservoirs.

Permeability is essentially dominated by the structure of connected pores in the porous media. For conventional oil and gas reservoirs, their permeabilities in core-scale can be directly measured by injecting fluids, such as nitrogen gas and water, through core samples (core flooding). Permeability can also be estimated using data from capillary pressure and fluid saturation measurements. The relationship between capillary pressure and permeability can be described using Swanson parameter and Capillary Parachor [1]. Permeability in field-scale can be estimated through pressure-transient and/or rate-transient data analyses from well testing.

Characterizing the permeability of tight rocks and shale is challenging due to nano-darcy range flows, high capillary pressures, and strong fluid-rock interactions. The traditional Steady-State Core Flooding method establishes a constant fluid flow rate or a constant pressure differential across a core sample until the inflow rate matches the outflow rate. Permeability is calculated directly using Darcy’s Law once stable conditions are reached. The advantages of the method include 1) direct measurement of permeability defined in Darcy’s law, relying on fewer analytical assumptions and model simplifications compared to transient methods, 2) multi-phase data, allowing for measuring relative permeability for two-phase or three-phase systems, and 3) simulated reservoir conditions, allowing for high-pressure and high-temperature (HPHT) setups that accurately mimic in-situ reservoir stress. The limitations of the method include 1) extremely time-consuming required for achieving a true steady state (weeks or months per sample), 2) equipment limitations, requiring highly precise pumps capable of maintaining ultra-low flow rates and ultra-sensitive differential pressure transducers, and 3) artifact risks due to prolonged exposure to fluids, causing core swelling, clay migration, or chemical degradation before data collection finishes.

Pulse-Decay/Pressure-Pulse Method is a transient technique using a small pressure pulse at the upstream end, and the rate at which the pressure decays upstream and builds up downstream is measured over time. The advantages of this method include 1) rapid testing window, cutting down measurement time from weeks to a few days for low-permeability samples, 2) high accuracy, and 3) minimal rock alteration, minimizing chemical interactions between the fluid and reactive shale clays. The limitations of the method include 1) complex mathematical inversion required for solving non-linear, transient flow equations that are highly sensitive to accurate storage capacity inputs, 2) leak sensitivity where any minor leak in the upstream or downstream system heavily skews the resulting data, and 3) single-phase restriction.

Pressure-Transient Method utilizes crushed core samples placed inside a sealed cell. Gas is rapidly introduced, and permeability is calculated by tracking the pressure decay history as the gas invades the matrix pores of the grains. Its advantages include 1) ultra-fast throughput (minutes to hours), 2) matrix focus, eliminating the artificial connectivity caused by drilling-induced or stress-relief microfractures, and 3) extremely low range of permeability (nano-darcy). This method has limitations including 1) loss of macro-features, 2) lack of confining stress, 3) and grain size dependency.

Pore-Structure-Based Permeability Correlation is an indirect approach to estimating permeability by inserting microstructural data from images or intrusion tests into theoretical models (such as the Kozeny-Carman equation or bundle-of-capillary-tubes models). Input data is typically gathered via Mercury Injection Capillary Pressure (MICP), Scanning Electron Microscopy (SEM)/Focused Ion Beam SEM (FIB-SEM), Micro-Computed Tomography. This method has advantages including 1) rich microstructural insight (direct visual and physical data on pore-throat size distributions, tortuosity, and aspect ratios), 2) non-destructive screening, and 3) anisotropy tracking. Limitations of this method include 1) upscaling errors, 2) MICP damage, and 3) simplified physics.

Spontaneous-Imbibition-Based method monitors the capillary-driven uptake of a fluid into an unconfined or semi-confined core sample over time. By tracking either the weight gain of the sample or the volume of expelled fluid, analytical or numerical models are used to back-calculate the matrix permeability. This method has advantages including 1) inexpensive and accessible, 2) wettability inclusion, capturing fluid-rock interactions, surface tension, and wettability characteristics specific to the rock matrix, and 3) Scalable geometry successfully applied to both regularly shaped core plugs and irregular drill cuttings. The limitations of this method include 1) strong idealizations (equations assume a completely uniform pore structure and constant wettability throughout the process), 2) clay stability issues, and 3) multi-parameter dependency (difficult to isolate permeability from the data, because the imbibition rate is simultaneously controlled by fluid viscosity, relative permeability curves, and capillary pressure).

The recent advancement in high-quality optical imaging and measurements provides a means of monitoring water imbibition in shale formation cores. For example, the Krüss Drop Shape Analyzer can continuously capture the shape data of sessile drop on the surface of shale cores, allowing for determination of time-dependent geometry/configuration of sessile drop [2]. The time-dependent geometry/configuration data can be translated to time-dependent volume of the sessile drop. It is logical to believe that the time-dependent change of the volume of the sessile drop is due to fluid imbibition into the shale and fluid evaporation into the air. If the amount of fluid evaporation into the air can be determined, the amount of fluid imbibition into the shale can be determined. The fluid imbibition rate data, together with a rigorous mathematical model for imbibition, should make it possible to determine shale permeability. This work tested the idea and proved it valid.

2. Theoretical Background

Permeability measures how easily fluids can move through a porous medium. Figure 1 illustrates three types of pores in shale formations. They are a) closed pores, b) dead end pores, and c) passing pores [3]. The first two types of pores do not contribute to shale transmissibility. The fluid transport capacity of a passing pore is theoretically limited by its pore size, i.e., the lower the pore size is, the lower the fluid transport capacity is. Because the pore sizes in different passing pores are different, the distribution of the pore sizes of all passing pores affects permeability. On the other hand, small pore channels promote capillary force that drives the wetting phase to flow through the pores. This means that the fluid transmissibility of porous media is also promoted by the smaller pores.

Figure 1. Closed, dead end, and passing pores in shale formations.

Nelson [4] presented a relationship between permeability, porosity, and pore-throat size for sedimentary rocks. Several modifications were made later using different criteria for pore-throat radius as summarized by Nelson [5]. Nabawy et al. [6] presented a permeability correlation for highly porous and permeable sandstones. Nishiyama and Yokoyama [7] introduced pore tortuosity into the classic correlation. Schulz et al. [8] included the grain sphericity in the permeability, porosity, and pore size correlation. Most of these correlations show that the rock permeability is directly proportional to the square of pore-throat size, or the characteristic pore radius is directly proportional to the square root of rock permeability. The following correlation promoted by Nelson [5] is widely accepted by petroleum engineers:

$k_m=141.2r_m^2{\phi }^2$ (1)

where k_m is measured permeability in md, r_m is the median pore-throat radius from pore-size spectrum in mm, and f is fractional porosity. Because the median pore-throat radius r_m is not measurable for shale reservoir, Equation (1) cannot be used for estimating shale permeability. Fluid imbibition data required.

The following equation was obtained in this study to describe fluid imbibition (see Appendix A for derivation):

$x=\sqrt{\frac{\sigma \cos \theta +0.5r_c\left(P_F-P_R\right)}{245250\frac{{\mu }_w{r}_c}{k_{rw}k}}t}$ (2)

where the imbibition distance x is in cm, liquid-solid interfacial tension s is in Dyne/cm, the contact angle q is in degrees, the equivalent pore radius r_c is in cm, the hydraulic pressure P_F at fluid source is in atm, the reservoir pore pressure at the front of fluid imbibition P_R is in atm, the wetting phase viscosity m_w is in cp, the relative permeability k_rw is dimensionless, the absolute permeability k is in Darcy, and the imbibition time t is in sec. This imbibition equation degenerates to the following equation for spontaneous imbibition test in laboratories where P_F and P_R do not exist:

$x=\sqrt{\frac{\sigma \cos \theta k_{rw}k}{245250{\mu }_w{r}_c}t}$ (3)

This equation is similar to the one given by Handy [9] for spontaneous imbibition when k_rw is equal to the wetting phase saturation (linear relative permeability).

If we assume that the equivalent pore radius r_c can be approximated by the median pore-throat radius r_m, according to Equation (1), the equivalent pore radius r_c can be written as (after unit conversions):

$r_c=\frac{0.000266}{\phi }\sqrt{k}$ (4)

Substituting this relation into Equation (3) gives:

$x=\sqrt{\frac{\phi \sigma \cos \theta k_{rw}\sqrt{k}}{65.24{\mu }_w}t}$ (5)

Because the interfacial tension, contact angle, and fluid viscosity are measurable parameters, if the spontaneous imbibition profile x(t) is accurately measured, matching Equation (5) to the measured imbibition profile allows for determination of the k-value.

3. Workflow

Krüss Drop Shape Analyzer (Figure 2), a high-tech optical goniometer suite, was developed for measuring solid-liquid contact angle. Wortman [2] used the instrument to measure sessile drop geometry parameters, including the time-dependent contact angle q, height H, and wet-diameter 2S on the shale core surface under ambient laboratory conditions (Figure 3). The measurement procedure is outlined as follows.

1) Prepare shale core sample of 1-inch diameter and ½-inch thickness and polish its surface using aluminum oxide sandpaper stepwise up to 2000 grit with a small amount of deionized water (multiple samples were taken over 1.5-to-6-inch areas of the core depending on core intactness, averaging 15 samples per tested region).

2) Place the core sample on the sample platform of Krüss Drop Shape Analyzer with the polished surface facing up.

3) Place a sessile drop of deionized water with a 2 μl droplet volume on the horizontal surface at standard temperature (20˚C) and pressure (1 Atm).

4) Capture and record the sessile droplet image at a rate of 5 frames per second. Measure the time-dependent droplet contact angle(θ), height (H), and diameter (2S) data of the recorded sessile drop.

5) Calculate the time-dependent volume data of the recorded sessile droplet using the following equations (see Appendix B for derivation):

$V_m=\pi a^2\left(\frac{2b^3+D^3}{3b^2}-D\right)$ (6)

where

$b=\frac{HS\tan \left(\theta \right)-H^2}{S\tan \left(\theta \right)-2H}$ (7)

$a=\frac{Sb}{\sqrt{2bH-H^2}}$ (8)

$D=b–H$ (9)

$\tan \left(\theta \right)=\frac{b}{a^2}S{\left(1-\frac{S^2}{a^2}\right)}^{-\frac{1}{2}}$ (10)

Because Equation (7) is not valid for $\theta =\pi /2$ , $b=H$ should be used if $\theta =\pi /2$ .

6) Calculate the time-dependent surface data of the recorded sessile droplet using the following equations (see Appendix B for derivation):

$A=\pi a\left(\frac{b^2}{c}\left({\sinh }^{-1}\frac{c}{b}-{\sinh }^{-1}\frac{c\left(b-H\right)}{b^2}\right)+a-\frac{\left(b-H\right)}{b^2}\sqrt{c^2{\left(b-H\right)}^2+b^4}\right)$ (11)

where

$c=\sqrt{a^2-b^2}$ (12)

7) Calculate the time-dependent evaporation volume of the deionized water using the following equations (see Appendix C for derivation):

$\Delta V_E=-E\sqrt{A}\Delta t$ (13)

where E is solved for a pure evaporation test case on an impermeable plate by minimizing the difference between the measured volume ($$ V_m$$) to the fitted volume ($$ V_f$$) at all timesteps using the Generalized Reduced Gradient (GRG) Method. That is,

$V_{f_{i+1}}=V_{f_i}-E\sqrt{A}\Delta t$ (14)

and

${\displaystyle {\sum }_{i=1}^n\left(V_{m_i}-V_{f_i}\right)}=0$ (15)

where n is the number of time steps.

8) Calculate the time-dependent volume data of the water invasion into the core sample using

$V_{Ii}=V_{m_0}-{\displaystyle {\sum }_{i=1}^n\left(\Delta V_{E_i}\right)}-V_{m_i}$ (16)

9) Calculate the time-dependent Darcy invasion depth data using

$x_i=\frac{V_{Ii}}{\pi \varphi S^2}$ (17)

where f is shale porosity.

Match Equation (5) to the invasion depth data to determine permeability.

Figure 2. Side-view of a Krüss drop shape analyzer 100S.

Figure 3. Image of a typical water sessile drop.

4. Error Analysis

The matching procedure for determining permeability values relies heavily on the accuracy of the contact angle and droplet volume measurements. If the errors in the measurements of contact angle, sessile height, and sessile wet diameter are known from the instrument manufacturer, the error bounds in permeability value can be determined. In fact, the permeability k can be solved from Equation (5) to obtain:

$k={\left(\frac{65.24{\mu }_w{x}^2}{\sigma \cos \theta k_{rw }t}\right)}^2$ (18)

which gives an approximation of partial derivative with respect to q as:

$\frac{\partial k}{\partial \theta }=8512.52{\left(\frac{{\mu }_w{x}^2}{\sigma k_{rw }t}\right)}^2\frac{\text{sin}\left(\theta \right)}{{\cos }^3\left(\theta \right)}$ . (19)

This equation implies that, if the error in q is Dq, the error in k should be:

$\Delta k=8512.52{\left(\frac{{\mu }_w{x}^2}{\sigma k_{rw}t}\right)}^2\frac{\text{sin}\left(\theta \right)}{{\cos }^3\left(\theta \right)}\Delta \theta$ . (20)

Similarly, the error in k due to the errors in H and S can be estimated based on the partial derivatives. Taking a partial derivative of Equation (18) with respect to x gives:

$\frac{\partial k}{\partial x}=4{\left(\frac{65.24{\mu }_w{x}^2}{\sigma \cos \theta k_{rw }t}\right)}^3$ (21)

If the error in x is Dx, the error in k should be:

$\Delta k=4{\left(\frac{65.24{\mu }_w{x}^2}{\sigma \cos \theta k_{rw}t}\right)}^3\Delta x$ (22)

To evaluate Dx, neglecting the change in evaporation volume and substituting Equation (15) into Equation (16) gives:

$x_i=\frac{V_{m_0}-V_{m_i}}{\pi \varphi S^2}$ (23)

which gives

$\frac{\partial x_i}{\partial V_{m_0}}=\frac{1}{\pi \varphi S^2}$ (24)

and

$\frac{\partial x_i}{\partial V_{m_i}}=-\frac{1}{\pi \varphi S^2}$ (25)

which yields

$\Delta x_i=\frac{\Delta V_{m_0}}{\pi \varphi S^2}$ (26)

and

$\Delta x_i=-\frac{\Delta V_{m_i}}{\pi \varphi S^2}$ . (27)

The total error in x_i due to the errors in V_mo and V_mi is $\left(\frac{\Delta V_{m_0}}{\pi \varphi S^2}-\frac{\Delta V_{m_i}}{\pi \varphi S^2}\right)$ . During a short imbibition test, it is reasonable to assume that the total error in the initial sessile volume ($\Delta V_{m_0}$ ) due to the errors in H and S is equal to the total error in the sessile volume at time t_i ($\Delta V_{m_i}$ ) due to the errors in H and S. Therefore, $\Delta x\approx 0$ , and thus, the error in k caused by the errors in H and S is negligible.

5. Case Analysis

5.1. Control Test

Initial tests were conducted on the sample platform of the Krüss DSA 100S which is a flat, painted, aluminum surface which is assumed to be impermeable. A 2 μl droplet stopped spreading on the surface at 13 seconds and began to recede at 510 seconds which marks the beginning and end of the data used to fit the imbibition and evaporation measurements. The 2 μl droplet was measured at the beginning to be 2.03 μl using Equation (6) which verified the calculation model within 1.5% error. The droplet radius, droplet height, and the mean contact angle measured were imported into an Excel file for calculating the volume of droplet at each timestep. The evaporation constant was determined to be 7.6 × 10⁻⁶ cm²/s and the imbibition rate was determined to be 0 cm/s^0.5. The evaporation constant and imbibition rate were used to back calculate and project the volume of the spreading and receding regimes. The average volume difference between the measured and fitted volumes is 0.00012 μl.

5.2. Material Descriptions

The developed method was used to determine permeability of four U.S. shales namely Eagle Ford Shale (EFS), Tuscaloosa Marine Shale (TMS), Marcellus Shale (MS), and Green River Shale (GRS). The EFS has a porosity ranging from 5.30% to 9.79% with an average value of 7.55% [10]. The TMS has a porosity ranging from 3.86% to 9.86% with an average value 6.1% [11]. TMS is characterized by a total clay content of 40 - 80 wt%, quartz of 20 - 40 wt%, less than about 40 wt% of calcite, and a TOC of around 1.0 wt% [12]. The MS has a porosity ranging from 5% to 15% in the southwest region with an average value of 10% [13]. The GRS has a porosity of 10% for lean oil shale with ~1 wt% TOC, typically 5% - 15% carbonate and 20% - 40% illite by weight [14].

Water viscosity and interfacial tension at ambient pressure and temperature are 1.0 cp and 72 dynes/cm, respectively. The relative permeability to water is assumed to be 1.0 in the water-invaded region. The relative permeability is fixed to 1.0 because water saturation is 100% during one-phase imbibition. The end-point relative permeability is assumed to be 1.0 at 100% water saturation. Water contact angles at the surfaces of EFS, TMS, MS, and GRS were measured by the Krüss Drop Shape Analyzer in this study to be 66.68˚, 36.62˚, 52.78˚, and 84.73˚, respectively.

5.3. Eagle Ford Shale (EFS)

Without humidity control, a water droplet deposited on the surface stopped spreading on the surface of EFS at 83 seconds and began to recede at 375 seconds which marks the beginning and end of the data used to fit the imbibition and evaporation measurements. The evaporation constant was determined to be 7.77 × 10⁻⁶ cm²/s and the imbibition rate was determined to be 7.9 × 10⁻⁵ cm/s^0.5. The evaporation constant and imbibition rate were used to back calculate and project the volume of the spreading and receding regimes. The large changes in volume difference between 25 and 30 seconds and 370 and 400 seconds displays how the sensitivity of the diameter measurement affects the results. The precision of the diameter measurement was 0.01 mm, and the measured diameter fluctuates by 0.01 mm from the automatic fitting making small changes between frames during these time periods. For humidity control, a sample of the Eagle Ford Shale was tested in a clear acrylic chamber where the sample was suspended over 100 ml of water on a perforated platform. The chamber was sealed with a perforable plastic covering that was sealed on the corners to prevent airflow from entering or leaving the chamber. The chamber was left to rest for one week to allow the humidity to increase in the chamber. The droplet deposited on the surface stopped spreading on the surface at 140 seconds and began to recede at 396 seconds which marks the beginning and end of the data used to fit the imbibition and evaporation measurements. The evaporation constant was determined to be 3.5 × 10⁻⁶ cm²/s and the imbibition rate was determined to be 1.56 × 10⁻⁴ cm/s^0.5. The evaporation constant and imbibition rate were used to back calculate and project the volume of the spreading and receding regimes. Figure 4 shows the result of match of Equation (5) to the invasion depth data over 3.84 minutes of test time (1,152 points). The final match with R-value of 0.9426 (R² of 0.8885) was achieved with k = 0.00594 md.

5.4. Tuscaloosa Marine Shale (TMS)

The water droplet deposited on the TMS surface stopped spreading at 80 seconds and began to recede at 144 seconds which marks the beginning and end of the data used to fit the imbibition and evaporation measurements. The evaporation constant was determined to be 6.43 × 10⁻⁶ cm²/s and the imbibition rate was determined to be 9.87 × 10⁻⁵ cm/s^0.5. The evaporation constant and imbibition rate were used to back calculate and project the volume of the spreading and receding regimes. Figure 5 shows the result of match of Equation (5) to the invasion depth data over 7.54 minutes of test time (2,262 points). The final match with R-value of 0.9885 (R² = 0.9771) was achieved with k = 0.000205 md.

5.5. Marcellus Shale (MS)

The droplet deposited on the MS surface stopped spreading at 45 seconds and began to recede at 295 seconds which marks the beginning and end of the data used to fit the imbibition and evaporation measurements. The evaporation constant was determined to be 7.86 × 10⁻⁶ cm²/s and the imbibition rate was determined to be 5.52 × 10⁻⁵ cm/s^0.5. The evaporation constant and imbibition rate were used to back calculate and project the volume of the spreading and receding regimes. Figure 6 shows the result of match of Equation (5) to the invasion depth data over 9.74 minutes of test time (2,922 points). The final match with R-value of 0.9711 (R² = 0.9547) was achieved with k = 0.00107 md.

5.6. Green River Shale (GRS)

The droplet deposited on the GRS surface stopped spreading at 113 seconds and began to recede at 788 seconds which marks the beginning and end of the data used to fit the imbibition and evaporation measurements. The evaporation constant was determined to be 7.14 × 10⁻⁶ cm²/s and the imbibition rate was determined to be 9.98 × 10⁻⁶ cm/s^0.5. The evaporation constant and imbibition rate were used to back calculate and project the volume of the spreading and receding regimes. Figure 7 shows the result of match of Equation (5) to the invasion depth data over 16.45 minutes of test time (4,955 points). The final match with R-value of 0.9548 (R² = 0.9116) was achieved with k = 0.000443 md.

Figure 4. Match of model-calculated to the tested water imbibition profiles for a EFS core (test data are from Wortman (2025)).

Figure 5. Match of model-calculated to the tested water imbibition profiles for a TMS core (test data are from Wortman (2025)).

Figure 6. Match of model-calculated to the tested water imbibition profiles for an MS core (test data are from Wortman (2025)).

Figure 7. Match of model-calculated to the tested water imbibition profiles for a GRS core (test data are from Wortman (2025)).

6. Discussion

It is understood that shale permeability can vary substantially with effective stress. Measurements performed at ambient pressure and temperature are therefore interpreted as laboratory-condition values rather than estimates of in-situ reservoir permeability. The absolute permeability values derived from this study under ambient pressure and temperature are useful for shale quality comparison and evaluation purposes. To analyze fluid transport in real reservoir-scale, pressure-transient or rate transient test should be conducted to reveal the effective permeability of shale reservoir under the in-situ stress conditions.

This work proposed a method for estimating shale permeability through inversion of imbibition data with an analytical model. Because the accuracy of the permeability values from the method has not been verified with independent measurements, this method primarily provides an estimation of the actual intrinsic permeability of shale samples.

The new method has some limitations in applications. First, Equation (5) only describes imbibition of one liquid phase in porous media. If more than one liquid phases exist in the imbibition test, the test result cannot be matched with Equation (5) to determine permeability. Secondly, if the imbibition tests are conducted in surface conditions, the determined permeability values are only for comparison of shale transmissibility in surface conditions. However, Equation (5) and its source Equation (2) may still be used for scaling the result from lab condition to the in-situ shale condition.

It is worthy addressing the applicability issue of the method because of gravity. The involved mathematical model, i.e., Equation (5), was derived based on horizontal imbibition where the gravity effect was not considered. The applicability of this equation depends on Bond number which is defined as the ratio of gravitational force to the interfacial force. A high value of the Bond number indicates that the system is relatively unaffected by surface tension effects; a low value (typically less than one) indicates that surface tension dominates [15]. The pore size of shale typically exhibits a bimodal distribution. This means that there are two distinct types of pores present within the shale material: small pores and fractures. The small pores are in the range from 20 to 50 nm, which contributes significantly to the overall pore volume. The fractures have larger openings around 50 µm that also play a crucial role in fluid movement within the shale. For water imbibition in shale cores, it can be shown that the Bond number in the small pores is less than 1 × 10⁻⁸ and the Bond number in the fractures is less than 1 × 10⁻². These low values of Bond number suggest that the shale-water system is essentially unaffected by the gravity effect.

7. Conclusions

Shale permeability is difficult to measure due to the nature of ultra-low transmissibility. This study proposes an efficient method for measuring shale permeability based on optical-imaging spontaneous imbibition. Permeability is estimated by matching a theoretical spontaneous imbibition model to the time-dependent imbibition data recorded in spontaneous imbibition test. The new method was applied to four US shales namely the Eagle Ford Shale (EFS), the Tuscaloosa Marine Shale (TMS), the Marcellus Shale (MS), and the Green River Shale (GRS). The following conclusions are drawn.

1) The Krüss Drop Shape Analyzer initially developed for measuring the configuration of sessile drop can be used for measuring the dynamic spontaneous imbibition in shale cores. The measured dynamic sessile drop contact angle, sessile drop height, and sessile drop wet diameter can be used to calculate the dynamic volume and dynamic surface area of the sessile drop. The dynamic surface area data can be used for calculating the dynamic fluid evaporation volume. After correction for fluid evaporation, the change in the dynamic volume of the sessile drop can be translated into the time dependent depth of spontaneous imbibition based on material balance.

2) A closed form analytical equation was derived in this study to describe the time-dependent depth of spontaneous imbibition. The trend of the spontaneous imbibition profile predicted by this equation is very consistent with the trend of the measured spontaneous imbibition profile, suggesting that the semi-analytical equation reveals the mechanism of spontaneous imbibition.

3) Matching the spontaneous imbibition profile given by analytical equation to the measured spontaneous imbibition profile allowed for estimation of shale permeability. The permeabilities of the EFS. TMS, MS, and GRS were estimated to be 0.00594 md, 0.000205 md, 0.00107 md, and 0.000443 md, respectively.

4) The time required for estimating the shale permeability depends on the time of imbibition measurement, typically a few hours, mostly for humidity control. However, independent measurements are necessary in future studies to validate the accuracy of the new method for permeability determination.

Acknowledgments

The authors are grateful to the Energy Institute of Louisiana at the University of Louisiana at Lafayette for its support to this research.

Author Contributions

Boyun Guo—Mathematical modeling, data analysis, and manuscript writing. Philip Wortman—Experimental measurements.

Appendix A. Analytical Modeling of Imbibition in Porous Media

Horizontal mass transfer in porous media is driven by hydraulic pressure and capillary force and is resisted by viscous friction force. Figure A1 illustrates the fluid between the source plane with pressure P_F and the imbibition front with reservoir pressure P_R.

Figure A1. Sketch to illustrate 1-dimensional fluid flow.

The capillary force acting to the fluid over an area of porous media can be formulated based on interfacial tension, contact angle, perimeter of capillary, and porosity of porous media. It is customary to express the capillary force as a function of capillary pressure through a bulk area:

$F_c=\phi hw{p}_c$ (A1)

where F_c is capillary force in Dyne, $\phi$ is porosity in fraction, h is bulk height in cm, w is bulk width in cm, and p_c in capillary pressure in Dyne/cm². For a capillary having cross-section of equivalent circular shape with an equivalent radius, the capillary pressure is expressed as

$p_c=\frac{2\sigma \cos \theta }{r_c}$ (A2)

where s is interfacial tension in Dyne/cm, q is the contact angle measured in the wetting phase, and r_c is the equivalent capillary radius in cm.

The viscous friction force F_f acting on the porous area over an imbibition depth is expressed as

$F_f=\phi hw{p}_f$ (A3)

where p_f is friction pressure in Dyne/cm². For porous media, the frictional pressure can be expressed by Darcy’s law as

$p_f=981000\frac{{\mu }_{w}vx}{k_{rw}k}$ (A4)

where m_w is wetting fluid viscosity in cp, v is velocity in cm/s, x is penetration of imbibition in cm, k_rw is the relative permeability to the wetting phase, and k is the absolute permeability in Darcy.

For horizontal imbibition processes where gravity effect is zero, applying Newton’s second law of motion to the flowing fluid gives

$\phi hw{P}_F+F_c-F_f-\phi hw{P}_R=\rho \phi hwx\frac{\text{d}v}{\text{d}t}$ (A5)

where r is fluid density in g/cc. Substitutions of Equations (A1) and (A3) into Equation (A5) yield:

$\phi hw\frac{2\sigma \cos \theta }{r_c}-981000\phi hw\frac{{\mu }_{w}vx}{k_{rw}k}+\phi hw\left(P_F-P_R\right)=\rho \phi hwx\frac{\text{d}v}{\text{d}t}$ (A6)

which is simplified to yield:

$\frac{2\sigma \cos \theta }{r_c}-981000\frac{{\mu }_{w}vx}{k_{rw}k}+\left(P_F-P_R\right)=\rho x\frac{\text{d}v}{\text{d}t}$ (A7)

Because $v=\frac{ \text{d}x}{\text{d}t}$ , (A7) is rearranged to give

$\frac{2\sigma \cos \theta +r_c\left(P_F-P_R\right)}{\rho r_{c}x}-981000\frac{{\mu }_w}{k_{rw}k\rho }\frac{ \text{d}x}{\text{d}t}=\frac{{\text{d}}^{2}x}{\text{d}{t}^2}$ (A8)

or

$-A\frac{ \text{d}x}{\text{d}t}-\frac{B}{x}=\frac{{\text{d}}^{2}x}{\text{d}{t}^2}$ (A9)

where

$A=981000\frac{{\mu }_w}{k_{rw}k\rho }$ (A10)

$B=-\frac{2\sigma \cos \theta +r_c\left(P_F-P_R\right)}{\rho r_c}$ (A11)

The governing equation Equation (A9) can be solved using the following initial conditions:

$x=x_0$ at t = 0 (A12)

and

$\frac{\text{d}x}{\text{d}t}=v_0$ at t = 0. (A13)

For the governing equation Equation (A9) take the following form:

$\frac{{\text{d}}^{2}x}{\text{d}{t}^2}+A\frac{\text{d}x}{\text{d}t}+\frac{B}{x}=0$ (A15)

This equation can be solved with the initial conditions expressed in Equations (A12) and (A13). The solution has a singularity at the initial point (t = 0). It can be shown that the acceleration effect represented by second order derivative in the governing equation diminishes in a brief period. If the second order derivative term is negligible, the governing equation degenerates to

$\frac{\text{d}x}{\text{d}t}+\frac{B}{Ax}=0$ (A16)

which takes an integration form of

${\displaystyle {\int }_0^{x}x\text{d}x}=-\frac{B}{A}{\displaystyle {\int }_0^t\text{d}t}$ (A17)

which is integrated to yield

$x=\sqrt{-\frac{2B}{A}t}$ . (A18)

Substituting Equations (A10) and (A11) into Equation (A18) results in a solution for imbibition in porous media:

$x=\sqrt{\frac{\sigma \cos \theta +0.5r_c\left(P_F-P_R\right)}{245250\frac{{\mu }_w{r}_c}{k_{rw}k}}t}$ . (A19)

Appendix B. Mathematical Modeling of Sessile Drop Volume and Surface Area

A mathematical model was developed to calculate the volume and surface area of a sessile drop using the output data from the optical measurements. The height (H), radius (S), and contact angle (θ) are measured from the droplet image. They are defined in Figure B1.

Figure B1. Geometry of a sessile drop (shaded area) on a horizontal surface of a shale sample.

The shaded area in Figure B1 represents a sessile drop. If the surface of a sessile drop assumes an elliptical shape described by:

$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$ (B1)

the upper half of ellipse is expressed as:

$y=b{\left(1-\frac{x^2}{a^2}\right)}^{\frac{1}{2}}$ (B2)

which has a derivative of:

$\frac{\text{d}y}{\text{d}x}=-\frac{b}{a^2}x{\left(1-\frac{x^2}{a^2}\right)}^{-\frac{1}{2}}$ (B3)

Therefore, at the fluid contact point (−S, D), the following equation holds:

$\tan \left(\theta \right)=\frac{b}{a^2}S{\left(1-\frac{S^2}{a^2}\right)}^{-\frac{1}{2}}$ (B4)

At the contact point, Equation (B1) degenerates to:

$\frac{S^2}{a^2}+\frac{D^2}{b^2}=1$ (B5)

Because D = b − H, substituting Equation, (B5) into Equation (B4) gives:

$b=\frac{HS\tan \left(\theta \right)-H^2}{S\tan \left(\theta \right)-2H}$ (B6)

Because Equation (B6) is not valid for $\theta =\pi /2$ , b = H should be used at $\theta =\pi /2$ .

Submitting the b-value from Equation (B6) into Equation (B5) gives the following expression for a:

$a=\frac{Sb}{\sqrt{2bH-H^2}}$ (B7)

The volume of the sessile drop is expressed as:

$V_m={\displaystyle {\int }_D^b\pi x^2\text{d}y}$ (B8)

Substituting Equation (B1) into Equation (B8) gives

$V_m={\displaystyle {\int }_D^b\pi a^2\left(1-\frac{y^2}{b^2}\right)\text{d}y}$ (B9)

which is integrated to yield:

$V_m=\pi a^2\left(\frac{2b^3+D^3}{3b^2}-D\right)$ . (B10)

The area of the sessile drop is expressed as:

$A={\displaystyle {\int }_{\left(b-H\right)}^{b}2\pi x\sqrt{1+{x^{\prime }}^2}\text{d}y}$ (B11)

Substituting Equation (B1) into Equation (B11) gives:

$A={\displaystyle {\int }_{\left(b-H\right)}^b\frac{2\pi a}{b^2}\sqrt{b^4-b^2{y}^2+a^2{y}^2}\text{d}y}$ (B12)

Let

$c=\sqrt{a^2-b^2}$ (B13)

Substituting Equation (B13) into Equation (B12) give:

$A={\displaystyle {\int }_{\left(b-H\right)}^b\frac{2\pi a}{b^2}\sqrt{b^4+c^2{y}^2}\text{d}y}$ (B14)

which is integrated to yield:

$A=\pi a\left(\frac{b^2}{c}\left({\sinh }^{-1}\frac{c}{b}-{\sinh }^{-1}\frac{c\left(b-H\right)}{b^2}\right)+a-\frac{\left(b-H\right)}{b^2}\sqrt{c^2{\left(b-H\right)}^2+b^4}\right)$ (B15)

Appendix C. Formulation of Evaporation Volume of Sessile Drop

To calculate the evaporation during the sessile shrinking period, the equation from Hu et al. [16] is used:

$\frac{\text{d}{V}_E}{\text{d}t}=-\sqrt{2\pi }\frac{d\left(1-h\right)C_{sat}\left[T_a\right]}{\rho }\sqrt{A}$ (C1)

where d is the diffusivity of vapor in air, h is the relative humidity of the air, $C_{sat}$ is the saturated vapor concentration, $T_a$ is the ambient air temperature, $\rho$ is the fluid density, and A is the surface area for evaporation.

If we assume that d, h, $C_{sat}\left[T_a\right]$ , $$ \rho $$ in Equation (C1) are all constant, we can simplify it to a single constant:

$E=-\sqrt{2\pi }\frac{d\left(1-h\right)C_{sat}\left[T_a\right]}{\rho }$ (C2)

Substituting Equation (C2) into Equation (C1) gives

$\Delta V_E=-E\sqrt{A}\Delta t$ (C3)

where, for the pure evaporation case, E is solved by minimizing the difference between the measured volume ($V_m$ ) to the fitted volume ($V_f$ ) at all timesteps using the Generalized Reduced Gradient (GRG) Method. That is,

$V_{f_{i+1}}=V_{f_i}-E\sqrt{A}\Delta t$ (C4)

and

${\displaystyle {\sum }_{i=1}^n\left(V_{m_i}-V_{f_i}\right)}=0$ (C5)

where n is the number of time-steps.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Guo, B., Ghalambor, A. and Duan, S. (2004) Correlation between Sandstone Permeability and Capillary Pressure Curves. Journal of Petroleum Science and Engineering, 43, 239-246.[CrossRef]
[2]	Wortman, P.B. (2025) Well Log Correlation with Experimental Sessile Drop Contact Angle Measurements in the Tuscaloosa Marine Shale. Ph.D. Dissertation, University of Louisiana at Lafayette.
[3]	Zhang, R., Jiang, S., Zhang, L., Wang, H., Zhang, T., Yu, R., et al. (2023) Microscopic Pore Structures and Their Controlling Factors of the Lower Carboniferous Luzhai Shale in Guizhong Depression, China. Geofluids, 2023, 1-18.[CrossRef]
[4]	Nelson, P.H. (1994) Permeability-Porosity Relationships in Sedimentary Rocks. The Log Analyst, 35, 38-62.
[5]	Nelson, P.H. (2005) Permeability, Porosity, and Pore-Throat Size—A Three-Dimensional Perspective. Petrophysics, 46, 452-455.
[6]	Nabawy, B.S. (2013) Estimating Porosity and Permeability Using Digital Image Analysis (DIA) Technique for Highly Porous Sandstones. Arabian Journal of Geosciences, 7, 889-898.[CrossRef]
[7]	Nishiyama, N. and Yokoyama, T. (2017) Permeability of Porous Media: Role of the Critical Pore Size. Journal of Geophysical Research: Solid Earth, 122, 6955-6971.[CrossRef]
[8]	Schulz, R., Ray, N., Zech, S., Rupp, A. and Knabner, P. (2019) Beyond Kozeny-Carman: Predicting the Permeability in Porous Media. Transport in Porous Media, 130, 487-512.[CrossRef]
[9]	Handy, L.L. (1960) Determination of Effective Capillary Pressures for Porous Media from Imbibition Data. Transactions of the AIME, 219, 75-80.[CrossRef]
[10]	Finger, L.J., Godet, A. and Billingsley, L. (2017) Porosity Evaluation in the Lower Eagle Ford Shale Formation, Atascosa County, South Texas. Gulf Coast Association of Geological Societies Transactions, 67, 587.
[11]	Lohr, C.D. and Hackley, P.C. (2018) Using Mercury Injection Pressure Analyses to Estimate Sealing Capacity of the Tuscaloosa Marine Shale in Mississippi, USA: Implications for Carbon Dioxide Sequestration. International Journal of Greenhouse Gas Control, 78, 375-387.[CrossRef]
[12]	Borrok, D.M., Yang, W., Wei, M. and Mokhtari, M. (2019) Heterogeneity of the Mineralogy and Organic Content of the Tuscaloosa Marine Shale. Marine and Petroleum Geology, 109, 717-731.[CrossRef]
[13]	Zamirian, M., Aminian, K. and Ameri, S. (2016) Measuring Marcellus Shale Petrophysical Properties. SPE Western Regional Meeting, Anchorage, 23-26 May 2016, SPE-180366-MS.[CrossRef]
[14]	Burnham, A.K. (2017) Porosity and Permeability of Green River Oil Shale and Their Changes during Retorting. Fuel, 203, 208-213.[CrossRef]
[15]	Li, S., Liu, M., Hanaor, D. and Gan, Y. (2018) Dynamics of Viscous Entrapped Saturated Zones in Partially Wetted Porous Media. Transport in Porous Media, 125, 193-210.[CrossRef]
[16]	Hu, D., Wu, H. and Liu, Z. (2014) Effect of Liquid-Vapor Interface Area on the Evaporation Rate of Small Sessile Droplets. International Journal of Thermal Sciences, 84, 300-308.[CrossRef]