Predicting Surface Diffusivities of Gas Molecules in Shale

Carbon dioxide injection can be utilised as a means of both enhancing gas recovery from shales and sequestering carbon, and thereby simultaneously addressing the growing worldwide gas demand, as well as the challenge of greenhouse gas emissions. Greater mobility of CO 2 within the shale improves the displacement efficiency of the originally present CH 4 , as well as, increasing the CO 2 penetration of the shale formation. Previous investigations have indicated that surface diffusion is much more significant than the bulk gas transport in shale gas reservoirs because of the larger fraction of adsorbed phase found in the nanopores of shales. The surface diffusivities of CO 2 on different shales, at various temperatures, have been measured. A fractal theory for predicting the Arrhenius parameters of the surface diffusivity of molecules on heterogeneous surfaces has been applied to the surface diffusion of CO 2 in shales. In line with the theory, it was found that both the pre-exponential factor and the activation energy are functions of the surface fractal dimension. Hence, the surface diffusivity, around a monolayer coverage, on shales could be established from an equilibrium gas adsorption isotherm, once the Arrhenius parameters have been calibrated for the specific chemical species. To the best of our knowledge this study is the first to apply the fractal theory and effectively predict, a priori, surface diffusivity parameters for such structurally and chemically heterogeneous natural samples as shales. This theory now enables the optimization of the designs of CO 2 injection in field applications since surface diffusion is of major importance in the apparent permeability, and, thus, in the gas flow mechanisms.


Introduction
Gas shales are an increasingly exploited resource across the world.Gas reserves in unconventional shale are estimated at nearly 719 trillion cubic metres 1 .Gas recovery from shales can be enhanced greatly by injection of carbon dioxide [2][3][4] .CO2 injection to recover methane also has the advantage of simultaneously sequestering carbon, thereby concurrently addressing the issue of increased greenhouse gas emissions from the use of fossil fuels.The displacement efficiency of the CH4 and the carbon sequestration potential is much enhanced by greater mobility of the CO2.Previous work 5 has shown that the mass transport flux in shales is dominated by the surface diffusion mechanism because of the large internal surface area-to-volume ratio of shale rocks.Hence, a better understanding of the structure-transport relationship for surface diffusion in shale rocks will greatly improve the assessment of the production and storage potential of shale gas reservoirs.
In gas-phase mass transfer, surface diffusion probably plays the major role because of the greater amount of adsorbed gas, particularly within the abundant nanopores within the organic matter of shale gas reservoirs.Surface diffusion is a complex physical phenomenon, which is characterized by an activated process 6,7 .It is a physical process that entails random hopping, as the adsorbed particles move between adsorption sites, that requires a minimum activation energy and experiences an activated transition state.During transport, the adsorbed gas is characterised by a large concentration gradient, and the occupation of a large specific surface area gives rise to a large flux 8,9 .In the presence of surface diffusion, the apparent permeability can be ten times higher than when compared to without 10 .The high magnitude of the surface diffusion contribution to overall mass transport, leads to challenges in the prediction of long-term production for shale gas reservoirs since the structure transport relationship for surface diffusion in shales is not well understood 11,12 .Some experimental investigations have also indicated that, in comparison to the bulk gas transport, surface diffusion, is more significant in particular circumstances 13 , as in the case where the pore network is not yet well-developed within shale gas reservoirs 14 .Hence, it reasonable to say that surface diffusion is considered as an essential mechanism for transport in shale gas resources.
Surface diffusivity depends upon the concentration of the adsorbate on the porous solid and the nonlinearity of the isotherm.The particular value of surface diffusivity derived from experimental measurements depends upon how the surface flux is defined, and this is usually achieved by subtraction, whereby the pore diffusion rate is isolated from the overall rate observed.One procedure, for isolating the pore diffusion contribution, is usually conducted by increasing the temperature until it reaches a point whereby the effects of surface diffusion are greatly reduced.Surface diffusion is rendered negligible at high temperatures because the surface diffusion flux is a product of the surface diffusivity and the surface loading.An increase in temperature results in an increase in surface diffusivity but this is, relatively, much less than the decrease of the surface loading; hence, overall, the surface flux decreases 15,16 .However, this particular technique, has several disadvantages.Firstly, it may not be practical to attain the temperature at which surface diffusion is no longer of any matter since this value may be very high.Secondly, when the pore diffusion rate is compared with the surface diffusion rate, it is usually found that the latter only becomes negligible as the temperature is increased if the adsorption isotherm is linear.This problematic high-temperature method can be avoided by using the half-time method that allows for the extraction of surface diffusivity through a simple physical analysis without the need to resort to the use of the mass balance equation, which is always complicated and computationally intensive 17,18 .The half-time is the period that it takes for the quantity of adsorbate on the porous solid to attain half of the equilibrium amount.
A simple diagram of the inverse of half time of adsorption against concentration factor can be plotted 18 .The pore diffusivity can be obtained from the intercept of the anticipated linear plot, while the surface diffusivity can be calculated from the slope.The application of this particular technique to the mass transport of butane within Ajax carbon has been demonstrated.This paper utilized the theory proposed by Do for the derivation of surface and pore diffusivities, to obtain these parameters for highly heterogeneous adsorbents, namely shales.It will be seen that this technique reveals that surface diffusion is the predominant mass transport mechanism for CO2 in a series of Marcellus shales.
The relationship between surface mass transport and the nature of a surface is poorly understood.Further, the surface of shale rocks is particularly complex, possessing both chemical and geometric heterogeneities, given they are composed of a variety of mineral types, including clays, quartz and carbonaceous materials.The use of fractal models to describe the structural heterogeneity in shales has become increasingly common 19,20 , since fractals enable the discernment of hidden patterns in the face of seemingly intractable disorder.Indeed, it has been found that fractals can provide good structural models for a number of different types of shales [21][22][23] .However, the implications for mass transport have been much less studied 23 , and for surface diffusion not at all.In the past, a fractal theory for surface diffusion, applicable to a variety of molecular species, has been found to be successful for predicting surface diffusivity on relatively homogeneous materials like activated carbons, precipitated silica and porous glasses 24 .It was found that both of the Arrhenius parameters, characterising the variation of the diffusivity on these surfaces with temperature, were directly related to the surface fractal dimension and some other structural parameters of the pore network of the material.The fractal theory was also found to directly predict the compensation effect observed experimentally for surface diffusivity 24,25 .It is the purpose of this work to determine whether this fractal theory for surface diffusion can be applied successfully to more heterogeneous materials like shales, and, thereby, offer a way to predict the variation in surface diffusion flux found in a series of different shales.
In this paper, low-pressure gas adsorption isotherms and helium pycnometer experiments will be primarily conducted to investigate pore size distribution and estimate porosity for a series of Marcellus shales from different depths.However, the low-pressure nitrogen adsorption-desorption measurements are also analysed to estimate surface fractal dimension, according to the Frenkel-Halsey-Hill method.In addition, Rock Eval Analysis and Mineral Liberation Analyser experiments are also applied to determine the TOC and mineralogy, respectively.Furthermore, gravimetric experiments are performed to effectively measure the half-time, and, thereby, pore-surface diffusivities could be derived according to Do's theory.Subsequently, gravimetric gas uptake experiments are reported at three different temperatures, in order to estimate the Arrhenius parameters for surface diffusivity on each shale sample.These data will be used to test the applicability of the fractal theory to highly structurally and chemically heterogeneous natural materials.

Fractal theory for surface diffusion
In an activated process of surface diffusion, where the rate of diffusion varies with temperature, the diffusivity can be represented using the Arrhenius expression: where D0 is the pre-exponential factor and ED the activation energy for diffusion, and the surface diffusion proceeds by a series of activated jumps of range λ occurring on a characteristic timescale of τ, that each have Arrhenius dependence such that τo is the pre-exponential factor and Eτ is the activation energy for the correlation time and, similarly, λo is the pre-exponential factor and Eλ is the activation energy for jump length.
The pre-exponential factor in the Arrhenius expression is the entropic term related to the difference in entropy between the initial and final states of a diffusional hop.This entropy change relates to the number of possible configurations of the migrating molecule in the initial and final states.Past studies 26,27 have determined that the pre-exponential factor of the correlation time has an inverse relationship with the number of available sites to which a molecule can hop, where this space was estimated to be within a jump range of R. For a fractally rough surface, the quantity of accessible destinations to which it is feasible for a molecule of linear extent r to hop to is equivalent to the quantity of molecular sized boxes, N, expected to fill the surface within the characteristic upper length scale R: where d is the scaling law exponent referred to as the fractal dimension of the surface 28 .Therefore, (R/r) d is the number of boxes of size r 2 needed to cover an area A within an upper length scale R. By using eq.2, to adapt an expression previously 26 obtained for zeolites to fractal surfaces, Rigby 29 showed that the pre-exponential factor for the correlation time for the motion of a molecule on a surface characterised by a fractal dimension d is given by the expression: where r is the cross-sectional area of the molecule, R is the apparent limiting upper length scale cutoff for the area as the temperature tends to infinity (related to the upper limit of the jump length), and dr and τor are the fractal dimension and pre-exponential factor, respectively, for a reference material.The correlation times in eq. 3 are predicted correctly for the surface diffusion of benzene adsorbed on a variety of silica surfaces 29 .Combining eqs. 1 and 3 means that: Thus, a linear relationship exists between the natural logarithm of pre-exponential factor for the surface diffusivity and the fractal dimension.In deriving eq. 4 it is assumed that the Arrhenius parameters for the jump length are independent of the surface fractal dimension 30 .
For a fractal dimension to be valid it must hold over a wide range of length-scales.In this surface diffusion model the fractal dimension referred to above, which holds at the jump range of the molecule, must also be the same fractal dimension that holds for shorter length-scales, just above that of the single molecule, occupied by its nearest neighbours.Hence, the condition for a valid fractal dimension is also essential to this theory of surface diffusion.The theory behind the derivations below will be given briefly, since it is described in more detail in a previous paper 24 .The activation energy is the enthalpy term and can also be directly related to surface fractal dimension.In prior studies, it was proposed that the total interaction energy, comprising such as the activation energy for the correlation time or the heat of adsorption, originates from the individual contributions from each of the adjacent adsorption sites and from directly beneath the adsorbed molecule.The convolutions of rougher surfaces mean that they have a higher connectivity, which results in more contributions from the nearby sites.On a fractal surface, the particular number of nearby sites is a function of the surface fractal dimension.
In past work [29][30][31] , it was accepted that just closest neighbor interactions were significant.Hence, it has been shown that the surface fractal dimension is related to total interaction energy Ei by the expression: where the subscript i denotes the particular surface of interest, and ESi the contribution from the surface site directly below the adsorbed molecule, εi the contribution from a single neighbouring site, and Rn a characteristic length-scale.In cases where only the interactions from the closest neighbors are deemed significant, then the distance from the middle of a molecule to the furthest edge of an immediately adjacent site can be denoted by Rn.In this particular case, Rn/r then equals 1.5.From eq. 5, it can be observed that a linear relationship exists between Ei and (Rn/r) d .Eqs. 2 and 5 can be combined to give: where (= 2  −  S r +  r ) and (= −  ) are terms composed only of constants.Eq. 6 demonstrates that the activation energy for the surface diffusivity is a linear function of the group (Rn/r) d . .Rigby 25 showed that a compensation effect results when both the natural logarithm of the pre-exponential factor and the activation energy depend on surface fractal dimension as described above, such that: where m and h are constants.The surface fractal dimension factor in eqs. 4 and 6 can be measured independently using a fractal version of the Frenkel-Halsey-Hill (FHH) equation for the analysis of gas adsorption data 32 .

Theory for adsorption on fractal surfaces
The FHH model has been used extensively by researchers to quantitively characterize pore structure of shales based on N2 adsorption isotherms [33][34][35] .The hypothesis used in the determination of the surface fractal dimension by this strategy is described extensively in other texts 32 , and will only be given briefly here.The procedure depends upon an expression involving multilayer adsorption to a fractal surface such that: where V is the equivalent volume of adsorbed gas at equilibrium pressure P, P0 is the saturation pressure and Vm is the equivalent volume of gas in a monolayer.S is a power law exponent dependent on the surface fractal dimension (d), whereas the constant CF is a pre-exponential factor.
Two limiting cases arise 32 : It has been determined that at the beginning of the multilayer build-up, the film-gas interface is subjected to van der Waal's forces, which act between the solid and the gas, thereby causing the film-gas interface to assume the same shape as the surface roughness.The value of the constant S is then given by: However, generally for thicker surface films, the form of the interface is influenced by the gas-liquid surface tension, thereby causing it to migrate away from the surface, which eventually results in the reduction of the upper external surface area of the film.In this case, S is given by:  =  − 3 (10)  It is noted that the heat of adsorption and the surface fractal dimension are independent parameters, since they are obtained from different parts of the isotherm data.The adsorption heat is derived from the sub-monolayer region, while the surface fractal dimension is determined from the multilayer region.

Theory for obtaining surface diffusivity
As mentioned above, Do proposed a theory whereby surface diffusivity may be obtained from gas uptake measurements 18 .The diffusion fluxes of the free (J) and adsorbed (Js) species are expressed as 36 : where εΜ is the voidage of the particle (cc void volume/cc of total particle envelope volume), Cμ is the concentration in the adsorbed phase (mole/cc), C is the bulk fluid concentration (mole/cc), Ds is the surface diffusivity, and Dp is the pore diffusivity.It is noted that the Do 18 analysis uses equations for standard Fickian diffusion, rather than those of anomalous diffusion associated with pore fractals 37 .The use of the standard Fickian equations for data analysis here is justified as follows.The theory presented in Section 2.1 is for the rate-controlling step in surface diffusion, namely the rate of individual molecular jumps, and thus is sensitive to the characteristic length-scales of this process (ie of the order of the molecular jump length).In contrast, the CO2 uptake experiments, by which the surface diffusivity will be measured, have characteristic timescales and length-scales (ie diffusion path length is ~size of shale particles) that are much bigger than the correlation time of individual molecular jumps and the size of pores in shale (<100 nm), respectively.The shale is not a pore fractal over these length-scales (ie >100 nm to 100s microns) so diffusion is not anomalous in uptake experiments.For example, mercury intrusion porosimetry data (not shown) exhibits no macroporosity.Hence, the random walk of the surface diffusing molecules is not anomalous over long times, and normal diffusion equations are applicable.
The fluxes are determined from the overall surface area of the cross-section, while their units are based on the number of moles transported across a unit section within a unit time.The material balance equations can then be written as follows: The symmetry condition always applies at the centre of the particle, while continuity of flux at the exterior boundary of the pellet is maintained, such that: If the pore diffusion is the only diffusion mechanism and the adsorption isotherm is linear, then the half time of adsorption is given by 38 : = 0.03055 (16)   By carrying out numerical simulation of intermediate cases for a Langmuir isotherm (eq.16), Do 17 derived the following expression for the half time for the pore diffusion model with local adsorption equilibrium: where Cμo((mole/volume of solid), is the adsorbed amount in equilibrium with the bulk concentration Co(mole/volume of fluid).The overall adsorption rate is considered to happen faster than pore and surface diffusion.Local adsorption equilibrium is thus attained within a short time.Therefore, the equilibrium data can be described using the local Langmuir isotherm.
where Cμs is its maximum concentration, and b is the Langmuir constant.The justification for the continued use of the Langmuir isotherm, as employed by Do 18 , in the analysis of CO2 uptake data for the fractal shales studied here is as follows.It should be noted that all isotherm model predictions (Langmuir, BET, fractal BET and fractal FHH) look very similar up to around point B (the first knee of the isotherm corresponding roughly to a statistical monolayer coverage) since the influence of multilayer adsorption (which distinguishes the others from Langmuir) only becomes evident beyond that point.
Since we are only considering the surface diffusion of CO2 up to the region of point B, then the Langmuir model is sufficient to characterise the CO2 isotherm in this region.The surface fractality enters the adsorption isotherm models via two separate effects.First, the surface roughness at short length-scales affects the heat of adsorption for individual molecules, which determines such model parameters as the Langmuir and BET constants.Second, multi-layer build-up is affected by the decline in number of adsorption sites in successive adsorbed layers on fractal surfaces, which is the effect also incorporated into the fractal BET and fractal FHH models, but which is only manifest in experimental data well beyond point B. Further, it is noted that if the surface fractal dimension tends towards a value of 3 (as it does in this work), then the overall forms of the fractal BET and fractal FHH isotherms tend towards that of the Langmuir isotherm.Now, allowing for the pore and surface diffusion being in parallel, and that the isotherm is linear, the half time of adsorption is then given by: = 0.03055 (19)   Finally, when the surface diffusion is the controlling mechanism, the half time will be: R 2 = 0.03055 (20)   no matter what the isotherm nonlinearity is, and the Langmuir isotherm can be used.Thus, by combining the behaviour of the half time at various limits (eqs.16-20), the following general equation for the half time for parallel pore and surface diffusion and any nonlinearity of the isotherm is obtained 17 : where: (22)   and ) (23)   This equation suggests that if one plots Ω versus Ξ, a straight line with the slope (1-εM)Ds and the intercept εMDp is expected.Eq.21 is valid for any of the three possible shapes (as expressed via the geometrical parameters α, β, γ) of the particle.The only difference between the three particle shapes is the values of the parameters α, β, and γ given in Table 1.The samples used in this paper have a spherical shape and will be analysed in Section 4.

Materials
Three core shale samples from the Marcellus Formation were obtained during exploratory drilling of a borehole located in the Appalachian basin, Ohio, USA.The core samples were collected from three different depths: 7804-7807ft, 7834-7837ft and 7864-7867ft.Due to commercial confidentiality reasons, a more detailed location cannot be disclosed.The Marcellus Shale lies within a total area of greater than 100,000 miles, and its depth ranges between 4000 and 8500 ft, having an average thickness of 50-200 ft 37 .The formation contains 1500 trillion cubic feet (tcf) of original gas in place (OGIP) and has 141 TCF of technically recoverable gas 39 .
Each shale specimen was ground and sieved using 0.15-0.18mm metal sifters and placed in a drying oven at 110 °C for 16 h to dehydrate.Subsequently, the prepared sample was stored in a desiccator prior to adsorption measurements.The physical parameters of the shale samples are shown in Table 2.

Volumetric Analysis
A Micromeritics 3Flex volumetric analyzer was used to automatically determine the gas sorption isotherms.The sample preparation process involved the crushing of the samples to powder particle sizes between 150 and 180 mesh, followed by degassing, using a VacPrep Degasser, for 16 hours at 110 o C.The adsorbates used were nitrogen (N2) at 77K, and carbon dioxide (CO2) at 273K.These gases were used to determine the overall pore volume from the Micromeritics 3Flex volumetric analyser.2-2.5 grams of shale samples were used for each N2 isotherm, and a filler rod was utilized in all the experiments.1-1.5 grams of each shale sample were used for the CO2 isotherms.The sample tubes containing these shale samples were then submerged in 50% ethylene glycol solution.These solutions were contained in an isothermal controller maintained at 0 o C. The pressure measurement approach was used to measure the equilibrium isotherm of each sample tested; the partial pressure fluctuation was within 0.1%.Non-Local Density Functional Theory (NLDFT) was used to obtain the pore size distribution 40 .This technique was adopted since it is considered that this is the most appropriate for shales because the pore sizes are so small, and various carbon DFT kernels have been successively used previously in experiments involving organic materials like carbons 40,41 .In this paper, the carbon slit pore model of NLDFT kernel was applied for meso-and macro-porosity determination using N2 adsorption data, and the CO2-DFT model was applied for microporosity using CO2 adsorption data.The determination of the bulk density for each sample was performed using helium pycnometry.

Heat of Adsorption
The measurement of the heats of adsorption (such as in Figure 1) has been done using simultaneous calorimetry, as opposed to the method of isosteres and the Clausius-Clapeyron equation.The Clausius-Clapeyron (CAC) method cannot be used because it may lead to an overestimated isosteric adsorption heat because it utilizes the ideal gas law and assumes that the heat does not depend upon temperature.Both of these common assumptions do not apply in this case.
The adsorption was performed in a Quantachrome Autosorb iQ.In order to precisely measure the real heats of adsorption for the duration of equilibration and the adsorption process, a Sensys Evo Tian-Calvet heat flux calorimeter was used.The calorimeter had to be placed as close as possible to the analysis station of the Autosorb iQ, with the aim of minimizing the length of the tube between the calorimeter cell and the analysis station.All the calorimeter results have been pre-calibrated with an indium standard.

Gravimetric Analysis
A Hiden XEMIS (gravimetric analyser) was used to obtain kinetic gas uptake data, and a schematic diagram of this apparatus is shown in Figure 2. In the case of each sample, a series of experiments was conducted with a range of sizes of steps in exterior bulk pressure from vacuum to different ultimate pressures.In the uptake experiments the set of values chosen for ultimate pressures in the pressure steps were close to that required to achieve a statistical monolayer according to the carbon dioxide adsorption isotherms.Even though experiments were carried out at different bulk concentrations over the isotherm, the range of bulk gas concentrations steps was consistent for all shale samples.
Gravimetric measurements were conducted using a sensitive microbalance which measured the change in mass of an adsorbent sample subjected to a step change in adsorbate concentration.This represents a direct indicator for the adsorption rate onto the solid.In the kinetic measurements case, only a small sample (usually of the order of ten milligrams) was used.
The adsorbents to be used were initially degassed overnight to remove the excess moisture.On completion of degassing the adsorbent is brought to the adsorption temperature by immersing in an isothermal water bath.At this point the system is ready for commencing adsorption.At time zero adsorbate gas passes through the ceramic tube (see Figure 2) and the sample weight change was monitored until constant mass was observed indicating that equilibrium had been attained.The adsorbent was then degassed again until constant mass was achieved and the gas that was initially injected had been desorbed.This procedure was repeated for all the points.Between twelve to eighteen data points were normally acquired to characterise an uptake curve.Lastly, data correction for buoyancy effects was also made.

Results and Discussion
In order to characterise the inorganic mineral content of the shale samples, Mineral Liberation Analysis (MLA) was used.The information thereby obtained for the samples is listed in Table 2. From Table 2 it can be seen that the samples are heterogeneous, being predominantly composed of illite, quartz, and carbonaceous phases.It is further noted that no kaolinite or smectite were detected in the Marcellus shale samples, which shows that these particular clay minerals must have undergone a complete transformation.Typically, the overall transformation of the clay takes place in two stages; kaolinite first transforms to smectite, and, then, smectite transforms into illite.The densities for the shale samples were obtained from helium pycnometry, and the values given in Table 2 are similar to those reported previously 42 for Marcellus shale of 2.63 g/cc, and shales in general 43 of 2.06-2.75g/cc.Figures 3 and 4 show the adsorption isotherms for nitrogen and carbon dioxide, respectively, on Marcellus Shale samples from three different depths.It is highlighted that the pore volume of the microporosity decreased with depth, which also resulted in lower surface area (Figure 5).The reason for this decrease in the micropore volume was a combination of the change in content of both TOC and illite, as can be seen in Table 2. shows an example of a fractal FHH plot obtained using one of the sets of nitrogen adsorption isotherm data given in Figure 3.The parameters obtained from such fits (using eqs.8 and 10) to all the nitrogen isotherm data above a statistical monolayer coverage are given in Table 3.It can be seen that the surface fractal dimension tends to decline with increasing depth of the shale sample.Examples of the isotherm data for carbon dioxide on the three shale samples, obtained via the gravimetric method, can be seen in Figure 7(a-c).Since the local adsorption equilibrium is reached quickly, these isotherm data were reasonably fitted to a Langmuir isotherm expression by using a nonlinear regression technique for the selected temperatures of 10 o C,20 o C, and 30 o C. The agreement between the experimental data and Langmuir isotherm model was good, and the parameters, Cμs and b, thereby obtained from this non-linear regression are given in Table 4 for the three temperatures.From Figure 7(a-c), it can be seen that the Langmuir isotherm fits the experimental data well for the three temperatures.Figure 8 shows examples of the fractional uptake of CO2 with time which was obtained via the gravimetric method at various different ultimate bulk concentrations of CO2.When the initial bulk gas concentration is lower, the adsorption equilibrium time will also be shorter.A single exponential Linear Driving Force (LDF) model was used in order to fit the experimental data and thereby obtain the rate constant to find the half time of adsorption 44 .The measured half-times were ultimately used in eq.24 to estimate the surface diffusion.Figure 9(a-c) show the experimental fractional uptakes at three similar bulk concentrations of CO2 for samples of spherical particles of shale having a radius of 0.00825 cm.It is noted that the half-time is decreasing with temperature for all three bulk concentrations and for all samples; that is the time to attain equilibrium is reduced at higher temperatures.These results are in agreement with past findings of Do for a spherical particle 44 .At this point, it should be emphasized that higher gas mobility systems do not necessarily mean that they will reach equilibrium at a much faster rate.The time to reach equilibrium is also dependent upon the ultimate adsorbed quantity that the solid can accommodate at the equilibrium state.The speed with which a system approaches equilibrium is determined by two factors: capacity and mobility.Given the same concentration of the bulk gas at the initial stages, it is considered that when the temperature is higher, the adsorption equilibrium time will be shorter.The effectiveness of the Do technique 18 for surface diffusion parameter determination was tested with the sorption data for CO2 into the shale samples.Adsorption dynamics were measured at different bulk gas concentration steps but for the same range as explained in section 3.3.Plotting the parameter Ω from eq.21 versus the parameter Ξ, as shown in Figures 10(a-c) for all the samples gave rise to straight lines, which were then used to estimate the slopes and the intercepts.From eq.21 the intercept and the slope correspond to εΜDp and (1-εΜ)Ds respectively.Figures 10(a-c) show that the data gave rise to good fits to the various expressions for the Do technique 18 .
It should be noted that, if there was no surface diffusion occurring in the system, then a linear plot of the eq.24 must then have zero slope, with the intercept being the pore diffusivity.The carbon dioxide sorption data on the three Marcellus shale samples in the laboratory have indicated a good fitting to the theoretical uptake model.Therefore, this finding supports the implementation of the Do's theory 18 in heterogeneous systems like shales in which pore diffusion and surface diffusion may be determined within the system.Figure 11 shows Arrhenius plots for the surface diffusivities of CO2 on the three shale samples.The values of the surface diffusivities obtained in this work are similar in size to those obtained in previous studies in the literature, such as those found, by Karacan and Mitchell 45 , for CO2 in coal.The values of pore diffusivities obtained here for CO2 in the shales are ~10 -6 cm 2 /s.It should be noted that these results are consistent with previous findings 45 where pore diffusivity was not several orders of magnitude higher than surface diffusivity but of similar order of magnitude or just of the order of ten times larger.
The reason is probably due to the existence of small pores restricting the entrance of CO2 molecules into the shale.The resultant fitted Arrhenius parameters for the surface diffusivity at a monolayer coverage are shown in Table 5.
The characteristic isosteric heat of adsorption, ΔH, for each shale listed in Table 5 was calculated from the iQ-calorimeter.From Table 2 and 5, it can be seen that the heat of adsorption and surface area had a positive correlation with illite and the TOC content, while their correlation with quartz content was negative.Figure 12 shows the correlation of the characteristic heat of adsorption of CO2 with the group (Rn/r) d for the various depths of Marcellus Shale.The coefficient of determination for the fit to the data shown in Figure 12 was 0.999, and thus, a good fit between the fractal parameter (Rn/r) d and the heat of adsorption was obtained.The good quality of fit to a straight line shows that the data are consistent with the theoretical prediction given in eq. 5. Figure 13 shows a plot of the natural logarithm of the pre-exponential factor for the surface diffusivity at a statistical monolayer coverage against the fractal dimension of the surface of the shale from the three depths of the Marcellus field described above.Figure 13 also shows a fit of the data to eq.4.It can be seen that the quality of fit is high (r 2 =0.994) and thus the data are consistent with the fractal theory described in Section 2.1.From Figure 20 it can be seen that a good fit of the experimental data was obtained to eq. 6.It is noted that both the data for activation energy shown in Figure 13 and that for heat of adsorption in Figure 11 show a similar form of behaviour when plotted.
Various past authors 46 have suggested that there is a linear relationship between these two enthalpy parameters, and the above findings are consistent with this proposal.This means that, if heat of adsorption has a linear relationship with (Rn/r) d and activation energy of surface diffusion, it then follows that the surface diffusion activation energy also has a linear relationship with (Rn/r) d .The results shown in Figure 15 suggest the occurrence of the theoretically predicted 24 compensation effect for surface diffusion for CO2.Moreover, it was determined that the activation energy for surface diffusivity for CO2 produced a good fit for all three depths of Marcellus shale in accordance with the theoretical prediction, and hence the model is probably correct.In summary, the results indicate that the data obtained for CO2 surface diffusion on the Marcellus shale samples from different depths give outcomes consistent with the theoretical predictions.The fact that only one single fractal dimension for each shale was necessary to predict the surface diffusivity of CO2 for such heterogeneous samples, may be a combination of two factors.
First, critical path analysis 30 suggests that the observed rate of mass transport processes in networklike structures, such as shale rocks, is controlled by a particular set of critical conductances.The critical conductance is the lowest value in the network of pathways through which the mass transport flux actually migrates.In a shale rock, these critical conductances would correspond to particular patches of the internal pore-space surface that had the critical value of surface diffusivity.These critical surface patches would be the regions of the rock through which the surface diffusion flux is necessarily funnelled.This is because conductances above the critical value have the most rapid mass transport and, therefore, they are not rate-limiting, while conductances below the critical value would not contribute appreciably to the flux and they are effectively bypassed.Hence, the observed surface diffusivity is that characteristic of the set of critical surface patches with intermediate surface diffusivity.
Second, previous work suggests that nitrogen is much more of a specific adsorbate than is commonly suspected 47 .Even for supposedly relatively homogeneous materials, such as industrial sol-gel silicas, nitrogen has some tendency towards preferential adsorption, rather than being completely pervasively distributed across the whole surface 47 .Therefore, pore structural parameters obtained from nitrogen sorption can be more heavily-weighted towards certain regions of the void space.
Hence, if the critical patches for surface diffusion are part of the pore space surface predominantly characterised by nitrogen sorption, then even for a relatively heterogeneous material the structural parameters obtained from gas sorption will be predictive of the observed surface diffusivity.These two pore space regions are likely to overlap because the above findings suggest regions with intermediate heat of adsorption will correspond to regions of intermediate surface diffusion activation energy.Further, it is frequently found for shales that most of the accessible void space is predominantly associated with one component, namely the carbonaceous pores, with only a relatively small fraction of accessible porosity associated with illite or quartz phases 48 .Further, previous work 30 has indicated that for surface coverages around a statistical monolayer, due to surface heterogeneities and/or intermolecular interactions, adsorption on the surface is patchwise heterogeneous and surface coverages of the critical patches controlling mass transport approaches unity, irrespective of overall average surface coverage (at least for moderate partial coverages and above).
The combined diffusivity was calculated using the pore diameters attained through gas adsorption.The Bosanquet equation was applied in this calculation.The particle tortuosity was then calculated.Different values of tortuosity, ranging between 1.74 and 2.69, were obtained which are reasonable values for shale tortuosity given past results in the literature 49 .

Conclusions
It has been found that Do's theory for combined pore and surface diffusion gives rise to good fits to data even for highly heterogeneous adsorbents such as shales, and, thence, the surface diffusivity can be effectively estimated.The technique of parameter determination was demonstrated using sorption data for CO2 onto Marcellus shale.It has also been found that the experimental data shows that the surface geometry of the adsorbent determines the activation energy for CO2 surface diffusion and the heat of adsorption.The activation energy and pre-exponential factor were both found to be dependent on the surface fractal dimension, and this led to an expected compensation effect.It has been found that the surface diffusion data for CO2 on Marcellus shale from a variety of depths is consistent with the fractal theory developed by Rigby.The results confirm that the degree of structural heterogeneity of a shale surface determines the value of the Arrhenius parameters for surface diffusivity.Hence, the fractal theory for the structure-transport relation for surface diffusion can be used even for highly heterogeneous natural materials like shales.

Figure 1
Figure 1 Variation of heat of Adsorption of CO2 with coverage for Marcellus Shale 7804-7807ft

Figure 2
Figure 2 Schematic diagram of the Xemis gravimetric analyser

Figure 3 Figure 4 Figure 5
Figure 3 Adsorption isotherms for N2 on Marcellus shale

Figure 6
Figure 6 Fit of Fractal FHH eq.8(dashed line) to the adsorption isotherm data(symbols) for nitrogen on Marcellus shale 7804-7807ft

Figure 7
Figure 7 Isotherms for CO2 adsorption onto Marcellus Shale from a) 7804-7807ft b) 7834-7837ft c) 7864-7867ft measured at 10 0 C ,20 0 C, 30 0 C using Xemis apparatus.The lines shown are fits to the Langmuir isotherm model using parameters given in Table 4

Figure 8
Figure 8 Plot of fractional uptake of CO2 onto Marcellus shale 7804-7807ft.at three different bulk concentrations at 20 o C

Figure 11
Figure 11 Arrhenius plot of the surface diffusivities for CO2 on Marcellus shale from various depths

Figure 12
Figure 12 Variation of heat of adsorption of CO2 versus (Rn/r) d for various depths of Marcellus shale.The solid line is a fit of the data to eq. 6.

Figure 13 A
Figure 13 A plot of natural logarithm of the pre-exponential factor for the surace diffusivities for CO2 against the fractal dimension on a variety of depths of Marcellus shale Figure 14 shows a plot of the activation energy for the surface diffusivity at statistical monolayer coverage against the group (Rn/r) d for the shale sample from three depths of the Marcellus shale core.From Figure20it can be seen that a good fit of the experimental data was obtained to eq. 6.It is noted that both the data for activation energy shown in Figure13and that for heat of adsorption in Figure11show a similar form of behaviour when plotted.

Figure 14 A
Figure 14 A plot of activation energies for the surface diffusivities of CO2 against (Rn/r) d on different depths of Marcellus shale Figure 15 depicts a plot of the natural logarithm of the pre-exponential factor for the surface diffusivity at statistical monolayer coverage against the corresponding activation energy for surface diffusion on the shale samples.Figure 15 also shows a fit of the experimental data to a straight line of the form of eq. 7.

Figure 15 A
Figure 15 A plot of compensation effect of the Arrhenius parameters for the surface diffusivities of CO2 on Marcellus shale

Table 1
Parameters values for shape particle

Table 2
Results of Marcellus Shale characterisation

Table 3
Parameters obtained from fractal FHH analyses of the N2 gas adsorption isotherms

Table 4
Langmuir isotherm parameters for CO2 on Marcellus shale derived from isotherms measured at the indicated different temperatures

Table 5
Arrhenius parameters for the surface diffusivity of CO2 on Marcellus shale samples