Proposed Model for Shale Compaction Kinetics

: Shales are the most abundant class of sedimentary rocks, distinguished by being very ﬁne-grained, clayey, and compressible. Their physical and chemical properties are important in widely different enterprises such as civil engineering, ceramics, and petroleum exploration. One characteristic, which is studied here, is a systematic reduction of porosity with depth of burial. This is due increases in grain-to-grain stress and temperature. Vertical stress 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 fractional porosities less than 0.5. Single activation energies for each well are obtained in the range 15–33 kJ/mole, close to the approximate pressure solution of quartz, 24 kJ/mol. Values of m and n are in the range 1 to 0.8, indicating nearly fractal water-wet pore-to-matrix interfaces at pressure solution locations. Results are independent of over- or under-pressure of pore water. This model attempts to explain shale compaction quantitatively. For the petoleum industry, given porosity-depth data for uneroded sections and accurate activation energy, E, paleo-geothermal-gradient can be inferred and from that organic maturity, indicating better drilling prospects.


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 previously published porosity data for shale sections a chemical-type expression describing the time evolution of the vertical grain-to-grain pressure, or frame pressure.

Proposed Kinetic Equation
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.314 J/(mol·deg K), T is temperature in Kelvin. The fractional change in Λ thus proposed and supported is 1.33n Ae −E/(RT) (4) or equivalently The Arrhenius frequency factor, A, is not expected to depend strongly 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-to-grain 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 more inclusive approaches see [4]). Fractional water saturation, S w , might reduce compaction but is set to 1 and omitted hereafter.

Data
Porosity versus depth profiles from the literature [5][6][7][8][9] for six shale sections are used here. The shales are unrelated and not chosen because of similarities. Sample-based profiles were included in a recent comprehensive survey of 21 geologic sections and laboratory studies of clays [10]. All profiles are included in an engineering correlation [11] which included two seismicly-inferred porosity profiles [6]. They were chosen because geologic age ranges, basin types, and basin names were available [11]. All porosity-depth profiles will be referred to as 'wells' for brevity. Map coordinates were known only for the two Makran seismically-inferred porosity-depth profiles. For actual wells, no formal well names were given, and data from several local wells were combined. Temperatures are not available for any wells. Paleo surface temperature estimates are made on the basis of paleoenvironments, paleolatitude, and geologic age, Table 1. Sediment surface temperatures of 3 • C are assigned to wells from abyssal plains. Shallow sea [8] or lake environments [9] are assigned expected land paleotemperatures, corrected for paleolatitudes, and for any intervening high temperature events such as the Eocene high [12]. If two estimates are made they are averaged. The sediment surface temperature used for the Sulu Sea, 13 • C, was known for the deepest point [13]. The paleo-surface-temperature estimates are less important to this study than are paleo-geothermal-gradients, as the compaction mechanism discussed hereafter becomes noticable deeper at the higher well temperatures and at 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 0.03 • C/m were initially assigned. The well surface and temperature gradient estimates used here are well within present day experience [14]. A paleo-geothermal-gradient of 0.05 • C/m was assigned to the Akita well, and to the Maracaibo back-arc well, on the basis of necessary organic maturity [5]. In Northeast Oklahoma nearby sedimentary rock mineralization and granite outcrop [15] suggest a higher paleo-geothermal-gradient, perhaps 0.06 • C/m as assigned.   It is perspicuous to view porosities versus Λ, Figures 7-12, which is the force pressing the grains together, and is a more fundamental variable than depth, common to all wells. The Λ axis is approximate. (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 upper surface. The grain density and water density are ρ g and ρ w with (ρ g − ρ w ) = 1.67g/cc. 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. Parameters m and n are fixed so that the well porosity profile in the deeper portion of a well matches the surface reactant terms on the right of Equation (5), where c is a constant. In chemical terms, this means that the ratio of the reacting surfaces terms to pore volume is constant, and indicates that water is a necessary reactant. These figures indicate that the proposed reaction, Equation (4), is not rate-limiting at temperatures below 30-40 • C, and that longer geologic times (Table 1) or lateral compression can produce greater compaction effects. Table 2 gives m and n, and shows that the scatter of the porosities is reduced using the ratio of the entropy of the porosities after the 'porosity correction' to the entropy before the 'correction'. The convergence of the transformed porosities is good above 30-40 • C and at porosities below 0.5. 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. The example wells are scattered around the World. The main task remaining is to find whether Equation (5) predicts the existance of rate-limiting reaction(s) charactorized by one or more activation energies, E, for these wells. Moving the porosity terms to the left hand side (LHS) of Equation (5) gives The left side needs to be evaluated for the well porosity profile in order to obtain E and A. The time derivative-or of equal practical use the derivative of porosity with respect to time-can potentially be obtained from laboratory experiments. Alternatively the time derivative can be obtained if the sediment geologic ages are known at close depth intervals along with the known porosity profile. However only very broad bounding geologic ages are known for these six wells, e.g., the Oklahoma sediments were deposited from early Pennsylvanian to late Permian times which start 323 Myr ago and end 252 Myr ago, Table 1. As these approaches are not currently available, an approximation to the time derivative is used.
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 The subscript '0' indicates near the top of the well but at temperatures above ≈40 • C and porosities less than 0.5. (A 60 • C cut-off produces similar results.) Because only a maximum range is known for ∆t, the average ∆t for the entire section is used in the initial estimates of E in Table 1. Thus the section is supposed to be deposited in a time interval short compared to ∆t with ln(∆t/∆t 0) zero in this approximation although it is positive for sediments as normally deposited. Plotting the log of the left hand side versus −1/RT, gives E and the log of A for the six 'wells', Figures 19-24

Results
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 [16] give 24 ± 15 kJ/mol, suggesting that this is the rate-limiting reaction occurring in shales.

Discussion
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. The higher activation energies obtained for Akita and Oklahoma wells are well within the 44 ± 4 kJ/mol activation energy upper bound for solubility of quartz with no load [6], and the Oklahoma well is known to have lost at least 360 m of overburden [8].
With this E av an Oklahoma paleo-geothermal-gradient of 0.12 • C/m was inferred although this is too high if unloading by erosion has been the reason for the higher E. More example wells are needed to sort out this problem. Using an average E 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 Macran 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. The Makran data, Figures 2 and 3, show that lateral forces in the accretionary wedge (Makran2) produced porosities reduced from nearby abyssal plain (Makran1) 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 also given by Equation (4).

Conclusions
Equation (4) explains shale compaction at temperatures above ≈35 • C with activation energies E ≈ 23 ± 6 kJ/mol, close to 24 ± 15 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 0.8 and 1.0 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. For obtaining activation energy it is assumed that deposition occurs over a comparatively short time period compared to the geologic age. Porosities greater than 0.5 are not used in computing activation energy as near-surface sediments have chaotic porosities. Provisional paleotemperature gradients are used to obtain activation energies which are averaged to bootstrap to hopefully better paleotemperature gradients. Given porosity data for an uneroded section and an average activation energy, E av , one can infer or rank maximum paleo-geothermal-gradients and from that organic maturity, thus avoiding unnecessary drilling.