1. Introduction

Defined as the void space in the structure of any given material [1,2] porosity is a metric that can impact many physical and chemical aspects of solid materials [3]. Metal Organic Frameworks (MOF) are an example of important porous materials, where porosity is engineered to achieve massive surface areas [3,4], thus offering the opportunity to store or filter out greenhouse gases such as CO₂ [4,5]. Porosity is critical in catalysis, where surface areas are directly related to reaction yields and conversions [6], as well as in separation, where specific pore shapes and sizes in zeolite materials are exploited as molecular sieves [5,7]. Porosity can also point at the origin and formation processes of geological samples such as rocks. Pores (i.e. the voids), are classified based on their nominal average diameter ‘a’ to micropores (a<1μm) and macropores (a>1μm), although the threshold diameter may vary based on agreements [8], and field of applications [9]. Flow behavior and micro porosity in porous materials is equally important since it regulates the transport of nutrients and pollutants in soils, and determines the extraction of oil and natural gas from host rocks and geological matrices [10]. Opposite to the conventional sandstone reservoirs, carbonate reservoirs such as those present in the Middle East contain a distribution of pores sizes (both macro and micro) of which a significant portion can be in the sub-micron regime. These micropores can hold on to hydrocarbons tightly and thus lead to large quantities of bypassed reserves [11–13]. Traditional experimental methods for the determination of porosity include gas pycnometry and sorption [2], mercury intrusion [1,14], optical [15,16] and scanning electron microscopy [17,18] (generally limited to two-dimension-porosity observations). Advanced X-ray [19,20] and neutron scattering [21,22] techniques offer unprecedented spatial and diameter distribution resolution that not only provide porosity profiles, but also other related parameters such as pore coordination number and throat radii. An overview of common commercial techniques is given in the Appendix section. Unfortunately, in addition to cost, the majority of these methods are invasive, hazardous, or incapable of accurately distinguishing pore size distributions. Therefore, in recent years there has been some research toward less invasive and less hazardous non-contact optical methods that can offer an estimation of micro porosity or total porosity [16]. These methods can also suffer from diffraction limit especially for pores smaller than the wavelength.

Terahertz time-domain spectroscopy (THz-TDS) is a relatively young technique that has been improving rapidly over the last three decades [23–26]. The technique gives access to the elusive band between ~100 GHz and ~10 THz (30 μm and 3 mm) which can be used to inspect chemicals (Fig. 1 (a)) or even everyday structures such as that of closed book [24,27]. The spectrometer generates a THz pulse (Fig. 1(b), (c)) with 0.5 ps pulse width using photoconductive switches [28]; the pulse then travels through the sample. Depending on the sample structure and chemical composition the pulse is reflected, scattered, and absorbed by the sample. The scattered or transmitted pulse is then detected by the spectrometer and the Fourier transform of the transmitted pulse shows the absorbed or scattered lines in the spectral band of the pulse. Due to flexibility and safety of THz-TDS, recently there has been a notable interest in estimation of porosity through THz spectroscopy from pulverized samples especially in pharmaceutical industry [29–33]. It has been shown that the THz pulse transmission and THz pulse delay in the porous media can estimate the macro porosity through effective medium theory [29–31] based on conventional Maxwell-Garnett and the Bruggeman formalisms [30]. Unfortunately, the Maxwell-Garnett approximation is only valid for a low-volume fraction of inclusions, and although the Bruggeman formalism can be used to approximate the effective medium for highly porous systems, none of these methods provide any distinction between micro porosity and macro porosity. Recent attempts have proposed using polarization dependant parameters to investigate the directionality of the pores [32].

 

Fig. 1 Beam geometry and specifications. (a) A slice of limestone sedimentary rock with thickness d is illuminated with THz beam. The scattered light (THz pulse) is measured in the THz receiver. (b) Typical THz pulse in time domain and its spectrum (c) in the frequency domain. The green region is where the signal to noise ratio (SNR) is reliable. Between 2 and 2.2 THz the SNR starts to decline rapidly. 2.2-2.5 THz range can be improper for highly absorptive samples. Commonplace commercial THz-TDS systems usually don’t provide reliable SNR above 2.5 THz.

Download Full Size | PDF

In this study, we assess the use of THz scattering of the macropores in dry samples for total porosity estimation through optical diffusion theory approximation. We further exploit strong THz absorption of water to estimate the micro to macro porosity ratio for limestone sedimentary rock samples. Such enabling sensitivity to water is unique to THz domain and it allows us to analyze subwavelength micro porosity that is not visible at other optical frequencies. A comparison of the macro to micro pore ratios from THz scattering to that of conventional mercury intrusion capillary pressure (MICP) for limestone sedimentary rock samples shows that THz scattering and absorption have notable promise for use as a mercury-free, noninvasive porosimetry method.

2. Theoretical framework

In recent years, THz scattering from random media has been a topic of interest for many applications [34,35]. Transmission mode measurements from such media can be used to infer the spatial distribution of the refractive index and provide noninvasive measurements for solids, aqueous solutions, or even solid-liquid mixtures [34,36]. There are numerous approaches that all are different estimations of Radiative Transfer Equations (RTE). Some of the simplifying assumptions which lead to different forms of RTEs are dense media approximation [37], quasi crystalline approximation [38,39], effective media approximation [40], and diffusion theory approximations [41]. Each approximation holds for certain sample types.

To estimate the total porosity, we start with the diffusion theory assumption, which assumes negligible absorption of the sample at the incident wavelength. If a diffusive medium with thickness d and scattering coefficient μ’_s, is illuminated with an intensity I₀(λ), the transmitted intensity, I(λ), after scattering is found as:

In Eq. (2), k is the wave number, g is the scattering anisotropy term, and Q is the scattering efficiency term. These two functions are both composed of infinite series of Hankel and Bessel based functions (see appendix section) that can be estimated numerically. V_p is the particle or pore volume, V_total is the total sample volume under illumination, n_p is the pore refractive index and n_m is the host medium refractive index. It must be noted that volume fraction is the estimate manifestation of total porosity for our samples. This is because our Mie particles are essentially void pores with air refractive index inside a rock host. If we separate f in Eq. (2), we find volume fraction as a function of wave number and pore diameter as in Eq. (3),

The volume fraction that is found through diffusion theory assumption however can have errors due to a) the absorption of the sample, b) irregularities of the pores (i.e. the pores are not spherical) and c) the presence of pores with diameter much smaller (a<0.1λ) than the THz wavelength.

The volume fraction found here approximates the material as a diffusive medium with spherical pores of radii in the Mie scattering regime. Equation (3) indicates an equation set for the two unknown variables f and a. In ideal conditions, by using only two data points at two different wavelengths of the THz pulse Eq. (3) can be accurately solved to find f and a as in Fig. 2, however, as we will explain in the results section, using more data points can result in a better estimation.

 

Fig. 2 Setup and signal specifications. (a) THz setup. Samples are held on top of posts. Inset figure shows the top view of the sample, the Gaussian beam spot has 1.0 mm diameter at FWHM (b) Temporal profile of THz reference pulse before transmitting through the sample in blue and after the transmission through an arbitrary (4.5 mm) dry limestone sample in red. (c) Fourier transform of the THz pulse before the sample (blue) and after the dry sample (red).

Download Full Size | PDF

As noted in Eqs. (1)-(3) there is no means of estimating micro porosity through diffusion theory as this theory is based on Mie scattering for particles or pores larger than 0.1λ which corresponds to a 13 μm minimum radii at 2.5 THz. Unfortunately, micropores have radii profile smaller than 1 μm therefore we should use a different parameter to be able to estimate or study these smaller pores. One dominant parameter that can be exploited is the spectral absorption at THz bandwidth. Since there is negligible absorption by air-filled pores in the rock we have to saturate these pores with a fluid that shows notable absorption in the THz band. Water at liquid or vapor form is known to have strong THz absorption, in fact water vapor absorption lines in THz band are so dominant that normal spectroscopy for dry samples is usually performed in nitrogen-purged chambers. Since liquid H₂O also has a strong broadband absorption [42] in THz range, water profilometry is a well-known tool in THz analysis of different samples [36,43]. As we will explain further in the results section, we use water as an absorption agent to differentiate between micro and macro porosity in the rock.

3. Experimental setup

The efficacy of our technique was tested on two samples of limestone sedimentary rocks. Limestone was obtained from Caldwell, Texas, USA. The rock also known as Indiana limestone is originated during the Mississippian period also known as lower carboniferous (~360-320 million years ago) and consists of about 98% CaCO₃ with trace amounts (less than 2%) of MgCO₃, SiO₂ and, Fe₃O₄. Limestone is an important material [44,45] in oil and gas industry, that is typically heterogeneous in meso-scales. Limestone can form complex porosity networks due to intrinsic solubility of CaCO₃ in acidic conditions. Carbonate reservoirs, including limestone reservoirs, constitute approximately 60% of world oil resources, with additional potential for gas reservoirs [46]. We did two set of experiments one with dry samples and one with hydrated samples. Each set of experiments was done for six thicknesses and two rock samples with different permeabilities (A≈70 mD and B ≈10 mD).

Starting from two cylindrical rock cores with 38.1mm diameter we prepared 12 sample pieces (6 from each type of porosity) with varying thicknesses from 2 mm to 10 mm. To saturate the samples with water we first degassed the distilled water (1000 mL) in a vacuumed desiccator with magnetic stirrer for one day. Then, we submerged the rock samples in the water and kept them under low vacuum (0.22 atm) for three days. This time is necessary to saturate the majority of the micropores with water with minimal invasion to the sample structure [44,45], and may not be compatible with water-sensitive materials. Our optical setup is a conventional THz-time domain spectrometer that is running in transmission mode (Fig. 2(a))

The setup consisted of an Er fiber laser (Menlo Systems), which emits infrared light pulses at a wavelength of 1550 nm with pulse lengths of 90 fs with average power of 120 mW and with a repetition rate of 100 MHz. These laser pulses were split into two arms and guided through an optical fiber to two THz antennas. These antennas were mounted on a probe head that was fixed (Fig. 2(a)). The upper antenna acted as an emitter and the lower one as a detector. The radiated THz beam was guided through the probe head via four high-density polyethylene lenses (only the two inner ones are shown in the schematic figure), thus focusing the beam on a rock. Holding clamps were used to define a fixed measurement spot for every probed rock sample (Fig. 2(a) inset). For dehydration profile measurements, all of the 12 wet rock samples were placed over a translational stage that periodically placed each rock on the focus of the terahertz. The measurements were made for more than 1200 min until the samples were dried, with each measuring cycle lasting about 3 min. This automatization helped us to multiplex the measurements of our samples. The setup was in a nitrogen purged box to avoid air humidity contributions. A typical measured result is shown in Figs. 2(b) and 2(c). The wet rock notably reduces the peak-to-peak amplitude of the pulse and there is also a delay due to higher refractive index of the CaCO₃ (n_CaCO3 = 1.8) rock compared to air (n_air = 1) [47]. The results will be further explained in the following section.

4. Measurements and estimations

Based on Eq. (3) the volume fraction can be found as a function of pore diameter and wavenumber (Fig. 3). This means if we have the measurements of the scattered THz intensity at two different wavelengths we should be able to form an equation set with two equations and two unknowns that would be precisely solvable to find volume fraction f and pore radius a. This is also indicated in the inset curve at Fig. 3(a1). As noted, the two curves cross each other at only one point, which would be the answer to the equation set. The volume fraction curves are dependent upon the refractive index of the medium (as indicated by Eq. (3)), however, this does not mean that a medium with different refractive index but the same porosity would result in a different incorrect porosity estimation. This is because, while the curves would change with changing refractive index, the points where the curve cross (or the solutions) would still be at the same point for media with different refractive indices. This is true in the Mie scattering range for refractive index. For example, if the refractive index of the medium matches that of the pores then this method would not be able to estimate the porosity because the denominator in Eq. (3) will become zero and there will be no more scattering [48]. For most rocks with air-filled pores such refractive index matching conditions are very unlikely [48].

 

Fig. 3 Volume fraction plots. For each sample the volume fraction function f is a 2D function of the pore radii (0.5a) shown by vertical axis and wavelength or frequency shown by horizontal axis. The color bar indicates the value of the volume fraction between 0 and 1. (a1) Mean volume fraction profile for sample A across varying sample thicknesses between 2.76 and 10.55 mm. The inset figure shows two volume fraction sample curves (the values on 2 vertical lines on the 2D color map) (a2) standard deviation of sample A mean volume fraction (b1) Mean volume fraction profile for sample B averaged across varying sample thicknesses between 2.15 and 6.9 mm (b2) standard deviation of sample B mean volume fraction. Darker blue regions indicate the regions where the standard deviation is lower, and therefore, less oscillations is present in the volume fraction value among all samples. These data points are more probable to be close to the solution of the equation set (or crossing points of the volume fraction curves at different wavelengths) as they are showing less variation across all the sample thicknesses and the nearby frequency ranges. The top left corner flat-yellow regions in (a1) and (b1) are close to the poles of Eq. (3) with f >1 that are then thresholded to 1.

Download Full Size | PDF

There are two practical sources of error for the estimated volume fraction: a) the measurements have a spectral noise that is not considered in Eq. (3) and b) the g and Q functions, which are composed of infinite series of Bessel and Hankel functions can be estimated only to a limited number of terms. These two factors induce inaccuracy into the solution; therefore, instead of all the volume fraction curves ($f\left(k,a\right)$) crossing a single solution point $\left(f^{\prime },a\text{'}\right)$ for different wavelengths (or equivalently different k) as the root of equations set, there is a range of solutions (a range of $\left(f^{\prime },a\text{'}\right)$) which tend to aggregate around the real solution. Therefore, to reduce the estimation error we find the crossing points for all the different wavelengths in the THz spectrum and average them within a threshold variance. Additionally, to better estimate the pore radius and volume fraction for each sample type the measured results were averaged across all six different thicknesses for each rock sample. As seen in Figs. 3(a1) and (b1), there is a large range of possible volume fractions and diameters for each sample. Thus, the deviations in results (Figs. 3(a2) and 3(b2)) were used as a metric to find which results varied least across rock samples and different wavelengths. These variations tend to increase as the frequency increases due to reduction of SNR in the higher frequency ranges of the THz band. Variation also increases with increase in radius size due to oscillations caused by the poles in the scattering efficiency terms (see appendix section Eq. (3)-(10)).

Considering the notable (over 2 orders of magnitude) drop in SNR level after 0.6 THz indicated in Fig. 2(c) we applied an upper limit on the frequency range at 0.6 THz, and an upper limit on the volume fraction variance allowed at 0.04. Points suiting these conditions were localized and averaged to estimate the volume fraction and mean pore diameter $\left(f^{\prime },a\text{'}\right)$ for each sample type. Sample A was estimated to have a porosity of 15.3% and sample B a porosity of 32.3% with average pore radius of 255 µm and 332 μm respectively. This result agrees with our expectations as sample B has more macropores than sample A, however, it has a large 12% average error as will be discussed and compared in the discussion section.

The next set of measurements is done with the water-saturated samples to investigate the micro porosity. The results are presented in Fig. 4(a) for sample A and Fig. 4(b) for sample B. Figure 4 shows that the peak-to-peak of the pulse notably changes with the amount of water that is held in the sample. At the start (point P1 at Fig. 4(a)) when the sample is fully hydrated, the measured amplitude is attenuated notably. The transmission starts to rise rapidly as the macropores are drying out up to the point P2. The dynamics of water drying is dependent upon the sample volume, pore connectivity, and on the size of the pores [10,49]. Micropores (a<1μm) are expected to dry at a slower rate than macropores, since capillary forces tend to hold the water [50,51]. Although a detailed study of capillarity in the pore network of the present samples was not performed, based on the rates of drying observed our data suggests that micropores dehydration starts to merge in the profile at P2 [52–54]. Here, this is found when the slope of the curve starts to deviate by ~10% from the tangent line that crosses the inflection point of the macropore dehydration regime (fast drying rate). An example of such tangent line is shown by a dotted pink line for sample A21 in Fig. 4(a). After P2 the sample transitions to dehydration of micropores. This process is much slower than dehydration of macropores, confirming the hypothesis. The horizontal dashed lines are the peak-to-peak of the transmitted E-field for dry samples before water saturation. For example, this level is shown by P4 for sample A21 in Fig. 4(a). As indicated, the curves for saturated samples do not reach the dashed horizontal lines even after 1200 minutes of dehydration (P3); this means that the full dehydration of the samples down to smallest pores can take days or even weeks. Tables 1 and 2 summarize the measured points for different samples. To estimate the macro porosity we have used the net dehydration at the fast regime (P2-P1) divided by the total attenuation induced by the fully saturated sample (P4-P1); for micro porosity estimation we use the residual slow dehydration range (P4-P2) divided by the same total attenuation. We have then averaged the obtained number across all the sample thicknesses for each type of samples.

 

Fig. 4 (a) Dehydration profiles of sample A. The vertical axis shows the peak-to-peak THz E-field amplitude and the horizontal axis shows the time in minutes. The dashed horizontal lines are the peak-to-peak value measured for the dry sample before the saturation. (b) Dehydration profile for sample B.

Download Full Size | PDF

Table 1. Sample A dehydration profile values.

View Table | View all tables in this article

Table 2. Sample B dehydration profile values.

View Table | View all tables in this article

Other than the porosity ratio, the tables show that the rate of the dehydration for the different sample types is notably different. For example, while sample A22 has almost the same thickness as B22 (Fig. 4), the dehydration time is longer which can be interpreted as the presence of more micropores. The drying time per thickness is consistently longer for all A samples. In addition, thicker samples tend to have a longer dehydration time per thickness, which is expected considering that the inner pores require more time to access the outside surface as evaporation happens. Water evaporation is a macroscopic manifestation of a stochastic molecular process, i.e. a random change in state (liquid to vapor) of water molecules. Furthermore, numbers indicate that the drying time increases super-linearly with sample thickness. Such an observation can be explained by considering that the ratio of the surface available for water evaporation to the total volume of water absorbed into the limestone sample, decreases when sample thickness is increased. Therefore, it is expected that the rate of evaporation of water decreases with increase in sample thickness ratio super linearly. It must be mentioned that the radiation power is too small (p<1µW) to induce notable heat and interfere with the natural surface evaporation process. Also, note that since we are using ratio values or normalized values rather than absolute values, e.g. (90% drying time) the effect of exponential Lambertian absorption with thickness is not directly affecting our calculations. The Lambertian attenuation is however noticeable from absolute peak-to-peak values relative to thickness.

An average additive noise of ~ ± 0.01 is present in the measured data. Therefore, since two values need to be subtracted to calculate the numerator (num) and denominator (den) the error follows ${\Delta }^2={\left(num/den\right)}^2\left({\left(\Delta den/den\right)}^2+{\left(\Delta num/num\right)}^2\right)$ and is ~5% in average. The values are certainly lower for thinner samples (~3%) with stronger signals and larger (~7%) for thicker samples with lower SNR.

There are some fluctuations in the measurements of sample A12 and B22, we believe that these occasional fluctuations (happening for few data points every few hours or so) may have been caused by sample slight movements due to stage movements or occasional vibrations of the optical table. It is noteworthy that the absorption value p4-p3 relates to a fraction of deep micropores that can’t access the surface easily and dehydrate with extremely slow rate (days) through water diffusion processes. This value is expected to be affected by the thickness of the samples with larger relative value for thicker samples.

5. Results and evaluation

To assess our results, we compare them with results from other benchmark porosimetry methods. For this purpose, we studied two methods namely the Barrett-Joyner-Halenda (BJH) method and the mercury injection capillary pressure method. BJH uses N₂ surface adsorption to estimate micro porosity. This method is best suited to provide a high-resolution pore diameter profile for pore diameters less than 200 nm (see appendix section). Mercury intrusion is another invasive method that uses mercury uptake by the porous media at discrete pressure steps up to 413,685 kPa (60,000 psia) to determine pore size distributions. Although hazardous and invasive, this method provides both the macro and micro porosity profiles for a broad range of pore diameters and is well suited for comparison with our results. Figure 5 shows the results for our Indiana limestone samples with mercury intrusion porosimetry. As noted, there are two dominant peaks in Indiana limestone pore profile and there is a significant percentage of micropores in each sample. These micropores are specifically important to the oil industry as they are difficult to detect with other porosimetry techniques and tend to form larger percentage of cumulative volume than what was conventionally predicted [55]. Therefore, extracting such information with non-invasive or optical means is of notable value. Table 3 summarizes the results from this work.

 

Fig. 5 Mercury intrusion results for CaCO₃ samples. (a) Pore volume percent versus pore diameter. (b) Cumulative pore volume percentage versus pore diameter.

Download Full Size | PDF

Table 3. Mercury intrusion data summary for our two Indiana limestone samples

View Table | View all tables in this article

As noted from Table 3, sample B not only has larger pores but is also more porous. In Table 3 the first row “Total intrusion volume at 413,685 kPa” is the total sum of all pore volume intruded by the Mercury at max pressure of the test. This may not be directly representative of porosity since there can be fracturing, crumbling, etc of the rock structure at ~413,685 kPa.

The second row (Total pore area at 413,444 kPa) is like the first row except for the surface area, there is the risk of ‘new area’ being created at such high pressures. “Median pore diameter (volume) at 190 kPa” shows average pore diameter at a given pressure per volume of intrusion and the “Median pore diameter (area) at 126,23 kPa and 0.881 m²/g” shows the average pore diameter at a given pressure per surface area of intrusion. “Average pore diameter” is the average pore diameter from the whole spectrum of pressures, this value is therefore very broad, mean oriented value. “Bulk density at 3.6 kPa” is the density as a function of packing and grain density. This parameter is the density of the rock core with the pores included. It takes into account the porosity and mineral density; generally, lower values correlate with more space occupied by air (or less dense minerals). For example, this value is ~2.0 g/mL for porous sandstone. “Apparent (skeletal) density at 413,444 kPa” is grain density, a very important data point that indicates the density of rock matrix without pores. The “Porosity” row shows the porosity that is obtained from the broad range of pressure all the way up to 413,685 kPa.

Table 4 compares our results with mercury intrusion (M.I.) method and indicates the error rate. On a general qualitative level, the results follow the mercury intrusion results for determining the separation between micropores and macropores. The THz scattering clearly distinguishes the more porous rock (sample B) from the less porous rock (sample A) based on the diffusion theory estimation. Quantitatively, the results from diffusion theory tend to be notably off from the actual values of pore diameter and total volume for sample B. This error is generated because the diffusion theory does not take into account the pores smaller than 0.1λ and the absorption of the sample. The diffusion theory therefore cannot give valid results with pore diameter smaller than Mie scattering regime.

Table 4. Comparison of mercury intrusion with THz scattering and dehydration method.

View Table | View all tables in this article

Thus, the diffusion theory estimation models the medium as a host rock with large spherical bubbles inside. For better results a numerical random media model or effective medium theory might be used in future studies [34,35]. The results for dehydration profile, however, is consistent with the macro to micro porosity from mercury injection with about 1% relative error. This shows that THz absorption of water can be used to estimate micro porosity in carbonate rocks. It must be also mentioned that the results from THz dehydration profile have a rather large variance depending on the thickness of the samples. Table 1 and 2, indicates that samples from 4 to 6 mm tend to provide the best results. This thickness range can be used as a guide for future studies. For thicker samples, the SNR is low and for thinner samples the thickness is nonuniform due to larger pores. Finally, in our study we used the peak-to-peak value, which indicates an average across 0.1 to 2.5 THz band. This is solely to enhance the SNR of our metric signal. As will be discussed further, the liquid water absorption profile is rather broad and follows the same dehydration profile as the overall peak-to-peak.

6. Discussion

In our study, we use the peak-to-peak value as an estimate indicator of the amount of water in the sample by assuming a broadband absorption for liquid water, therefore we are assuming that there is a linear relation between the peak-to-peak amplitude of the pulse and water amount. This assumption induces increasingly larger error for thicker and thicker samples. Also, it must be noted that as the water dries out, it also varies the scattering of the pores as some portion of the pores will be filled with air while some others are filled with water. Figure 6 shows that the thicker samples have larger water scattering and absorption which is expected. Figure 6 also shows that the broadband water absorption assumption is correct but there is higher level of absorption at larger frequencies (f>200GHz). Therefore, integration of these higher frequency components may help the accuracy of the results.

 

Fig. 6 Fourier transform of the pulses measured during dehydration profile. (a) Sample set A. (b) Sample set B. Dark areas circled by green dashed curve indicate the contribution of water absorption and scattering. Comparing the two samples indicate that A samples have less water compared to B samples which indicates higher porosity for B. Horizontal axis is time in minutes and vertical axis is frequency in THz.

Download Full Size | PDF

It is noteworthy that while the dehydration profile can indicate the micro to macro porosity ratio for our rock samples with notable microporosity peaks, it does not indicate any information about the geometrical structure of the pores (e.g. anisotropy and network structure), and therefore, leveraging polarization sensitive parameters after transmission of the pulse may seem as an appealing approach. However, in practice for rock core samples such as those used in our study, making samples with thickness less than 2 mm is rather challenging and pulverizing the sample can affect the intrinsic porosity; additionally, the polarization information is most likely lost after notable scattering through thick diffusive sample. Another reason that makes use of effective medium theory challenging for the case of rock samples as opposed to pulverized pharmaceutical samples is the roughness of the surface itself. The roughness of the surface (especially for thin dry samples) can notably affect the delay for THz pulse where in case of water saturated absorption the thickness of the sample is kept rather uniform across the THz beam profile on the surface.

7. Conclusion

We find THz absorption to be a powerful tool to probe micro porosity through water absorption. These pores can be in deep subwavelength region (a<0.01λ). Our results show that the THz temporal dehydration profile follows the benchmark mercury intrusion results within the 5% error induced by measurement noise. This indicates the potential of THz scattering and absorption for porosimetry of carbonate rocks. We found the 4-6mm thickness to be the most suitable thickness for transmission measurements of limestone sedimentary rock samples with 10² to 10³ SNR level (<1μWatt) THz TDS system. On the other hand, based on our study, THz scattering alone only provides a qualitative mean to assess porosity differences in samples through diffusion theory and it fails to provide accurate estimation of total porosity. Our study also indicates that deeper penetration depth of lower frequency ranges of EM waves can be used to assess mechanical features that are in the subwavelength region of the given wavelength. This promotes the use of microwave and far infrared radiations for noncontact porosimetry.

Appendix

In Mie theory, the scattering efficiency and scattering anisotropy is found by solving the Maxwell’s equations in spherical coordinates by considering the boundary conditions at the particle-medium or in our case pore-rock interface. The scattering cross section ${\text{σ}}_{\text{s}}$ and scattering anisotropy $\text{g}$ are found for spherical pores in the Mie scattering regime as a combination of first kind spherical Bessel functions and spherical Hankel functions [56]. The scattering efficiency${\text{ Q}}_{\text{s}}$, ${\text{Q}}_{\text{s}}={\text{σ}}_{\text{s}}/{\text{πr}}^2$ with a being the pore radius (r=0.5a) and the scattering anisotropy $\text{g}$ is expressed as below with $x=\text{kr}$:

(4)$$Q_s=\frac{2}{x^2}{\displaystyle \sum }_{l=1}^{\infty }\left(2l+1\right)\left({\left|a_l\right|}^2+{\left|b_l\right|}^2\right),$$

(5)$$g=\frac{4}{Q_s{x}^2}{\displaystyle \sum }_{l=1}^{\infty }\left[\frac{l \left(l+1\right)}{l+1}\text{Re}\left(a_l{a}_{l+1}+b_l{b}_{l+1}^{*}\right)+\frac{2l+1}{l \left(l+1\right)}\text{Re}\left(a_l{b}_l^{*}\right)\right],$$

here $a_l$ and $b_l$ coefficients are expressed in terms of $x$ and $y=mx$ (m is the refractive index ratio $m=n_p/n_m$) and,

(6)$$a_l= \frac{{\psi }_l^{\text{'}}\left(y\right){\psi }_l\left(x\right)- m{\psi }_l\left(y\right){\psi }_l^{\text{'}}\left(x\right)}{{\psi }_l^{\text{'}}\left(y\right){\xi }_l\left(x\right)- m{\psi }_l\left(y\right){\xi }_l^{\text{'}}\left(x\right)},$$

(7)$$b_l= \frac{m{\psi }_l^{\text{'}}\left(y\right){\psi }_l\left(x\right)- {\psi }_l\left(y\right){\psi }_l^{\text{'}}\left(x\right)}{m{\psi }_l^{\text{'}}\left(y\right){\xi }_l\left(x\right)- {\psi }_l\left(y\right){\xi }_l^{\text{'}}\left(x\right)},$$

(8)$${\psi }_l\left(z\right)= z j_l\left(z\right)={\left(\frac{\pi z}{2}\right)}^{0.5}{J}_{l+1}\left(z\right),$$

(9)$${\chi }_l\left(z\right)= -z y_l\left(z\right)=-{\left(\frac{\pi z}{2}\right)}^{0.5}{Y}_{l+0.5}\left(z\right),$$

(10)$${\xi }_l\left(z\right)={\psi }_l\left(z\right)+i{\chi }_l\left(z\right)=z{h}_l^{\left(2\right)}\left(z\right).$$

l and l + 0.5 indicate the orders, $j_l\left(z\right)$and $y_l\left(z\right)$ are the spherical Bessel function of the first kind and second kind respectively, $J_{l+1}\left(z\right)$and $Y_{l+0.5}\left(z\right)$denote the Bessel function of the first and second kind respectively,$h_l^{\left(2\right)}\left(z\right)$ represent the spherical Hankel function of second kind [48,57]. z is just an arbitrary scalar variable.

An overview of experimental methods for the determination of porosity is given in Table 5. Mercury intrusion method exploits the non-wetting properties of mercury, which will not enter the pores of most substances by capillary action. Rather, mercury requires an external pressure to penetrate into the porosity of a given material. Such a pressure is inversely proportional to pores size: larger macropores require smaller external pressures to be filled, whereas micropores require higher external pressure.

Table 5. Overview of experimental methods for porosity determination

View Table | View all tables in this article

The main drawback is the use of mercury which makes operation in the laboratory complicated. An additional popular method for surface area determination and consequently porosity estimation is the BET method, which is based on the theory of gas adsorption on solid surfaces by Brunauer, Emmet, and Teller (BET) [58]. The theory assumes that i) sites at the surface of the solid with higher energy will be the sites more likely for gas adsorption, and ii) by increasing the pressure of the gas adsorbing on a sample, adsorption will be more likely on sites of previous adsorption, thus generating a molecular multilayer of gas [58]. A mathematical model is then derived to interpret gas adsorption isotherms, to determine the specific surface area (surface area normalized by mass unit) of a given sample. Several mathematical models have been proposed to amend the BET theory, or extending it to the determination of porosity in various samples. The t-plot and the Barrett-Joyner-Halenda (BJH) method are widely adopted [59]. Unlike the Mie theory in optics, the BJH method assumes that pores have a cylindrical shape and that the radius of the pore corresponds to the sum of i) Kelvin radius plus ii) thickness of the film adsorbed on the pore wall [60]. Overall, gas adsorption models work well in porosity ranges where the pores are smaller than ~1µm, which makes them responsive to the Kelvin equation [59]. The derivation and specificities of each model are discussed in details elsewhere [59].

Limestone samples were pulverized using XRD–Mill McCrone microning mill. The micronized sample was then sieved using a 120 µm sieve. Gas adsorption measurements were performed with a Micromeritics ASAP 2020 surface area and porosity analyzer. Before the analysis, samples were degassed overnight at 200 °C. Nitrogen adsorption isotherms for both sample A and sample B are shown in Fig. 7 (a). This figure shows the amount of gas adsorbed as a function of the pressure of the adsorbing gas; p is the pressure of nitrogen gas adsorbing at the surface of the solid powder whereas p₀ is the vapor pressure of nitrogen gas at the temperature of liquefaction. Therefore, p/p₀=1 means that the adsorbing nitrogen gas is liquefying at the surface of the solid. At any given relative pressure (p/p₀), sample A adsorbs more gas than sample B, indicating a higher surface area for the former sample. This is consistent with the mercury intrusion results used in this study. BET specific surface area calculated with the multipoint method are 9.7±0.1 and 12.8±0.1 m² g⁻¹, respectively. Porosity determinations according to the BJH method are shown in Fig. 7 (b). By comparing the pore size distribution determined via gas adsorption (Fig. 7 (b)) with that obtained with mercury porosimetry (see main text) it is clearly shown that pores ~10 nm wide are detected with gas adsorption, but are not seen with the mercury intrusion. Measurements on core samples cannot be directly related to what it is been shown here since the milling technique might introduce unwanted modifications of the original sample porosity.

 

Fig. 7 N₂ adsorption isotherms for a powder of sample A (specific surface area 9.7 m² g⁻¹) and a powder of sample B (specific surface area 12.8 m² g⁻¹). (b) BJH incremental pore volume versus average pore width.

Download Full Size | PDF

Funding

Massachusetts Institute of Technology (MIT) Media Lab Consortia (2746038).

Acknowledgments

The authors acknowledge technical support from T. Dawn Jobe at Aramco Services Company: Aramco Research Centers – Houston, Prof. Keith Nelson and Jiaojian Shi from Massachusetts Institute of Technology. Thanks to Micromeritics Analytical Services, 4356 Communications Drive, Norcross, GA for MICP analysis.

References and links

1. F. Porcheron, P. A. Monson, and M. Thommes, “Modeling mercury porosimetry using statistical mechanics,” Langmuir 20(15), 6482–6489 (2004). [PubMed]  

2. W. C. Conner, J. F. Cevallos-Candau, E. L. Weist, J. Pajares, S. Mendioroz, and A. Cortes, “Characterization of pore structure: porosimetry and sorption,” Langmuir 2(2), 151–154 (1986).

3. H. K. Chae, D. Y. Siberio-Pérez, J. Kim, Y. Go, M. Eddaoudi, A. J. Matzger, M. O’Keeffe, and O. M. Yaghi, “A route to high surface area, porosity and inclusion of large molecules in crystals,” Nature 427(6974), 523–527 (2004). [PubMed]  

4. O. K. Farha, I. Eryazici, N. C. Jeong, B. G. Hauser, C. E. Wilmer, A. A. Sarjeant, R. Q. Snurr, S. T. Nguyen, A. Ö. Yazaydın, and J. T. Hupp, “Metal-organic framework materials with ultrahigh surface areas: is the sky the limit?” J. Am. Chem. Soc. 134(36), 15016–15021 (2012). [PubMed]  

5. A. Samanta, A. Zhao, G. K. H. Shimizu, P. Sarkar, and R. Gupta, “Post-combustion co₂ capture using solid sorbents: a review,” Ind. Eng. Chem. Res. 51(4), 1438–1463 (2012).

6. A. Fischer, J. Jindra, and H. Wendt, “Porosity and catalyst utilization of thin layer cathodes in air operated pem-fuel cells,” J. Appl. Electrochem. 28(3), 277–282 (1998).

7. A. W. Adamson and A. P. Gast, Physical chemistry of surfaces (John Wiley & Sons, Inc., 1967).

8. M. H. Rahman, B. J. Pierson, and W. I. Wan Yusoff, “Classification of microporosity in carbonates: examples from miocene carbonate reservoirs of central luconia, offshore sarawak, malaysia,”in International Petroleum Technology Conference (2011).

9. S. M. Fullmer, S. A. Guidry, J. Gournay, E. Bowlin, G. Ottinger, A. M. Al Neyadi, G. Gupta, B. Gao, and E. Edwards, ” Microporosity: characterization, distribution, and influence on oil recovery,” in International Petroleum Technology Conference (2014).

10. J. D. Jansen, A Systems Description of Flow Through Porous Media (Springer, 2013).

11. R. Stanley, S. Guidry, and H. Al-Ansi, “Microporosity spatial modeling in a giant carbonate reservoir,” in International Petroleum Technology Conference (2015).

12. E. A. Clerke, “Permeability, relative permeability, microscopic displacement efficiency and pore geometry of m_1 bimodal pore systems in arab-d limestone,” SPE J. 14(3), 524–531 (2009).

13. D. L. Cantrell and R. M. Hagerty, “Microporosity in arab formation carbonates, saudi arabia,” GeoArabia 4(2), 129–154 (1999).

14. F. Porcheron, M. Thommes, R. Ahmad, and P. A. Monson, “Mercury porosimetry in mesoporous glasses: a comparison of experiments with results from a molecular model,” Langmuir 23(6), 3372–3380 (2007). [PubMed]  

15. L. Mei, S. Svanberg, and G. Somesfalean, “Combined optical porosimetry and gas absorption spectroscopy in gas-filled porous media using diode-laser-based frequency domain photon migration,” Opt. Express 20(15), 16942 (2012).

16. T. Svensson, E. Alerstam, J. Johansson, and S. Andersson-Engels, “Optical porosimetry and investigations of the porosity experienced by light interacting with porous media,” Opt. Lett. 35(11), 1740–1742 (2010). [PubMed]  

17. S. Eslava, M. R. Baklanov, C. E. Kirschhock, F. Iacopi, S. Aldea, K. Maex, and J. A. Martens, “Characterization of a molecular sieve coating using ellipsometric porosimetry,” Langmuir 23(26), 12811–12816 (2007). [PubMed]  

18. S. H. Kim and C.-C. Chu, “Pore structure analysis of swollen dextran-methacrylate hydrogels by SEM and mercury intrusion porosimetry,” J. Biomed. Mater. Res. 53(3), 258–266 (2000). [PubMed]  

19. C. C. Egger, C. du Fresne, V. I. Raman, V. Schädler, T. Frechen, S. V. Roth, and P. Müller-Buschbaum, “Characterization of highly porous polymeric materials with pore diameters larger than 100 nm by mercury porosimetry and x-ray scattering methods,” Langmuir 24(11), 5877–5887 (2008). [PubMed]  

20. M. Hartmann and A. Vinu, “Mechanical stability and porosity analysis of large-pore sba-15 mesoporous molecular sieves by mercury porosimetry and organics adsorption,” Langmuir 18, 8010–8016 (2002).

21. C. R. Clarkson, M. Freeman, L. He, M. Agamalian, Y. B. Melnichenko, M. Mastalerz, R. M. Bustin, A. P. Radliński, and T. P. Blach, “Characterization of tight gas reservoir pore structure using usans/sans and gas adsorption analysis,” Fuel 95, 371–385 (2012).

22. B. D. Vogt, R. A. Pai, H. J. Lee, R. C. Hedden, C. L. Soles, W. Wu, E. K. Lin, B. J. Bauer, and J. J. Watkins, “Characterization of ordered mesoporous silica films using small-angle neutron scattering and x-ray porosimetry,” Chem. Mater. 17(6), 1398–1408 (2005).

23. H. Pahlevaninezhad, B. Heshmat, and T. E. Darcie, “Advances in terahertz waveguides and sources,” IEEE Photon. J. 3(2), 307–310 (2011).

24. A. Redo-Sanchez, B. Heshmat, A. Aghasi, S. Naqvi, M. Zhang, J. Romberg, and R. Raskar, “Terahertz time-gated spectral imaging for content extraction through layered structures,” Nat. Commun. 7, 12665 (2016). [PubMed]  

25. B. Heshmat, H. Pahlevaninezhad, Y. Pang, M. Masnadi-Shirazi, R. Burton Lewis, T. Tiedje, R. Gordon, and T. E. Darcie, “Nanoplasmonic terahertz photoconductive switch on GaAs,” Nano Lett. 12(12), 6255–6259 (2012). [PubMed]  

26. B. Heshmat, H. Pahlevaninezhad, and T. E. Darcie, “Carbon nanotube-based photoconductive switches for thz detection: an assessment of capabilities and limitations,” IEEE Photon. J. 4(3), 970–985 (2012).

27. A. Aghasi, B. Heshmat, A. Redo-Sanchez, J. Romberg, and R. Raskar, “Sweep distortion removal from terahertz images via blind demodulation,” Optica 3(7), 754 (2016).

28. B. Heshmat, H. Pahlevaninezhad, T. E. Darcie, and C. Papadopoulos, in 2010 IEEE Radar Conference (IEEE, 2010), pp. 1176–1179.

29. K.-E. Peiponen, P. Bawuah, M. Chakraborty, M. Juuti, J. A. Zeitler, and J. Ketolainen, “estimation of young’s modulus of pharmaceutical tablet obtained by terahertz time-delay measurement,” Int. J. Pharm. 489(1-2), 100–105 (2015). [PubMed]  

30. P. Bawuah, A. Pierotic Mendia, P. Silfsten, P. Pääkkönen, T. Ervasti, J. Ketolainen, J. A. Zeitler, and K.-E. Peiponen, “Detection of porosity of pharmaceutical compacts by terahertz radiation transmission and light reflection measurement techniques,” Int. J. Pharm. 465(1-2), 70–76 (2014). [PubMed]  

31. P. Bawuah, T. Ervasti, N. Tan, J. A. Zeitler, J. Ketolainen, and K.-E. Peiponen, “noninvasive porosity measurement of biconvex tablets using terahertz pulses,” Int. J. Pharm. 509(1-2), 439–443 (2016). [PubMed]  

32. D. Markl, P. Wang, C. Ridgway, A.-P. Karttunen, M. Chakraborty, P. Bawuah, P. Pääkkönen, P. Gane, J. Ketolainen, K.-E. Peiponen, and J. A. Zeitler, “characterization of the pore structure of functionalized calcium carbonate tablets by terahertz time-domain spectroscopy and x-ray computed microtomography,” J. Pharm. Sci. 106(6), 1586–1595 (2017). [PubMed]  

33. D. Markl, J. Sauerwein, D. J. Goodwin, S. van den Ban, and J. A. Zeitler, “non-destructive determination of disintegration time and dissolution in immediate release tablets by terahertz transmission measurements,” Pharm. Res. 34(5), 1012–1022 (2017). [PubMed]  

34. L. M. Zurk, B. Orlowski, D. P. Winebrenner, E. I. Thorsos, M. R. Leahy-Hoppa, and L. M. Hayden, “Terahertz scattering from granular material,” J. Opt. Soc. Am. B 24(9), 2238 (2007).

35. K. M. Nam, L. M. Zurk, and S. Schecklman, “Modeling terahertz diffuse scattering from granular media using radiative transfer theory,” Prog. Electromagn. Res. B 38, 205–223 (2012).

36. E. Castro-Camus, M. Palomar, and A. A. Covarrubias, “Leaf water dynamics of arabidopsis thaliana monitored in-vivo using terahertz time-domain spectroscopy,” Sci. Rep. 3, 2910 (2013). [PubMed]  

37. M. Kaushik, B. W. H. Ng, B. M. Fischer, and D. Abbott, “Terahertz scattering by dense media,” Appl. Phys. Lett. 100(24), 241110 (2012).

38. V. V. Varadan and V. K. Varadan, “The quasicrystalline approximation and multiple scattering of waves in random media,” J. Acoust. Soc. Am. 77, S3 (1985).

39. L. Tsang, C. Te Chen, A. T. C. Chang, J. Guo, and K. H. Ding, “Dense media radiative transfer theory based on quasicrystalline approximation with applications to passive microwave remote sensing of snow,” Radio Sci. 35(3), 731–749 (2000).

40. F. Zhou, Y. Bao, W. Cao, C. T. Stuart, J. Gu, W. Zhang, and C. Sun, “Hiding a realistic object using a broadband terahertz invisibility cloak,” Sci. Rep. 1, 78 (2011). [PubMed]  

41. T. Tarvainen, M. Vauhkonen, V. Kolehmainen, S. R. Arridge, and J. P. Kaipio, “Coupled radiative transfer equation and diffusion approximation model for photon migration in turbid medium with low-scattering and non-scattering regions,” Phys. Med. Biol. 50(20), 4913–4930 (2005). [PubMed]  

42. H. J. Liebe, G. A. Hufford, and T. Manabe, “A model for the complex permittivity of water at frequencies below 1 thz,” Int. J. Infrared Millim. Waves 12(7), 659–675 (1991).

43. M. Exter, C. Fattinger, and D. Grischkowsky, “Terahertz time-domain spectroscopy of water vapor,” Opt. Lett. 14(20), 1128–1130 (1989). [PubMed]  

44. O. Gharbi and M. J. Blunt, “The impact of wettability and connectivity on relative permeability in carbonates: a pore network modeling analysis,” Water Resour. Res. 48(12), 1–14 (2012).

45. B. Bijeljic, P. Mostaghimi, and M. J. Blunt, “Insights into non-fickian solute transport in carbonates,” Water Resour. Res. 49(5), 2714–2728 (2013). [PubMed]  

46. C. H. Moore, Carbonate Reservoirs: Porosity, Evolution & Diagenesis in a Sequence Stratigraphic Framework (Elsevier, 2011).

47. H. Zhan, S. Wu, K. Zhao, R. Bao, and L. Xiao, “CaCO₃, its reaction and carbonate rocks: terahertz spectroscopy investigation,” J. Geophys. Eng. 13(5), 768–774 (2016).

48. B. Groever, B. Heshmat, and R. Raskar, “Tyndall windows: tunable scattering of disordered solid–liquid matching mixtures,” ACS Photonics 3(6), 930–935 (2016).

49. J. C. Melrose, “Role of capillary forces in detennining microscopic displacement efficiency for oil recovery by waterflooding,” J. Can. Pet. Technol. 13(4), 9 (1974).

50. R. Chandler, J. Koplik, K. Lerman, and J. F. Willemsen, “Capillary displacement and percolation in porous media,” J. Fluid Mech. 119(1), 249 (1982).

51. J. Van Brakel, “Pore space models for transport phenomena in porous media review and evaluation with special emphasis on capillary liquid transport,” Powder Technol. 11(3), 205–236 (1975).

52. E. Tsotsas, T. Metzger, V. Gnielinski, E. U. Schlünder, E. Tsotsas, T. Metzger, V. Gnielinski, and E. Schlünder, in Ullmann’s Encyclopedia of Industrial Chemistry (Wiley-VCH Verlag GmbH & Co. KGaA, 2010).

53. T. Metzger, A. Irawan, and E. Tsotsas, “Influence of pore structure on drying kinetics: a pore network study,” AIChE J. 53(12), 3029–3041 (2007).

54. M. Prat, “Percolation model of drying under isothermal conditions in porous media,” Int. J. Multiph. Flow 19(4), 691–704 (1993).

55. J. R. Dixon and B. F. Marek, ” The effect of bimodal- pore size distribution on electrical properties of some middle eastern limestones,” in SPE Annual Technical Conference and Exhibition (Society of Petroleum Engineers, 1990), SPE 20601-MS.

56. L. V. Wang and H.-I. Wu, Biomedical Optics (John Wiley & Sons, Inc., Hoboken, 2009).

57. A. J. Cox, A. J. DeWeerd, and J. Linden, “An experiment to measure mie and rayleigh total scattering cross sections,” Am. J. Phys. 70(6), 620 (2002).

58. S. Brunauer, P. H. Emmett, and E. Teller, “Adsorption of gases in multimolecular layers,” J. Am. Chem. Soc. 60(1), 309–319 (1938).

59. S. Lowell and J. E. Shields, Powder Surface Area and Porosity (Springer Science & Business Media, 2013).

60. E. P. Barrett, L. G. Joyner, and P. P. Halenda, “The determination of pore volume and area distributions in porous substances; computations from nitrogen isotherms,” J. Am. Chem. Soc. 73(1), 373–380 (1951).