A non-invasive method to directly quantify surface heterogeneity of porous materials

It is extremely challenging to measure the variation of pore surface properties in complex porous systems even though many porous materials have widely differing pore surface properties at microscopic levels. The surface heterogeneity results in different adsorption/desorption behaviors and storage capacity of guest molecules in pores. Built upon the conventional Porod’s law scattering theory applicable mainly to porous materials with relatively homogeneous matrices, here we develop a generalized Porod’s scattering law method (GPSLM) to study heterogeneous porous materials and directly obtain the variation of scattering length density (SLD) of pore surfaces. As SLD is a function of the chemical formula and density of the matrix, the non-invasive GPSLM provides a way to probe surface compositional heterogeneity, and can be applied to a wide range of heterogeneous materials especially, but not limited to, porous media and colloids, using either neutron or X-ray scattering techniques.

S urface heterogeneity is ubiquitous in both natural and manmade materials. It represents the coexistence of different chemical and structural properties on the surfaces of a system. The surface heterogeneity significantly influences interactions 1 , mechanical properties 2,3 , surface reaction rate 4 , and storage and transport phenomena 5,6 of the materials. Tuning surface heterogeneity can therefore have wide applications, including in the pharmaceutical industry [7][8][9] , catalysis 10,11 , and microfluids 5,6 . Moreover, heterogeneous porous materials have major economic impacts in hydrocarbon extraction, water exploitation, and radioactive waste disposal. The local matrix heterogeneity of coalbed methane and shale reservoirs has been shown to have a large impact on gas adsorption and transportation in pores imbedded in the matrix which determine the ability of natural gas extraction and greenhouse gas sequestration in these systems 12 . Thus there is a great need to understand and quantify the surface properties of these heterogeneous porous materials.
Various techniques have been used to determine the surface heterogeneity of materials, such as atomic force microscopy (AFM) 13 , scanning electron microscopy 2,3 , energy dispersive spectroscopy (EDS) 2,3 , X-ray photoelectron spectroscopy (XPS) 13 , auger electron spectroscopy (AES) 14 , and secondary ion mass spectrometry (SIMS) 15 . Among these, AFM is useful to obtain surface heterogeneity with the atomic-level resolution. EDS, XPS, AES, and SIMS are powerful to extract the chemical properties on the surface of the materials which directly affect other properties such as biological adhesion 14,15 and mechanical strength 2,3 and are critical for material performance. However, all the techniques mentioned above are, in general, invasive methods and very difficult to be used to probe the surface properties of pores inside porous materials. They usually require careful pre-processing of sample surfaces which may change the original surface properties of the materials. Inverse gas chromatography (IGC) 1,7 is one of the few non-invasive techniques and is useful for characterizing energetic heterogeneity of surfaces, that is, the distribution of surface sites of different energetic levels. However, the results of IGC depend on the probe molecules and assumptions 16,17 . Thus IGC provides information on 'relative' heterogeneity, which can only be used as fingerprints for comparing different materials. Even though the aforementioned surface characterization techniques are powerful and broadly useful, it is still highly desirable to have a new method to characterize the surface heterogeneity that is non-invasive, model-independent, and able to provide compositional properties and statistically reliable values used for comparison between different materials.
In 1951, Porod established the well-known Porod's scattering law for scattering data of materials measured either by X-rays or neutrons 18,19 , which is now widely used for extracting the surface area and the average scattering length density (SLD) of materials in various relatively homogeneous two-phase systems such as porous materials 20 , biological macromolecules 21 , and colloids 22 . SLD is only a function of molecular formula and material density and therefore is an intrinsic property of a material. Despite its wide applications, the conventional Porod's law is, however, intrinsically not applicable to extract the variation of surface properties for heterogeneous systems such as natural rocks, cement pastes, and multi-phase alloys. Owing to the ubiquity of heterogeneity in natural and engineered materials, a new scattering method is needed for obtaining the distribution of the surface properties of the materials.
In this study, we generalize the Porod's scattering law and develop a rigorous new scattering analysis method called generalized Porod's scattering law method (GPSLM) for extracting the surface property variation of heterogeneous systems. Moreover, GPSLM allows to determine the total surface area and surface averaged SLD more accurately than traditional Porod analysis when there is a large variation of pore surface properties. This novel method is non-destructive, model-independent, and applicable to bulk materials. It gives a dimensionless surface heterogeneity parameter that can be used to compare different materials. The obtained surface heterogeneity can be related to the compositional properties of materials through the calculation of SLD.
We first outline our method with the rigorous proof given in Supplementary Note 1. Then, ideal model systems with known heterogeneous surface properties are used to demonstrate the effectiveness of GPSLM. We further apply GPSLM to natural heterogeneous porous materials: kerogens isolated from shale rocks of different maturities. Kerogen is the organic component in shales that cannot be dissolved by any solvents. It is very important to understand the structure, especially the surface properties, of the pores in kerogens because the pores in kerogens are the major storage locations of most produced gas and the pore network controls the matrix flow of gas for production [23][24][25][26] . A strong correlation between the surface heterogeneity and maturity of isolated kerogen is observed using GPSLM for the first time.

Results
Generalized Porod's scattering law method. The conventional Porod's scattering law 18,19 states that the asymptotic term of the coherent scattering intensity, I, for a two-phase system can be formulated as I Q ð Þ ¼ 2πðΔρÞ 2 Q À4 S V at high Q when the interface between these two phases is relatively smooth. Q is the scattering wave vector and Δρ ¼ ρ 1 À ρ 2 is the contrast of the SLD, where ρ 1 and ρ 2 are the SLD of phase 1 and phase 2 in a two-phase system, respectively. S and V are the total surface area and total volume seen by a neutron or X-ray beam, respectively. Many experimental systems show clear Porod's law scattering patterns ðI Q ð Þ / Q À4 Þ 21,27-30 that can be measured by either small-angle neutron scattering (SANS) or small-angle X-ray scattering. Even though the conventional Porod's law has been widely applied to porous systems to determine the average SLD and the total surface area, it cannot provide any information of the variation of SLD in pore surfaces that is closely linked to the variation of compositional properties of pore matrix. In fact, it is still very challenging to obtain the distribution of surface properties of pores for heterogeneous porous systems using many existing techniques. Here, we develop the GPSLM that can be a powerful tool to provide essential surface heterogeneity information for heterogeneous materials.
Without losing the generality, a schematic picture of a heterogeneous porous material is shown in Fig. 1. Matrices with different SLDs are shown in different colors. In addition, there are both intra-matrix and inter-matrix pores with different shapes and sizes in the system that can be filled with probing fluids with the SLD as ρ f . The domain of each matrix is associated with its interface, S, and ρðSÞ is the SLD of the matrix component whose interface is S. Four parameters important for materials with heterogeneous surfaces can be extracted with GPSLM: the total surface area (S T ), the surface area averaged SLD (ρ A ), the surface area averaged second moment of SLD (ρ M 2 ), and the normalized mean square variation of SLD of the matrix (Δ 2 H ). The mathematical definitions of ρ A , ρ 2 M , and Δ 2 H are given as follows: It should be noted that integrals in Eqs. (1)-(3) are evaluated over the interface S, not the material volume V; ρ A , ρ 2 M , and Δ 2 H are surface averaged, not volume averaged parameters.
Δ H is called the normalized surface heterogeneity. It characterizes the degree that the SLD of matrix components along the interfaces (ρ S ð Þ) deviates from the surface averaged SLD (ρ A ). Because SLD is only a function of the chemical formula and material density, Δ H can be linked to the surface variation of chemical properties such as stoichiometry and composition. If the chemical formulae of the matrices are known, Δ H can be used to estimate the density variation of matrices. Furthermore, Δ H is a dimensionless parameter with absolute value and can be used to compare the surface heterogeneity between different materials. Larger Δ H means the surface properties over the entire interfaces in the material are more heterogeneous.
The complete mathematical derivation of GPSLM is provided in Supplementary Note 1. Here we only briefly outline the approach. For a system with heterogeneous matrices, we can write the general equation of the Porod's scattering law as The mean square deviation of SLD (MSD SLD ), or is the generalized Porod's scattering constant and it is a function of the SLD of guest fluid in the pores (ρ f ).
As aforementioned, for many porous systems, there is a Q region in the scattering pattern following the generalized Porod's law as shown in Eq. (4). The scattering patterns at this Q region change when different probing fluids such as liquid or gas are filled into pores. By varying ρ f , the parameters ρ A , ρ 2 M , and Δ H can be very straightforwardly determined. Before filling the pores with probing fluid, the pores are empty and ρ f ¼ 0. When loading fluid, ρ f can be easily tuned by changing either H/D ratio of fluids or pressure (and therefore density) of gases. The scattering contrasts between the fluid and the matrices are changed when tuning ρ f (Eq. (4)) and the scattering intensity of the material loaded with fluid changes accordingly.
The SANS intensity ratio, which is also the C GPS ratio, at the Q range following the Porod's scattering pattern (Eq. (4)) is defined to be a parabolic function of ρ f (Supplementary Note 1). The minimum of IR ρ f À Á as a function of ρ f is denoted as IR min ðρ f ;min Þ, where ρ f ;min is the fluid SLD when IR ρ f À Á reaches the minimum. Based on the proof in Supplementary Note 1, we have In the derivation, we assume that the solid matrix SLD does not change when filling the probing fluid in pores, and the Q region is high enough so that the generalized Porod's scattering is due to most of the pore surfaces in a system. Moreover, we assume all the pores are accessible to the guest fluid. Experimentally, it is very straightforward to determine ρ f ;min and IR min ρ f ;min by simply varying ρ f to find the minimum of IR ρ f À Á . Then by using algebra relations for Eqs. (1)-(6),ρ A , ρ 2 M , Δ H , and S T can be directly calculated for heterogeneous porous materials. We call this new method the GPSLM. It should be noted that even though the average SLD of various materials has been reported in the literature for decades using the scattering data at the Porod's scattering region, the correct physical meaning of ρ f ;min is finally clarified here using the GPSLM, that is, ρ f ;min is the surface averaged, rather than volume averaged, SLD of a system. Moreover, GPSLM can be also applied to some systems with highly correlated structure such as core-shell systems (see Supplementary Note 1 and Supplementary Fig. 1 for details).
Comparing to the traditional Porod's scattering law method, two new parameters, Δ H and ρ 2 M , can be extracted by GPSLM. GPSLM also determines ρ A and S T more accurately when Δ H is large, and can be used to quantitatively evaluate the effect of the surface heterogeneity on ρ A and S T obtained by the traditional Porod analysis which assumes materials are homogeneous (see brief discussion in the Supplementary Note 6). Application to theoretical model systems. In order to test the accuracy of GPSLM, three ideal model systems with known heterogeneous surface properties are constructed with different SLD distribution of the matrices. These model systems are composed of solid spherical particles with different sizes, SLDs, and number densities dispersed in space (details shown in Supplementary Table 1). Pores are formed between spherical particles. The SANS intensity of these systems filled with different pore fluids with various ρ f can be calculated (details described in Supplementary Note 2). Fig 2a plots Table 2). GPSLM correctly extracts the surface averaged SLD and surface heterogeneity of SLD from the heterogeneous materials. Moreover, Δ H obtained by GPSLM is an absolute value which can be directly compared between different materials. Application to shale kerogens. The GPSLM is applied here to study natural porous materials with heterogeneous surface properties: three kerogen samples isolated from shale rocks using acid digestion 31 . Kerogen is the organic component in shale rocks that is not soluble in any organic solvents. It is commonly agreed that the pores inside kerogen are the major storage locations for light hydrocarbons in most shale reservoirs 25,26 . Therefore, understanding the structure, especially the pore surface properties, of the kerogen in shales is essential to understand the total gas reserve and production rate of a shale gas reservoir. The BET surface area ðS BET Þ and pore volume are first characterized by the commonly used isotherm N 2 adsorption at 77 K for the three kerogen samples with different maturities described by vitrinite reflectance (R 0 ) (Supplementary Table 3 and Supplementary Fig. 5). S BET , pore volume, and R 0 are parameters averaged over the whole samples. Kerogen with higher maturity (higher R 0 ) has slightly higher accessible surface area and pore volume (Supplementary Table 3). However, S BET and pore volume are similar among the three kerogens even though their R 0 values are very different. Thus, other quantities are needed for better sample characterization.
To apply the GPSLM, gas is used as the probing fluid as an example and different ρ f is achieved by changing the gas pressure. We conduct SANS on the three kerogen samples loaded with      Figs 3b-d). Details of the kerogens can be found in Supplementary Table 3 deuterated methane CD 4 at different pressures ranging from 0 to 31.1 MPa (Fig. 3 and Supplementary Figs 3b-d). The scattering intensity at the low Q region follows the generalized Porod's scattering with Q −4 power law feature ( Fig. 3 and Supplementary  Fig. 3) and the solid matrix SLD is not expected to be significantly affected by the loading gas (see Supplementary Note 4). Thus the GPSLM can be applied. The density of methane inside kerogen pores at room temperature can be assumed to be bulk methane density in our cases [32][33][34] , i.e. ρ f ¼ ρ CD 4 (details are given in Supplementary Note 4). ρ CD4 is SLD of bulk CD 4 and is calculated and plotted as a function of pressure in Supplementary Figure 2. Separate SANS measurements of these kerogens loaded with helium demonstrate that the structure of solid kerogen matrix does not change in the pressure range being studied (Supplementary Fig. 3a and Supplementary Note 3). Therefore, the change of SANS curves for kerogens loaded with CD 4 (Fig. 3) is due to the density change of CD 4 in the samples as indicated by the Eq. (4). When increasing CD 4 pressure, ρ CD4 increases while the SLDs of the kerogen matrices maintain as constant positive values that depend on the maturity of the kerogen 35 . The SANS intensity, which is proportional to dS, decreases with ρ CD4 first and then increases when ρ CD 4 becomes larger than the surface area averaged SLD of the matrices, ρ A ( Fig. 3 and Supplementary  Fig. 3).
The intensity ratios, which are also the ratios of generalized C ðGPSÞ ðρ f ¼0Þ , obtained at different CD 4 pressures for the three kerogens are plotted as a function of ρ CD 4 (Fig. 4a). C GPS is extracted from fitting the SANS data in the Q range following generalized Porod's scattering using Eq. (4). For less mature samples, i.e. Sample 1 and Sample 2, the pressure range of CD 4 being used (0 to 31.1 MPa) allows IR ρ f À Á to reach the minimum, i.e. IR min ρ f ;min , and eventually IR ρ f À Á increases again when the pressure further increases. For the most mature sample, i.e. Sample 3, the current pressure range is not high enough to reach IR min ρ f ;min (Fig. 4a). However, since IR ρ f À Á measured at the highest pressure (31.1 MPa) for Sample 3 is very close to zero, we expect that this value is close to IR min ρ f ;min and is the upper limit of IR min ρ f ;min . Equation (5) indicates the fluid SLD, ρ f ¼ ρ CD 4 , at minimum IR ρ f À Á is equal to the surface average of kerogen SLD, ρ A . More mature kerogen is shown to have higher ρ A (Fig. 4a) that is consistent with the hydrogen carbon ratio obtained with the Prompt Gamma-ray Activation Analysis (PGAA) shown in Supplementary Table 5 and Supplementary Note 7. This indicates that less hydrogen is in the sample at this stage of maturation, which is consistent with the literature 35 . Using Eq. (6), the normalized surface heterogeneity, Δ H , is calculated for the kerogens and plotted as a function of vitrinite reflectance, R 0 (Fig. 4b). The parameters extracted from GPSLM are listed in Supplementary  Table 4. Our results for the first time demonstrate the direct experimental observation of the decrease in Δ H with the increase in R 0 for the isolated kerogens by acid digestion, suggesting that the surface properties of the shale kerogen become more homogeneous during the maturation.
After obtaining ρ A and Δ H , the total interfacial area, S T , can be easily calculated using Eq. (4) (see details in Supplementary Note 1 and Supplementary Note 5). The specific surface area, i.e. the total interfacial area per mass of dry kerogen, obtained from the GPSLM, S GPS , is compared with the BET surface area measured by isotherm N 2 adsorption at 77 K, S BET (Supplementary Table 3). Both S GPS and S BET have the same trend that kerogen with higher maturity has higher specific surface area. Since the SANS data following the Porod's law region from about 0.012 Å −1 to 0.03 Å −1 are used to calculate S GPS , only the surface area of pores with length scale approximately larger than L ≈ 2π Q ≈ 200 Å is included. Because S GPS is found to be very close to S BET for all the three kerogens (Supplementary Table 3), this indicates that most of the surface area measured by N 2 adsorption is contributed by large pores with size approximately larger than 200 Å and there are very few small pores in these kerogens.
For shales, the surface properties of pores are keys for determining the storage capacity and distribution of hydrocarbons, especially when hydrocarbon mixtures are considered 23 . In the literature, vitrinite reflectance (R 0 ) is often used as an index for quantifying the maturity of kerogen 36 . Other indices such as total organic carbon, hydrogen index, oxygen index, etc., are also used to characterize the quality and types of kerogen 37 . However, these indices are all average values for the samples and the variation of kerogen surface properties is difficult to obtain. The new GPSLM allows one to obtain Δ H (the width of distribution) and therefore provides essential details for the materials.

Discussion
The classic Porod's scattering law has been widely applied for extracting the surface-to-volume ratio in two-phase systems such as porous materials 20 and colloids 22 . The GPSLM developed in this work significantly extends the applications of original Porod's law from simple homogeneous systems to heterogeneous multiple-phase systems. The new GPSLM is based on the widely used contrast variation technique but is able to extract much more information such as the normalized surface heterogeneity (Δ H ), which is the deviation of SLD from the surface averaged SLD (ρ A ). Δ H directly links to the distribution of compositional properties in the materials. GPSLM can also obtain the values of the total surface area (S T ) and ρ A more accurately in systems with a large variation of surface properties than simply using the traditional Porod's analysis by assuming homogeneous materials. For systems with heterogeneous surface properties, the accuracy of S T and ρ A determined by the traditional Porod's analysis depends on the value of Δ H and the way the contrast variation method is conducted, which is discussed in details in Supplementary Note 6 and Supplementary Fig. 4. It is worth noting that the surface averaged properties are more useful than the volume averaged properties in systems involved in reactions on the surface, such as gas adsorption 12,38 , biological target 8,9 , and catalysis 10,11 , because the distribution of the surface properties can significantly influence the performance of materials. By using the novel GPSLM, we experimentally quantify surface heterogeneity of multiple isolated shale kerogen samples, and discover the interesting correlation between surface heterogeneity and kerogen maturity. To the best of our knowledge, this is the first time that quantitative values of the normalized surface heterogeneity, Δ H , are extracted through scattering method.
In addition to heterogeneous porous materials, GPSLM can be easily applied to other systems to extract ρ A , Δ H , and S T of the heterogeneous surface when the guest fluid SLD can be tuned. For colloidal and some biological suspensions where particles are dispersed in continuous medium such as water or other solvents, the fluid SLD can be tuned by mixing hydrogenated and deuterated solvents, such as H 2 O/D 2 O or C 2 H 5 OH/C 2 D 5 OH, with different ratios. For hydrophobic porous materials whose pores are not accessible to H 2 O/D 2 O, or other solids whose surface properties needed to be characterized without introducing liquid solvents, gas loading with different pressures can be used.

Methods
Acid digestion for shale kerogen isolation. Approximately 20 g of each shale sample was ground to pass a 20 mesh sieve. The samples were then extracted with dichloromethane to remove any soluble hydrocarbons. After drying, the samples were treated with hydrochloric acid (HCl) and left to stand in HCl for at least 2 h with occasional stirring. The samples were rinsed four times with water to remove any calcium or magnesium ions that were released by the HCl.
The decarbonated samples were treated with 70% hydrofluoric acid (HF). The acid was added slowly with stirring until there was no more reaction. After the samples were cooled to room temperature, the samples were transferred to a centrifuge tube, centrifuged, and fresh HF was added to tube. Samples were then placed in an ultrasonic and heated for at least 3 h at 50°C. Samples stayed in a tube with occasional shaking for at least 2 days.
Samples were rinsed three times and then concentrated HCl was added and samples were heated for at least 3 h at 50°C. Samples were rinsed twice with water and twice with distilled water and were ready for freeze drying.
Small-angle neutron scattering. SANS measurements were performed at nSoft-10m SANS and NGB-30 m SANS at the National Institute of Standards and Technology (NIST) Center of Neutron Research (NCNR). The incident neutron wavelength, λ, was chosen to be 5 or 6 Å and the sample-to-detector distances, SSDs, were selected to cover a scattering vector (Q) range from 0.0014 to 0.568 Å −1 . All SANS data were corrected for the sample transmission, the background scattering, and the detector sensitivity to obtain the absolute intensity based on a standard procedure described elsewhere 39 .
The kerogen samples were degassed under vacuum for 2 days before the SANS measurement. The degassed samples were loaded into the high pressure (HP) cells in the helium-filled glovebox. The HP cells are made of stainless steel and contain sapphire windows. The HP cells only allow the reliable data up to Q = 0.3 Å −1 . Helium (He) and deuterated methane (CD 4 ) pressure are controlled by a 100HLf hazardous location syringe pump with a gas loading line connecting to the syringe pump. SANS measurement was first conducted on the shale kerogen under vacuum using a turbo pump linked to the syringe pump. He with 31.1 MPa was then loaded into the kerogen for in situ SANS measurement. After the He measurement, the He was pumped out by the turbo pump. CD 4 with different pressures was loaded into the samples to perform in situ SANS study. All the measurements were maintained at 21°C by the cooling bath system. CD 4 was purchased from Cambridge Isotope Laboratories, Inc. Thermodynamic properties of the bulk methane were calculated using NIST standard reference database software REFPROP 40 .
Volumetric isotherm gas adsorption. Isotherm gas sorption measurements were performed on a carefully calibrated and high accuracy Sieverts apparatus under computer control. Instrument and measurement-protocol details have been published elsewhere 41 .
Prompt gamma-ray activation analysis. Cold neutron prompt gamma-ray activation analysis (PGAA) was performed at the NGD beamline at NCNR. The details of the facilities can be found elsewhere 42 . The samples were vacuumed for two days prior to the measurement to remove residual moisture..