Rock acoustics of diagenesis and cementation

We simulate the effects of diagenesis, cementation and compaction on the elastic properties of shales and sandstones with four different petro-elastical theories and a basin-evolution model, based on constant heating and sedimentation rates. We consider shales composed of clay minerals, mainly smectite and illite, depending on the burial depth, and the pore space is assumed to be saturated with water at hydrostatic conditions. Diagenesis in shale (smectite/illite transformation here) as a function of depth is described by a 5th-order kinetic equation, based on an Arrhenius reaction rate. On the other hand, quartz cementation in sandstones is based on a model that estimates the volume of precipitated quartz cement and the resulting porosity loss from the temperature history, using an equation relating the precipitation rate to temperature. Effective pressure effects (additional compaction) are accounted for by using Athy equation and the Hertz-Mindlin model. The petro-elastic models yield similar seismic velocities, despite the different level of complexity and physics approaches, with increasing density and seismic velocities as a function of depth. The methodology provides a simple procedure to obtain the velocity of shales and sandstones versus temperature and pressure due to the diagenesis-cementation-compaction process.


Introduction
Diagenesis is the chemical and physical process by which sediments which are buried in the Earth crust are compacted by a process where minerals precipitate from solution due to cementation, grain rearrangement, with a reduction in porosity (Pytte and Reynolds, 1989;Walderhaug, 1996;Lander and Walderhaug, 1999). Thus, sediments (clay and sand) become rocks by lithification, i.e., the precipitated mineral create bonds between grains, and sand becomes sandstone and clay-rich sediments become shale by this combined diagenesiscementation process. We do not consider the phenomenon by which metamorphic rocks are formed, which typically occurs after that process, beyond 200 • C, and before melting, roughly 800 • C, depending on the minerals . The diagenetic process requires the flow of water in geological time for the minerals to precipitate and crystallize. Typical cements include quartz, calcite and clay minerals.
The acoustics of diagenesis has scarcely been attacked. We can mention Draege et al. (2006), who presented a methodology for the estimation of the shale stiffnesses in the transition zone from mechanical compaction to chemical compaction/cementing. The model showed consistent results when comparing vertical P-and S-wave velocities with log data.
Another type of diagenesis process is the creation of kerogen and bitumen from which hydrocarbons (oil and gas) are formed by the thermal alteration of these materials. A kinetic model can describe this diagenetic process, to model the dissolution-precipitation mechanism (e.g., Luo and Vasseur, 1996;). We do not consider this mechanism in the present work, but it can easily be incorporated (Pinna et al., 2011;Carcione et al., 2011;Carcione and Avseth, 2015).
We consider two dissimilar theories of the diagenesis-cementation process, namely smectite/illite conversion based on a kinetic reaction given by the Arrhenius equation (Pytte and Reynolds, 1989), and sandstone compaction and cementation (Walderhaug, 1996;Avseth and Lehocki, 2016). The conversion from smectite to illite (clay diagenesis) with increasing depth occurs in all shales (Scotchman, 1987) and can be described by the widely accepted model proposed by Pytte and Reynolds (1989) based on a 5th-order kinetic reaction of the Arrhenius type. The result of the conversion is that the stiffnesses of the minerals composing the shale increase with depth. To study quartz cementation in sandstones, one has to estimate cement volumes mainly as a function of temperature. According to Walderhaug (1996), the first stage is precipitation, which occurs at temperatures above 60-80 • C, when the cementation starts to be effective. The process depends on the geothermal gradient and on the effective radii of the grains (and surface area), with small quartz grains producing more cement (Walderhaug, 1996, Fig. 4). Moreover, clay minerals or calcite cement covering the grains may reduce the surface area, and inhibit dissolution. Roughly, the porosity loss is equal to the volume of precipitated quartz, if the entire cement source is sourced from outside the sediment volume and a volume of water equal to the cement volume is removed from the sediment. An additional loss is accounted for compaction through the intergranular volume index (IGV) that depends on effective stress (Lander and Walderhaug, 1999). The model does not consider grain interpenetration due to intergranular pressure solution, which is assumed to be negligible, as supported by experimental evidence (Sippel, 1968;Paxton et al., 2002). However, if required, such an effect can be implemented by using the model of Weyl (1959) or that of Stephenson et al. (1992). Figures 1 and 2 show the microstructure of a shale and a sandstone, where the composition and cementation can be appreciated. Other illustrative SEM (scanning electron micrographs) images of several shales showing textures of montmorillonite (smectite) and illitic shale can be seen in Keller et al. (1986, Figs. 4-22). SEM images of several sandstones are shown in Walderhaug (2000, Fig 3), with high and minor quartz-cement, pores surrounded by almost totally quartz-cemented areas, and stylolite containing mica and detrital clay.
We assume a simple basin-evolution model with a constant sedimentation rate, heating rate and geothermal gradient. The diagenesis process starts at a given depth, pressure and temperature, with the pores filled with liquid water. Four petro-elastical models yield the seismic velocity of the rocks. The first model (Model 1) is based on an effective mineral moduli computed with the Hashin-Shtrikmann average, Krief equation to obtain the dryrock bulk and shear moduli and Gassmann equations to estimate the wet-rock moduli. Model 2 is based on the self-consistent approximation to obtain the properties of sand/clay mixtures (Gurevich and Carcione, 2000) and partially molten rocks . Model 3 is a generalization of Gassmann equation to multi-mineral media (Carcione et al., 2005), a theory also used to describe the properties gas-hydrate bearing sediments (Gei and Carcione, 2003). These three models implement Athy equation (Athy, 1930) to account for the pressure effects. The fourth approach (Model 4) is a modified version of the patchy cement model introduced by , based on the contact-cement theory (CCT) of Dvorkin and Nur (1996), assuming cement deposited at grain contacts (scheme 1 in Mavko et al., 2009, p. 257) (see also Avseth et al., 2010), and the Hertz-Mindlin theory and IGV concept to include pressure effects. Then, the Voigt-Reuss Hill averages yield the wet-rock moduli, with the self-consistent (SC) model used to obtain the properties of the mineral mixture.
2 Temperature-pressure conditions Let us assume a rock unit at depth z. The lithostatic (or confining) pressure for an average sediment densityρ is pc =ρgz, where g is the acceleration of gravity. On the other hand, the hydrostatic pore pressure is approximately p H =ρwgz, whereρw is the density of water (values ofρ = 2.5 g/cm 3 andρw = 1.04 g/cm 3 are assumed here, and g = 9.81 m/s 2 ).
Effective pressure is defined here as i.e., as differential pressure.
For a constant sediment burial rate, S, and a constant geothermal gradient, G, the temperature variation of a particular sediment volume is with a surface temperature T 0 at time t = 0, where t is deposition time and H is the heating rate. Typical values of G range from 20 to 40 o C/km, while S may range between 0.02 and 0.5 km/m.y. (m.y. = million years).

Smectite/illite conversion in shales
To evaluate the amount of smectite/illite ratio forming the shale matrix is important, since this ratio affects the density, stiffness moduli and wave velocities of the rock. Shale mineralogy may include kaolinite, montmorillonite-smectite, illite and chlorite, so the term smectite/illite as used in this study may be representative for a mixture of clay minerals (Mondol et al., 2008). The smectite/illite composite is subject to internal hydration, so its mechanical properties such as the stiffnesses can vary depending on the rock.
The conversion smectite/illite occurs in all shales with a general release of bound water into the pore space (Scotchman, 1987). Smectite dehydration implies a stiffer matrix due to the presence of more illite and therefore higher velocities. The conversion depends on temperature and sedimentation rate. A solution to this problem has been provided by Pytte and Reynolds (1989). The process is pictorially explained in Figure 3 and mathematical given in Appendix A.
Future works should consider another effect that could be important, i.e., quartz generation as a by-product of the diagenesis smectite/illite conversion. Thyberg et al. (2009) show that this is another factor to explain the velocity increase, due to micro-crystalline quartz-cementation of the rock frame.

Cementation and compaction in sandstones.
The physico-chemical process of diagenesis-cementation-compaction in sandstones is shown in Figure 4 and described mathematically in Appendix B by using the theory developed by Walderhaug (1996) for a constant heating rate. The model assumes that the quartz grains are spherical and have the same radius. Surface area decreases as a function of porosity to account for compaction and cementation. Compaction reduces quartz surface area by increasing grain contact area and by injecting matrix material into pore spaces (this last effect neglected here). Cementation can cause further reduction of surface area when quartz grains are encased by pore-filling cements.

Petro-elastical models
We consider four petro-elastical models to obtain the P-and S-wave velocities of the rocks as a function of depth, pressure and temperature. In the examples, Model 1, 2 and 3 are applied to shales and all four models to sandstone. Some of the equations, such as Voigt, Reuss and Hashin-Shtrikmann averages, are well-known in the rock-physics community. We refer to Mavko et al. (2009) for details.

Model 1. Hashin-Shtrikmann-Gassmann model.
In the first model, smectite and illite (or cement and quartz) are "mixed" by using the Hashin-Shtrikmann averages to obtain the bulk and shear moduli of the mineral composing the frame, Ks and µs (Appendix C), respectively. Then, Krief equation (Krief et al., 1990) yields the dry-rock moduli of the frame with porosity φ, where A is a dimensionless parameter (a value A = 3 is assumed here).
The effect of the pore fluid can be accounted for by using Gassmann equations (e.g., Mavko et al., 2009;Carcione, 2014). The Gassmann bulk and shear moduli are given by where and K f is the fluid modulus.
In shales, the reduction in volume due to release of bound water (see Figure 3) is a compaction effect that is modeled with Athy equation, where φ 0 is the initial porosity and β is an empirical constant (Athy, 1930;Rieke III and Chilingarian, 1974 ) (see also Appendix B) (we consider β = 0.01/MPa in this work).
In the SC approximation, the elastic moduli of an unknown effective medium have to be found implicitly. The model has been used by Gurevich and Carcione (2000) to obtain the stiffnesses of sand-clay mixtures, where the inclusions are spherical. Here, we consider spherical grains (aspect ratio γ = 1) and pores of aspect ratio γ < 1. In this case, it is n = 3, with smectite, illite and water in the case of shales, and quartz, cement and water in the case of sandstones. The proportion of phase 3 (φ 3 , the pore space) is given by Athy equation (7).
The effective bulk and shear moduli of the composite medium (K and µ), with n phases and proportion φ i , are obtained as the roots of the following system of equations where for the grains, i.e., i =1, 2 (Mavko et al., 2009, p. 187;Carcione et al., 2020) and P 3 = 1 3 T iijj and Q 3 = 1 5 (T ijij − P ) (for the pores), where T iijj and T ijij are given in Appendix A of Berryman (1980) or in page 189 of Mavko et al. (2009). If γ = 1, P 3 and Q 3 are given by equation (9). A limitation of this theory is that the inclusions are isolated, so that pore pressures are not equilibrated and the model computes high-frequency velocities.

Model 3. Generalized Gassmann (GG) theory
We also consider the composite model of Carcione et al. (2005), where the porous medium is composed of n = 2 two solids (smectite and illite in shales; quartz and cement in sandstones) and a fluid. If φ i is the fraction of the i-th solid and φ is the porosity, it is n φ i + φ = 1.
The Gassmann modulus is where β i is the fraction of solid i per unit volume of total solid. Here K i , i = 1, . . . n and K f are the solid and fluid bulk moduli, respectively, and K mi , i = 1, . . . n are the frame moduli.
A generalization of Krief's model for a multi-mineral porous medium is used to obtain the frame moduli, (Carcione et al, 2005), Similarly, where Vµ = n i=1 β i µ i is the shear Voigt average, and µs is given by equation (C.3). The dry-rock (and wet-rock) shear modulus of the composite is µm = n i µ mi .
The pressure effects are described by Athy equation (7).
4.4 Model 4. Patchy cement (PC) model. Dvorkin and Nur (1996) developed an elastic model [contact cement theory (CCT)] based on spherical grains (see also Dvorkin et al., 1999). Let us denote with subscripts 1 and 2 the properties of the grains and cement, respectively. The model assumes that initially the sandstone is a random pack of identical spherical grains with porosity near the critical one (φ 0 ≈ 0.36, a critical porosity) and an average number of grain contacts equal to C = 9, with bulk and shear moduli respectively, where C is the number of contact per grain and Sn and S t are related to the normal and shear stiffnesses, respectively, of a cemented two-grain combination. The explicit expressions of these quantities are given in Mavko et al. (2009, p. 256) and depend on the Poisson ratio of the grains ν 1 = (3K 1 − 2µ 1 )/(6K 1 + 2µ 1 ), the shear modulus µ 1 , and on the ratio of the radius of the cement layer to the grain radius: where we consider here thatφ where φp is given by equation (B.1) and φ p0 =0.05 is the maximum amount of cement at the grain contacts (see below and Figure 5). Above this value, the cement is free in the pore space with fraction φp − φ p0 . Alternatively, for a uniform distribution of the cement on the grain surface, A unified equation for α can be obtained by introducing a new parameter that allows to interpolate between the two cases represented by equations (16) and (18) (Allo, 2019; Table   1), but as we shall see in the examples, the velocities obtained with these two cases do not differ significantly from a practical viewpoint.
Finally, we interpolate between the effective high-porosity end member given by the HS upper bound and the mineral point (i.e., zero porosity) using the Voigt-Reuss-Hill average, i.e., an arithmetic average of the Voigt and Reuss moduli, to obtain the dry-rock bulk and shear moduli: where The Hill average is close in accuracy to more sophisticated techniques such as self-consistent schemes and are applicable to complex rheologies such as general anisotropy and arbitrary grain topologies (e.g., Man and Huang, 2011) and anelasticity (Picotti et al., 2018;Qadrouh et al., 2020).

Seismic velocities
The P-wave modulus and velocity are and respectively, where ρ is the composite density, given by where ρ i and ρ f are the densities of the i-th solid phase and fluid, respectively.
Similarly, the S-wave velocity is

Examples
The elastic properties of the minerals are given in Table 1. We assume G = 30 o C/km, S= 0.04 km/m.y, which give a heating rate H = 1.2 o C/m.y. and we take T 0 = 15 o C at z = 0. Figure 7 shows the relation between depth temperature and pressure, corresponding to this simple (linear) basin modeling, given by equations (1) and (2).

Shales
We consider an initial smectite/illite ratio r 0 = 0.99 and n = 5. The kinetic reaction corresponding to smectite/illite conversion assumes E = 36 kcal/mol and c = 1.217 × 10 23 / m.y. (Pytte and Reynolds, 1989). Figure 8 shows the effect of the geothermal gradient on the illite-smectite fraction, indicating that a higher value accelerates the conversion. The conversion ratio and density (a) and P-and S-wave velocities (b) of the mineral mixture, as a function of depth, are shown in Figure 9, where, in this case G = 30 o C/km. As can be seen, the elastic properties increases with depth due to the higher stiffness and density of illite compared to smectite.
Next, we consider that at z = 0, the initial porosity is φ 0 = 0.35, and β = 0.01/MPa in Athy equation (7). Then, we compute the porosity, density and wave velocities for Models 1, 2 and 3 as a function of depth, where γ = 0.17 in Model 2 (aspect ratio of the pores). Figure   10 shows the porosity and density, which have the same values for the three models, whereas Figure 11 displays the P-wave (a) and S-wave (b) velocities, where the three models yield similar results in practice. The results of Model 2 differ at lower depths, whose theory is based on the assumption of idealized geometries of the grains and pores. If we consider spheres, i.e., γ = 1, the P and S-wave velocities predicted by Model 2 are higher by approximately 0.5 km/s than those of Models 1 and 3.
To illustrate the effect of the geothermal gradient (and heating rate) on the velocities, we compare the results from two values, G = 20 and 30 o C/km (H = 0.8 and 1.2 o C/m.y.), keeping the same sedimentation rate. Figure 12 shows the velocities as a function of depth, where the black curves correspond to the higher value of the heating rate (more illite implies higher velocities). The lower velocities for 20 o C/km are due to the higher amount of smectite.

Sandstones
The dissolution and precipitation of dissolved silica are the source of cement, which is quan-  Figure 13 shows the cement fraction and porosity, with and without compaction. In the latter case, the theory holds down to approximately 4.5 km depth, below which the porosity becomes negative and the whole pore space is filled with cement.
Let us now consider two values of f , the amount of detrital quartz in the rock, i.e., f = 0.4 (0.6 feldspar) and f = 1. Figure 14 indicates that nonquartz grains inhibit the cementation and the porosity loss. Using equation (B.2) indicates that the exponent in (B.1) is proportional to f /D so that the effect of increasing the grain diameter is similar to that of decreasing the detrital quartz content, i.e., the amount of precipitated quartz cement is less in coarse-grain sandstones, because of the reduced surface area, compared to fine-grain sandstones.
The cementation only depends on the heating rate H = GS, so that different combinations of the geothermal gradient and sedimentation rate, but keeping the same heating rate, yield the same fraction of cement. However, it is not clear from equation ( and φ (fluid or pore). The pore-aspect ratio in Model 2 is taken γ = 0.108 and the bonding parameter in Model 4 is assumed to be δ = 100 MPa. Figure 16 shows the P-and S-wave velocities of the sandstones as a function of depth, corresponding to the four petro-elastical models. As can be observed, all the models are quite consistent with similar trends and values of the velocities, with Models 1 and 3 giving almost identical results. The research indicates that the use of simple models (e.g., Model 1) can be used in combination with the diagenesis-cementation process to make predictions.
However, Models 2 and 4 have additional parameters, i.e., the pore aspect ratio (γ) (Model 2), and the maximum contact cement φ p0 and bonding parameters δ (Model 4). It can be shown that taking φ p0 = 0.1 (instead of 0.05) does not affect noticeably the results.
On the other hand, the other two parameters (γ and δ) have an effect, as can be seen in Figure 17, where the velocities of Models 2 and 4 are also represented for γ = 1 (spherical pores) and δ = 40 MPa (Gangi and Carlson, 1996). While the aspect ratio highly affects the velocity (e.g., Cheng et al., 2020), the initial bonding of the grains has a lesser influence.
Next, we consider the more simple and complex models (1 and 4, respectively) and show the P-wave velocities for two values of the geothermal gradient G = 20 and 30 o C/km (H = 0.8 and 1.2 o C/m.y.) keeping the same sedimentation rate. The results are shown in Figure 18, where, as expected, a smaller temperature gradient (and heating rate) implies less diagenesis and cementation and lower velocities. Moreover, the two models predict the same trend and values of the velocities, indicating that the more simple Model 1 is suitable for predictions as Model 4, which is more complex from a mathematical point of view.
Model 4 allows to change the location of contact cement from the grain contact to a uniform distribution on the grain surface. In this case, we use equation (18) and assume G = 30 o C/km, φ p0 = 0.05 and δ = 100 MPa. Figure 19 shows the velocities when the bonding cement is deposited at grain contacts and uniformly distributed. The differences are minimal in practice.
The fact that the properties of the cement are close to those of quartz in Model 4 does not affect much the results. The CCT theory predicts that the type of cement has not a relevant influence on the the dry-rock moduli. Varying the cement stiffness one order of magnitude results in only an about 15 % increase in velocities, an effect much smaller than that of increasing cement content [ Fig. 11 in Dvorkin et al. (1994)].

Extensions of the models
The smectite/illite petro-elastical model can be extended to the anisotropic case, as in Carcione and Avseth (2015), and to the presence of overpressure by using the pore-volume balance equations developed by Carcione and Gangi (2000) for disequilibrium compaction, theories also implemented by Qin et al. (2019). Wangen (2000) developed an alternative theory of overpressure by cementation of sediments. There is no overpressure if the released water caused by the smectite/illite reaction mostly escapes from the pore space (due in part to a relatively slow sedimentation rate), which is the case in the present work.
Another type of diagenesis process that can easily be incorporated into the theory is the conversion of solid organic matter to hydrocarbons (Luo and Vasseur, 1996;Carcione et al., 2011), as mentioned in the introduction. Attenuation and velocity dispersion can be modeled with a theory based on physical principles (not developed yet, to our knowledge) or phenomenologically with the realistic assumption that cementation reduces seismic losses (Gei and Carcione, 2003), since attenuation decreases with depth and compaction. A completely different approach is the molecular dynamic modeling (Garcia and Medina, 2007), where at each step of the algorithm, cement is added at specific locations within the pores, in three different ways that model distinct origins and microgeometric features.

Conclusions
We have developed a procedure to combine diagenesis and cementation processes for the formation of shales and sandstones from non-consolidated clay and sand, respectively, to different petro-elastical models to obtain the seismic velocities as a function of depth, pressure and temperature, on the basis of a linear basin modeling. The amount of illite in shales and quartz cementation in sandstones and resulting porosity loss and compaction affect the velocities, obtained with Gassmann-like equations, models based on oblate spheroid inclusions, and theories that assume a random packing of identical spherical grains (granular media).
The proposed methods take into account in-situ properties such as the geothermal gradient, sedimentation and heating rates, elastic properties of the diagenetic material (cement), cement distribution in the pore space, pore aspect ratios, effective pressure and initial pre- The transformation of smectite to illite is probably the most important clay mineral reaction in sedimentary rocks (Perry and Hower, 1970). Pytte and Reynolds (1989) propose a model for the smectite/illite ratio r based on the nth-order Arrhenius-type reaction where E is the activation energy, R = 1.986 cal/( mol o K ) is the gas constant, c is a constant and T (t) is the absolute temperature.
The illite/smectite ratio in percent is 100 (1 − r). The solution of equation A.1 is where r 0 is the initial ratio, m = n − 1, Ei (x) is the exponential integral, and the dependence on the deposition time is given in the absolute temperature [see equation (2)]. To compute (A.2), we use the property Ei(x) = − E 1 (−x).

Appendix B Quartz cement precipitated in sandstones and compaction
It is assumed that dissolved silica is the source of cement (quartz dissolution) at stylolites or at grain contacts containing clay or mica, and diffuses short distances to the sandstone volume, so that no quartz dissolution or grain interpenetration takes place within this volume.
Following Walderhaug (1996), the amount of quartz cement (cm 3 ) precipitated in a 1 cm 3 of sandstone at time t n+1 is where φ n p is the amount of quartz cement present at time tn, φ 0 is the initial porosity, a and b are constants, which have units of mol/(cm 2 s) and • C −1 , respectively, H is the heating rate, ρq is the quartz density, Mq is the molar mass of quartz, and is the initial quartz surface area, where D is the grain diameter, f is the fraction of detrital quartz in the rock and V is a unit volume (1 m 3 if D is given in meters). Walderhaug (1996) shows examples with 65 and 50 % detrital quartzite and the rest is feldspar.
The porosity varies as (B.4) and the surface area as The effect of clay coating is to reduce the precipitation and is modeled as where C is a clay factor. C = 1 implies zero surface area and no precipitation. The clay factor affects the surface area in the same manner as the fraction of detrital quartz f and reciprocally with respect to the grain diameter D.
Equation (B.1) assumes a constant heating rate. An algorithm for a variable heating rate is gieven by Eq. (5) in Walderhaug (1996).
In the absence of matrix material, the porosity is given by which is an empirical equation of the type used by Athy (1930) ninety years ago, of the form Then, the bulk and shear (uncemented) moduli at the critical porosity φ 0 are given by where µ 1 is the shear modulus of the grains, ν 1 is the Poisson ratio of the grains and C is the average number of contacts per spherical grain.    Sandstone diagenesis and cementation. Compaction reduces quartz surface area by increasing grain contact area, as well as by the injection of matrix material into the pore space. Cementation causes surface area reduction when quartz grains are encased by pore-filling cements (modified from Lander and Walderhaug, 1999).
(contact cement) pore-lling cement 1 -φ 0 Fig. 5 Model 4. Sandstone after cementation: φ p0 is the maximum cement fraction bonded to the grains, φp is the total volume fraction of cement, and φp − φ p0 is the pore-filling cement fraction.   for two values of the pore-aspect ratio γ and bonding parameter δ, respectively.