Shale Compaction Kinetics Authors

The grain-to-grain stress vertically in sediments is given by the overburden less the pore ﬂuid pressure, σ , divided by the fraction of the horizontal area which is the supporting matrix , (1 − ϕ ), where ϕ is the porosity. It is proposed that the fractional reduction of this ratio, Λ, with time is given by the product of ϕ 4 m/ 3 , (1 − ϕ ) 4 n/ 3 , and one or more Arrhenius functions A exp( − E/RT ) with m and n close to 1. This proposal is tested for shale sections in six wells from around the world for which porosity-depth data are available. Good agreement is obtained above 30-40 C and porosities less than 0.5. Single activation energies for each well are obtained in the range 15-33 kJ/mole, close to pressure solution of quartz, 24 kJ/mol. Values of m and n are in the range 1 to 0.8, indicating nearly fractal pore-matrix spaces and water-wet interfaces. Results are independent of over- or under-pressure of pore water. This model explains shale compaction quantitatively. Given porosity-depth data and accurate activation energy, E, one can infer paleo-geothermal-gradient and from that organic maturity, thus avoiding unnecessary drilling.


Introduction
For overburden S and pore pressure p, the vertical stress difference can be changed by sedimentation, erosion or fluid flow through rock in this one-dimensional model. Similarly σ divided by the fractional area of solid matrix (1 − ϕ) supporting this stress can be relieved by pressure solution, breakage, relative movement of grains or fluid flow.
Differentiating with respect to time gives The first term deals with local porosity changes, while the second term deals with deposition, erosion, and fluid flow. The second term was treated previously [1] [2] [3]. Focus here is on the first term, proposing and supporting with data a chemical-type expression describing the time evolution of the vertical grain-to-grain pressure, or frame pressure.

Kinetic Equation and Experimental Tests
It is proposed that the reactions to reduce the frame pressure are proportional to the frame pressure Λ, to the pore surface area, ϕ 1.33 , to the frame surface area, (1 − ϕ) 1.33 , and to one (or a sum of) Arrhenius factors, A exp −E/RT . The exponent 1.33 = 4/3 is the fractal dimension of a percolation front and twice the surface dimension 2/3 of smooth solids. R is the ideal gas constant (8.314J/(mol degK), T is temperature in Kelvin. The fractional change in Λ thus proposed and supported is ∂(lnΛ) ∂t σ = −(ϕS w ) 1.33m (1 − ϕ) 1.33n Ae −E/(RT ) (4) or equivalently The Arrhenius frequency factor, A, is not expected to depend on T for the surface reaction. Here m and n are expected to be close to 1.0 and greater than 0.5, contact areas being less than fractal because of the finite size of atoms. m and n may be different because water-to-grain and grain-tograin interfaces differ. Notably this mechanism is primarily a dimensional argument and lacks terms relating to the sizes, shapes, and composition of pores and grains. ( For different approaches see [4] ). Fractional water saturation, S w might reduce compaction but is set to 1 and omitted hereafter. Equation 5 indicates that if both sides are multiplied by the 'surface correction', the right hand side (RHS) will be independent of porosity, ϕ av is the average over the 'well' which helps visualization by minimizing movement of data on individual plots and among plots.  [9]. The examples are referred to here as 'wells' although they often combine porosities from several local wells. The Macran porosity data sets are inferred entirely from seismic data, not well data. The "Sulu Sea" data are a composite of several mostly unspecified wells. the porosity data versus an approximate Λ horizontal axis. ( The approximation p w = gρ w z for normal water pressure is used to approximate the Λ axis as Here z is sample depth and depth 0 is the sediment-water or -air interface, briefly referred to as the mudline. The grain density and water density are ρ g and ρ w . The water pressures in these wells were not given.) Any mismatch of experimental porosity with approximate Λ is not important for illustrative purposes as porosity data points in non-normally-pressured zones are only shifted horizontally, e. g., overpresssured zones appear shifted to the right, and transformed porosity points are shown here next to be nearly constant. These figures indicate that the reaction is slow at temperatures below 30-40 C., and that longer geologic times (Table 1) produce greater compaction effects. At porosities greater than 0.5 Equation 4 can have multiple porosity values give the same number, as it has the logistic form. This agrees with deep sea studies which show porosities are chaotic [6] at porosities greater then 0.5. Porosities greater than 0.5 are excluded here in estimating activation energies.

Results
Equation 4 is supported in that, after the 'correction', porosities in each well have approached a common value except for shallow chaotic low-temperature portions. Table 1 gives the ratio of the scatter (entropy) of the porosities after the correction to the entropy before the corrections. The convergence of the transformed porosities is good above 30-40 C and at porosities below 0.5. An estimate of activation energy, E, can be made by plotting the logarithm of the left hand side of Equation 9 versus −1/RT . The common ϕ av terms in Equation 7 are cancelled here The partial derivative requires an approximation. Breaking geologic time τ into many, N, small time increments δt it is seen that each term in the sum below is the instantaneous time derivative ln((1 − ϕ i )/(1 − ϕ i−1 /δt for that time. The sum is the average over time of these instantaneous time derivatives and, like the pages of a book, add up to give the (large) finite-difference approximation on the right below, equivalent to only the starting and ending pages of the book as it were.
Here we have put ϕ N = ϕ, the current porosity. ϕ 0 would be the seabed porosity to prevent negative values of the logarithm, but is 0.5 to analyze deeper shales and avoid chaotic behavior close to the seabed.
Porosity changes due to lateral forces and vertical forces are obviously included in the summation of the average, Equation 12, which is used in place of the partial derivative, This partial derivative on the left of Equation 11 would have come up in the same way if σ, Equation 3, represented horizontal forces rather than vertical forces. This approximation, or replacement, thus represents overall porosity reduction due to the combined forces and temperatures acting, in natural order, and including changes in overburden and pore fluid pressure, to reduce the potential energy Λ. This replacement is more appropriate than the partial derivative for computing activation energy, and further cannot be avoided since the actual porosity data result from all these factors. Putting this approximation into Equation 9 gives Multiplying the LHS with a constant K to compensate for a possible mis-estimate would not change the estimate of E, as taking the slope of the logarithm of Equation 12 drops out K. Assuming deposition occurs over a time short compared the geologic age range, Table 1, with no important unconformities, τ is approximately constant also. Thus where the * superscript means at temperatures above ≈ 40C and porosities less than 0.5. (A 60C cutoff produces similar results.) With this E a rough approximate A is obtained, Table 1, This A estimate would contain poorly known τ and any K. K = 1 is used. The average of the top and bottom geologic time boundaries for a 'well' is taken for τ in computing A. Plotting the log of the LHS versus −1/RT gives E for the six 'wells', Figures 19-24. Arrhenius plot with slope E kj/mol As shown in Table 1 the E's vary from 15 to 33 kJ/mol, averaging 23+-6 kJ/mol. Laboratory measurements of the pressure solution of quartz aggregates [12] give 24 kJ/mol, suggesting that this is the rate-limiting reaction occurring in shales.

Discussion and Data
Map coordinates were known only for the Makran 'wells'. Temperatures were not available for any wells. Paleo mudline temperature estimates are made on the basis of paleoenvironments [5], paleolatitude, and geologic age, Table 2. Mudline temperatures of 3 C are assigned to wells from abyssal plains. Shallow sea or lake environments are assigned expected land paleotemperatures, corrected for paleolatitudes, and for any intervening high temperature events such as the Eocene high [13]. If two estimates are made they are averaged. The mudline temperature used for the Sulu Sea, 13 C, was known for the deepest point [14]. These paleo-mudline-temperature estimates are less important to this study than are paleo-geothermalgradients, as the compaction reaction found here becomes noticable deeper at the higher 'well' temperatures and porosities less than 0.5. Because the Earth is cooling as radioactivity and initial energy of assembly decline, and because these wells were drilled/investigated in areas where temperatures were supposed to be high enough to generate petroleum, minimum paleo-geothermal-gradients of .03 C/m were initially assigned. The 'well' mudline and temperature gradient estimates used here are well within present day experience [15]. A paleo-geothermal-gradient of .05 C/m was assigned to the Akita 'well', and to Maracaibo back-arc well, on the basis of necessary thermal maturity [5]. In Northeast Oklahoma nearby sedimentary rock mineralization and granite outcrop suggest a higher paleo-geothermal-gradient, perhaps .06 C/m as assigned, Table 1. Averaging the E's computed from all these provisional temperature gradient assignments, along with the one laboratory measurement, gives an average activation energy of 23+-6 kJ/mol. With this E av an Oklanoma paleo-geothermalgradient of 0.12 C/m was inferred. Using this average assumes the same mechanism occurs in all wells. With this E av , paleo-geothermal-gradients for the other wells were also obtained. These 'improved' estimates with other recalculated parameters are shown in Table 3 as far as data convergence allowed. Differences in porosity determination method might have impacted activation energy determinations, e.g., Oklahoma and Maracaibo porosities were measured on hand samples, while as mentioned the Macron porosities are inferred from seismic data. Current temperature profiles are known privately for thousands of wells. This information would help gain an accurate E av , especially for deep-water wells where the paleo-mudline temperature is probably ∼ 3 C. Laboratory data for pressure solution of quartz and other minerals can be expanded to higher temperatures and pressures to check current ideas. The Makran data, Figures 2 and 3, show that lateral forces in the accretionary wedge (Makron2) produced porosities reduced from nearby abyssal plain (Makron1) porosities. There the major lithology need not be clay throughout, as good reflectors were needed to produce data. The low Oklahoma porosities were also thought to have been reduced by lateral forces [8] [9]. Directionally variant lateral forces are often observed in horizontal drilling [16]. The compaction response mechanism to both vertical and horizontal forces is probably the same, Equation 4.

Conclusions
Equation 4 explains shale compaction at temperatures above ≈ 35 C with activation energies E ≈ 23+-6 kJ/mol, close to 24 kJ/mol for the pressure solution of quartz. Problems of mineralogy, grain shapes and sizes, permeability, overpressures and pore connectivity are by-passed. The results support nearly fractal pore interfaces with m and n between .8 and 1. The reduction in Λ is also modeled as proportional to geologic age to the first power. By Equation 12 the one-dimensional problem is recognized as three spatial dimensions. In obtaining activation energy the geologic age approximately drops out if there are no unconformities. The problem of chaotic porosity data at porosities greater than 0.5 is recognized and avioded. Provisional paleotemperature gradients are used to obtain activation energies which are averaged to bootstrap to hopefully better paleotemperature gradients. Given porosity data from the known or, in cases of later erosion, an estimated pre-erosional depth, along with an average activation energy, E av , one can infer or rank maximum paleo-geothermal-gradients and from that organic maturity, thus avoiding unnecessary drilling.

Acknowledgments
Gary M. Hoover clarified obscurities and limitations in the progression of this paper.

Author contributions
All work was carried out by corresponding author James Edward Smith 85, with physicist son Edward Millard Smith-Rowland as backup due to age and Covid 19 threat.

Funding
This research received no external funding.

Conflicts of interest
The authors declare no conflict of interest.