Proposed Model for Shale Compaction Kinetics: Forward Depositional Model

1. Introduction

A kinetic equation previously presented [1] allowed determining the time history of porosity of a given shale section in terms of four constant parameters, m,n,A and activation energy E :

$${\left.\frac{\partial \left(ln\right(1-\phi \left)\right)}{\partial t}\right|}_{\sigma }={\left(\phi \right)}^{4m/3}{(1-\phi )}^{4n/3}A{e}^{-E/\left(RT\right)}.$$

(1)

The right hand side is analogous to customary chemical kinetic equations, and anticipates that compaction occurs because some shale component either dissolves into pore water, or moves along the mineral-water interface, and then precipates to infill pore space. Constant ’A’ has units /time, and is a lumped constant which includes collision frequency, changing grain-to-grain contact areas, varying morphologies, geometical restrictions, and complications among inhogeneous solids, many of which may not participate.

The first two dimensionless parameters, m and n, were estimated by first noting that the time derivitive on the left would be zero at zero porosity, so that the first non-vanishing term in a Taylor series representation would be proportional to $\phi$. For a particular actual shale porosity profile, m and n can be varied to make $\phi *\left[{\left(\phi \right)}^{-4m/3}{(1-\phi )}^{-4n/3}\right]$ close to constant over the deeper part, as was illustrated for six sections. The exponents are not near 2/3, as expected for mathematically flat interfaces, but reflect the rough interfaces at the molecular level. The ’deeper part’ restriction arises because the compaction mechanism is pressure solution, and this necessarily disappears at the surface where there is no overburden. Other mechanisms like mechanical compaction and bioturbation contribute more at shallow depths. The six example shales used here are the same as in [1] and the same m, n and temperature profiles are used. The included temperatures are extended to 25 deg C [2] instead of the less inclusive 30-40 deg C as previously. This changes A and E slightly. The A values reported here for the largest deposition rate replace those in [1]. ’Deeper part’, as before also excludes porosities greater than 0.5 if measured or reported.

The activation energies to be derived here show some variation. This variation is, however, small compared to the range among common sedimentary minerals[3]. Current understanding and information is not complete enough to choose among minerals and conditions. For instance, there are six varieties of quartz for which different solubilities have been measured at 25 dec C, and it may be that derived pressure solution parameters will also vary for this one mineral. Mineral activation energies E for solubility also vary with pH [3]. Other mineral cements such as CaCO3 are potential candidates. Variations observed here could be a consequence of the model, a shale property, or experimental uncertainty. Deposition rates, discussed here, are another source of variation, as are temperature variations over geologic time. Temperatures may vary non-linearly with depth presently or in the past.

Probably the source of pressure solution/compaction is an abundant mineral fraction. $Si{O}_2$ qualifies, as it is usually 20-50% of shales, the other largest fraction being clay minerals which also contain $Si{O}_4^{-4}$ as framework. Reactions such as montmorillonite to illite can sometimes be a source of $Si{O}_2$. Additionally, quartz is the last major mineral to freeze out of igneous materials, and thus is briefly a pore fluid and coats minerals which solidified earlier in the Bowen’s reaction series.

As a footnote, a compendium of shale porosities [4] shows only greater than zero porosities for shales or slates.

If the actual deposition history of a shale section were known, it could be put into the left hand side (LHS) of equation 1 and the right hand side (RHS) parameters A and E could be determined. As mentioned, m and n will be unchanged here.

2. Development

As the deposition histories are not known for our examples, the missing link between deposition time and depth is required to continue. To proceed, a constant rate of solids deposition, r, in mm/y or km/Ma, is assumed for the entire section, allowing A and E to be obtained. This is done successfully for a large range, 0.2 to 20 km/Ma, for each of the six example shales, which insures that actual complicated depositional histories can be covered. The time-depth link is:

$$dt=(1/r)(1-\phi (z\left)\right)dz.$$

(2)

Intergrating from the section top, $z_0$ to the bottom $z_b$ gives the total time for deposition,

$$t_{z=b}-t_{z=0}=(1/r){\int }_{z=z_0}^{z=z_b}(1-{\phi }_z)dz.$$

(3)

Then the geologic age at any depth $z_i$ is

$$t_{z=z_i}=[t_{g1}+t_{g0}]/2-[t_{z=b}-t_{z=0}]/2+(1/r){\int }_{z=z_0}^{z=z_i}(1-{\phi }_z)dz,$$

(4)

Here $[t_{g1}+t_{g0}]/2$ is the middle of the known geologic time interval $t_{g0}$ through $t_{g1}$ in which deposition occurred. All times are reckoned positive backward from the present.

The integral (4) can be numerically evaluated by dividing $z=0$ through $z=b$ into many, N, equal small intervals ${\delta }_z=z_i-z_{i-1}$,

$${\int }_{z=0}^{z=z_i}(1-{\phi }_z)dz=\left({\delta }_z\right)\sum _{i=1}^{N_{z_i}}(1-\phi \left(z_i\right)).$$

(5)

Substituting Equation (5) into (4) gives the time $t_{z=z_i}$ since deposition, as needed in Equation (1). For each i = 1 through N,

$$[ln(1-{\phi }_i)-ln(1-{\phi }_0)]/t_i={\phi }_i^{4m/3}{(1-{\phi }_i)}^{4n/3}Aexp[-E/\left(R{T}_i\right)]$$

(6)

where subscripts containing ’z’ or ’$z_i$’ are shortened to ’i’.

This is the desired result, allowing parameters A and E to be obtained from the depositional history of the shale section on the LHS.

3. Results

For these calculations, for this mechanism, the porosity extrapolated to z=0is estimated for simplicity as 0.5 for all wells/sections [5,6,7,8,9]. Reported actual surface porosities are sometimes more or less, or were omitted.

Moving all porosity terms to the left of Equation (6) and taking the logarithm of each side yields a linear equation with intercept ln(A) and slope E. Table 1 gives the derived E and A for a broad range of realistic deposition rates, r, for each section or well studied. For each section the range of E is small compared to laboratory measurements of quartz[2], 24+-15kJ/mol, and all are in this range.

As contrasting examples which have abundant data, the Arrhenius plot is given for the Maracaibo well[9], Figure 1, and the Macran1 section [6], Figure 2. A low solids deposition rate of 0.5 km per million years is used. Figures for all six sections were previously given [1] for very high, >20 km/Ma, deposition rates such that $\Delta$ is maximum.

The Maracaibo porosities were measured on hand samples. Porosities were fitted with a straight line. The line omitted large near-surface porosities and made deeper porosities seem to go to zero too quickly, as was recognized[9]. Experimental porosities above about 36% were left out. The lake pH is slightly less than 7.0. The well was drilled in or near shallow lake Maracaibo in a tropical climate. In consrast,the Macran-1 data were derived from seismic data with the maximum porosity reported as 69%. The environment was deep cold ocean with pH 8.1.

This completes the forward model of deposition and compaction of these six shales. In conclusion, pressure solution of quartz can account for vertical compaction of shales quantitatively. Other examples may indicate carbonate or iron cements. The procedure is expected to work similarly for most shale sections, whether overpressured or not, if not strongly disturbed.

Table 1. Section deposition rates and parameter fits.

Example	r	E	A	${\mathit{\phi }}_{\mathit{e}\mathit{r}\mathit{r}}$
Properties	km/Ma	kJ/mol	$/s*{10}^{-12}$	std dev
Akita	0.2	30.8	43	0.04
5-149 Ma	0.3	31.1	49	0.04
m,n = 0.95,1.	1	31.6	59	0.04
t >= 25$°$C	4	31.8	62	0.04
${\phi }_{av}$=0.269	20	31.8	64	0.04
Makran1	0.2	17.2	1.4	0.01
2.6 - 66 Ma	0.3	18.1	2.1	0.01
m,n = 0.85,1.0	1	19.7	4.0	0.01
t >= 25$°$C	4	20.3	5.1	0.01
${\phi }_{av}$=0.168	20	20.4	5.4	0.01
Makran2	0.2	14.4	0.9	0.06
2.6 - 66 Ma	0.3	15.4	1.3	0.06
m,n = 0.85,0.95	1	17.8	2.5	0.06
t >= 25$°$C	4	17.5	3.2	0.06
${\phi }_{av}$=0.126	20	17.7	3.4	0.06
"SuluSea"	0.2	12.3	0.55	0.01
0 - 23 Ma	0.3	13.9	1.1	0.01
m,n = 0.9,0.8	1	17.4	4.6	0.01
t >= 25$°$C	4	19.0	8.9	0.01
${\phi }_{av}$=0.163	20	19.5	11.0	0.01
Oklahoma	0.2	31.4	44	0.04
254-323 Ma	0.3	31.5	46	0.04
m,n = 0.85,1.0	1	31.6	48	0.04
t >= 25$°$C	4	31.6	49	0.05
${\phi }_{av}$=0.130	20	31.6	49	0.05
Maracaibo	0.2	18.4	2.4	0.03
2.6-66 Ma	0.3	19.1	3.1	0.03
m,n = 0.9,0.9	1	20.0	4.5	0.03
t >= 25$°$C	4	20.3	5.2	0.03
${\phi }_{av}$= 0.212	20	20.4	5.3	0.03

Table 1. Section deposition rates and parameter fits.

Example	r	E	A	${\mathit{\phi }}_{\mathit{e}\mathit{r}\mathit{r}}$
Properties	km/Ma	kJ/mol	$/s*{10}^{-12}$	std dev
Akita	0.2	30.8	43	0.04
5-149 Ma	0.3	31.1	49	0.04
m,n = 0.95,1.	1	31.6	59	0.04
t >= 25$°$C	4	31.8	62	0.04
${\phi }_{av}$=0.269	20	31.8	64	0.04
Makran1	0.2	17.2	1.4	0.01
2.6 - 66 Ma	0.3	18.1	2.1	0.01
m,n = 0.85,1.0	1	19.7	4.0	0.01
t >= 25$°$C	4	20.3	5.1	0.01
${\phi }_{av}$=0.168	20	20.4	5.4	0.01
Makran2	0.2	14.4	0.9	0.06
2.6 - 66 Ma	0.3	15.4	1.3	0.06
m,n = 0.85,0.95	1	17.8	2.5	0.06
t >= 25$°$C	4	17.5	3.2	0.06
${\phi }_{av}$=0.126	20	17.7	3.4	0.06
"SuluSea"	0.2	12.3	0.55	0.01
0 - 23 Ma	0.3	13.9	1.1	0.01
m,n = 0.9,0.8	1	17.4	4.6	0.01
t >= 25$°$C	4	19.0	8.9	0.01
${\phi }_{av}$=0.163	20	19.5	11.0	0.01
Oklahoma	0.2	31.4	44	0.04
254-323 Ma	0.3	31.5	46	0.04
m,n = 0.85,1.0	1	31.6	48	0.04
t >= 25$°$C	4	31.6	49	0.05
${\phi }_{av}$=0.130	20	31.6	49	0.05
Maracaibo	0.2	18.4	2.4	0.03
2.6-66 Ma	0.3	19.1	3.1	0.03
m,n = 0.9,0.9	1	20.0	4.5	0.03
t >= 25$°$C	4	20.3	5.2	0.03
${\phi }_{av}$= 0.212	20	20.4	5.3	0.03

A check can be made to see how well the derived parameters m,n,A and E reproduce the known shale porosity profiles. Equation (6) rearranges to

$${\phi }_{i_{est}}=1-(1-{\phi }_0)exp\left[t_i{\phi }_i^{4m/3}{(1-{\phi }_i)}^{4n/3}\right)Aexp[-E/\left(R{T}_i\right)]$$

(7)

with all four fixed parameters on the right, along with the fitted porosities ${\phi }_i$, to give the estimated porosity at this depth on the left. The mean of the errors is zero percent for each well. The standard deviations are given in Table1. Values of A in Table 1 supercede values in [1], which were wrong but did not affect conclusions.

The Macron 2 section, in the Macron accretionary prism, has reduced porosities which are $60-70\%$[6] of the Macron 1 porosities in the nearby abyssal plain with similar sediment. This reduced porosity is thought due to greater horizontal stress. Usually the excess horizontal stress is not known. A quasi-equivalent vertical stress can be obtained as a multiplier, w, of the vertical stress, $\sigma$. The proposal is that if $\sigma$ produces the Macron 1 porosities,${\phi }_{M1}$, then $w\sigma$ produces the Macron2 porosities, ${\phi }_{M2}$ :

$$\Lambda =\sigma /(1-{\phi }_{M1})=w\sigma /(1-{\phi }_{M2})$$

(8)

which gives

$$w=(1-{\phi }_{M2})/(1-{\phi }_{M1}).$$

(9)

Using average porosities over the respective depth intervals 1 to 3.5 km gives w = 1.06 using the porosity curve fits, and w = 1.068 using the seismically inferred porosities. The infered quasi vertical stress is thus about 6.4% greater than $\sigma$ to include the excess horizontal stress with the vertical. This procedure qualitatively bridges a gap between the effects of vertical and horizontal stresses, ignoring anisotropy and timing factors.

Similar horizontal stress effects were thought to account for reduced porosities in the Oklahoma well [8,9]. There is no ’normal’ well for conparison in this case.

The Maracaibo location was thought to be free of excess horizontal stress[9]. This lake is fed by tropical rivers, with lake water having pH 7 or slightly lower, as compared to ocean waters of pH 8.1. Solubilities and rates are sensitive to pH [2,3], Experiments to determine whether reactions are endo- or exothermic could be made.

A three dimensional model is required to fully understand tectonically active areas.

4. Discussion

The present model brings up unanswered questions, which are listed as potential experimental problems.

1- Is more than one type of mineral solution occuring among different wells or in each well.

2- Variations of stress, both positive and negative, lead to increased solubility of shale minerals[2]. Are local vertical and horizonal stresses the only drivers. What contributions do quakes and crustal movements make. A linked question is, what is the minimum time that a change in stress can occur and produce porosity reduction.

3. How small can effective stresses be, and how small can a time increment be to produce pressure solution in shales.

The successful forward modeling here of a wide range of deposition rates for each shale section indicates that depositional rate for shales is of secondary importance. Temperature, stress, and stress variations are more important.

These questions highlight that tidal forces are possibly part of the proposed pressure solution mechanism and were initially ignored here. To state the problem again, the vertical force per unit area of supporting grains, $\Lambda$, is given by the ratio of the vertical force per unit area, $\sigma$, divided by the area of the supporting grains, $(1-\phi )$:

$$\Lambda \equiv (S-p)/(1-\phi )=\sigma /(1-\phi ).$$

(10)

Here the vertical force, $\sigma$, is the difference between the overburden, S, and the pore fluid pressure, p. Tidal forces are part of $\sigma$ and should be separated as the time derivative is intrinsically never zero, while the part of $\sigma$ due to Earth’s gravitational field alone, $g_E$ is essentially independent of time. Separating the tidal forces as $g/g_E$,

$$\sigma =[\sigma *g_E/g](g/g_E)$$

(11)

or

$$\sigma =\left[\overline{\sigma }\right](g/g_E).$$

(12)

Here $\overline{\sigma }$ is the vertical stress toward the center of the earth, without Moon and Earth tidal pulls, and $g/g_E$ is the scalar vertical component of these tidal forces which vary continuously. The latter should realistically be replaced by the vector pulls of Moon and Sun on every part of a section, these vectors rotating 360 degrees every lunar or solar day. Taking the time partial derivative with this in mind gives

$${\left.\frac{\partial \left(ln\Lambda \right)}{\partial t}\right|}_{\overline{\sigma }}=-{\left.\frac{\partial \left(ln\right(1-\phi \left)\right)}{\partial t}\right|}_{\overline{\sigma }}+"d(g/g_E)/dt"$$

(13)

The quote marks indicate this scalar term stands in for the three-dimensional vector derivative. If the shales are primarily elastic over small macro-sized distances but have imperfactions at sites where a few molecular diameters of extension collect, these sites might be a place to consider for pressure solution. Elaborating further, the ’equator’ of the Moon - Earth orbit can have a large vertical force component near the Moon-Earth axis, but as seen from the Moon, horizontal stretching forces toward the Earth’s edges predominate, especially toward the ’poles’. This difference might be experimentally observable. Mutatis mutandis for the Sun. Additionally, remembering that the Oklahoma well is Pennsylvanian age, the Moon was substantially closer to the Earth several hundred million years ago - with Oklahoma near the equator - and gravitational stretching forces vary inversly as the cube of the distance.