- Email: [email protected]
A new model of marine sediment compression
Earth and Planetary Science Letters 477 (2017) 21–26
Contents lists available at ScienceDirect
Earth and Planetary Science Letters www.elsevier.com/locate/epsl
A new model of marine sediment compression Kylara Martin a,∗ , Warren Wood b a b
University of Texas Institute for Geophysics, 10100 Burnet Rd (R2200), Austin, TX 78758-4445, United States Naval Research Laboratory, 1005 Balch Blvd, Stennis Space Center, MS 39529, United States
a r t i c l e
i n f o
Article history: Received 29 March 2017 Received in revised form 31 July 2017 Accepted 3 August 2017 Available online xxxx Editor: P. Shearer Keywords: compaction dewatering consolidation seismic inversion sediment properties
a b s t r a c t Marine sediments cover two-thirds of the earth, and porosity (or void ratio) is a major controlling parameter in virtually every model of seafloor properties, including strength, sound speed, hydrology, thermal conductivity, and electrical resistivity. Our new model of void ratio (e) is based on the proportional void ratio, [e p = (e − e r )/(e 0 − e r )], where e 0 is the depositional maximum at the sea floor, and e r is the minimum residual void ratio at depth. We assume the values of e 0 and e r are inherent characteristics of the sediment type. Our model further defines the compression index C c to be the square root of the proportional void ratio (C c (e ) = (e p )1/2 ). This new formulation establishes a direct relation between void ratio and effective stress: e = (e 0 − e r )−1 [log10 (σ0 /σ ) + 2(e 0 − e r )]2 /4 + e r and exhibits several advantages over previous models that we demonstrate with compression test data from the Gulf of Mexico and Nankai Trough. © 2017 Elsevier B.V. All rights reserved.
1. Introduction In order to understand sub-seafloor processes relevant to resource exploration, fluid cycling, slope stability and hazard analysis, scientists must first model the physical properties of the sediment column. Physical models of sub-seafloor properties, including density, sound speed, thermal conductivity, electrical resistivity and others, depend significantly on porosity (or equivalently, void ratio). Knowing the stress–strain relationship for each layer in a column of sediment allows estimation of the porosity from deposition at the seafloor to deep burial. From porosity, one can then estimate numerous sediment geophysical parameters and implement existing models, such as the sediment physics model of Dvorkin et al. (1999), the thermal property models of Goto and Matsubayashi (2009), and Waite et al. (2009), the permeability model of Revil and Cathles (1999) and the resistivity relationships of Archie (1942) and Collett and Ladd (2000), to create an extensive physical model of the sediment. We present here a new formulation for the reduction of sediment porosity (sediment compression) with increasing effective stress. The result is an equation relating void ratio (e) to the log10 of the vertical effective stress (σ ), which can be used to model void ratio as a function of depth. From the void ratio, and other inputs, we could use the models mentioned above to estimate
*
Corresponding author. E-mail addresses: [email protected] (K. Martin), [email protected] (W. Wood). http://dx.doi.org/10.1016/j.epsl.2017.08.008 0012-821X/© 2017 Elsevier B.V. All rights reserved.
pressure, temperature, density, sound speed, shear wave speed, and electrical resistivity. Because these models have roots in both marine geophysical and geotechnical literature, we use in this paper both porosity (φ ), and void ratio (e), to describe the volume fraction of void space, where e = φ/(1 − φ), and φ = e /(1 + e ). The expression for the compression of marine sediment during normal consolidation used in recent literature (e.g., Long et al., 2011; Dugan, 2012) dates back to Terzaghi and Peck (1948) and is based on the change in void ratio being proportional to the base 10 logarithm of the vertical effective stress:
e 2 = e 1 − C c log10
σ
(1)
where C c is the constant of proportionality, also called the compression index, and e 1 and e 2 are the initial and final void ratios. The vertical effective stress, σ (in kPa), is the load supported by the grains, equal to the difference between the total lithostatic pressure (σ ) and the pore pressure. In Eq. (1), C c is an empirical constant, equal to the slope of the e − log10 (σ ) curve for a given sediment type. Since C c is assumed to be constant, the virgin consolidation curve in e − log10 (σ ) space is modeled as a straight line. However, it has been observed for some time (e.g., Butterfield, 1979) that the stress–strain curve, or e − log10 (σ ) curve, is concave upward (as in Fig. 1), implying that C c is not constant, but actually a function of void ratio (Long et al., 2011). The variation in C c with void ratio is particularly notable in clay-rich marine sediments where void ratio varies over a larger range than it does in sand.
22
K. Martin, W. Wood / Earth and Planetary Science Letters 477 (2017) 21–26
Nomenclature Porosity: pore volume/total volume = e /(1 + e ) Void ratio: pore volume/solid volume = φ/(1 − φ) σ Total stress: assumed to be lithostatic pressure σ Vertical effective stress: load supported by the grains (difference between lithostatic and pore pressure) σ0 , e0 , φ0 Initial stress, void ratio, porosity: depositional values (at the seafloor) σr , er , φr Residual stress, void ratio, porosity: values at maximum compression (minimum pore volume) without crushing sediment grains
φ
φp
e
ep Cc v g a b
In order to build the best possible model of sediment parameters, the modeled e − log10 (σ ) curve should fit observed data as closely as possible at all stresses. Therefore, the modeled relationship must be concave upward, and produce only porosities which are physically possible. Butterfield (1979) achieved a concave upward e − log10 (σ ) curve by setting specific volume (v = 1 + e) proportional to change in effective stress, producing the equation v = v 0 (σ )C , where v 0 is the initial void ratio, and C is a fitting parameter. Long et al. (2011) showed that this equation fit consolidation tests on samples from IODP Sites U1322 and U1324 in the Gulf of Mexico better than the traditional geomechanical models (e.g., Eq. (1)) because it better mimicked the concave upward behavior of the e − log10 (σ ) curve. However, the Butterfield (1979) relationship still produces negative porosities at high effective stresses. Since porosity cannot reach negative values, the curvature of the true stress–strain curve must be fundamentally different from that represented by the Butterfield equation.
Our goal is to derive a better model for the virgin consolidation curve by allowing C c to vary with void ratio. Long et al. (2011) suggest a linear relationship between C c and void ratio (e). Such a relationship can theoretically produce negative porosities, which we wish to avoid. To determine a better functional form for C c , we first make some of the same assumptions as these previous models: that the dominant method of porosity reduction due to increased vertical effective stress is by the re-arrangement of grains through rotation, sliding, and bending with no significant mineral precipitation, crushing, or melting. We define the depositional porosity, φ0 , as the porosity at which the sediment falls out of suspension, which Dvorkin et al. (1999) refers to as the critical porosity. The vertical effective stress at φ0 is σ0 , which is not equal to zero, but is very small. The particular value of σ0 is important for determining the starting point of a stress–strain curve, and is discussed in Section 3.2. We define φr as the residual porosity, which is the porosity at which the grains have been re-arranged to their maximum packing efficiency. Beyond this limit, the dominant method of porosity reduction is by crushing, melting, or chemical alteration of individual grains. We then define the proportional porosity, φp , as the fraction of the way from φ0 to φr , or equivalently for proportional void ratio, e p , from e 0 to e r ;
1/2
de = −C c (e )dg = − (e − e r )/(e 0 − e r )
dg
Substituting a = 1.0/(e 0 − e r ) and b = −e r /(e 0 − e r ) produces
de = −(ae + b)1/2 dg . Accumulating terms of “e” on one side yields
dg = −(ae + b)−1/2 de .
1/2
C c (e ) = (e p )1/2 = (e − e r )/(e 0 − e r )
.
g = −2(ae + b)1/2 /a + C , where C is the constant of integration. To determine a value for C , we apply a boundary condition: at minimum void ratio e = e r , making ae r + b = 0. Therefore, C = g (e r ) = g r = log10 (σr ). In other words, the constant of integration, C , is the log10 of the vertical effective stress required to reach the minimum void ratio. Substituting g r for C ,
g = g r − 2(ae + b)1/2 /a or
log10
(2)
σ = log10 (σr ) − 2(ae + b)1/2 /a
(3a)
2 a log10 σr /σ /2 − b e =a 1 σr + er log210 e= 4(e 0 − e r ) σ −1
(3b) (3c)
The initial condition of the system is at maximum void ratio (e = e 0 ), where ae 0 + b = 1, g 0 = g r − 2/a, and therefore
e=
1 4(e 0 − e r )
log10
1 4(e 0 − e r )
log210
2
σ0 + 2(e 0 − e r ) σ
Expanding and simplifying yields
e=
e p = (e − e r )/(e 0 − e r ).
As our results will show, we have found empirically that the compression index is extremely well represented by the square root of the proportional void ratio;
We offer no physical justification for this functional form; it is entirely empirical. However, it offers a significant advantage over previous forms in that it is everywhere geologically reasonable. Neither negative nor infinite void ratios are ever encountered for any stress. (Equivalently, no porosity is ever less than 0 or greater than 100%.) Eq. (2) now gives us an expression for the slope at any point along the stress–strain curve. Letting g = log10 (σ ), and putting Eq. (1) in differential form yields
Integrating both sides yields
2. Model formulation
φp = (φ − φr )/(φ0 − φr ),
Proportional Porosity = (φ − φr )/(φ0 − φr ) Proportional Void Ratio = (e − e r )/(e 0 − e r ) Compression Index (e.g., Long et al., 2011) Specific Volume (Butterfield, 1979) = log10 (σ ) = 1.0/(e 0 − e r ) = −e r /(e 0 − e r )
+ er
(4a)
σ0 σ0 + log10 + e0 σ σ
(4b)
Eqs. (3) and (4) are closed form expressions for vertical effective stress as a function of void ratio and vice versa. Only three parameters are required to span the spaces described in these equations, namely e 0 , e r and either σ0 or σr . Relating e 0 to σ0 will reduce the necessary parameters to two.
K. Martin, W. Wood / Earth and Planetary Science Letters 477 (2017) 21–26
Fig. 1. Compression tests from the Gulf of Mexico (CRS, Flemings et al., 2012) and Nankai Margin (U, Saffer et al., 2010) are shown in gray. Tests shown were chosen to span the maximum range of vertical effective stress within the two datasets. Data fall within end members of φ0 , φr = 0.65, 0.31 (silty clay) and φ0 , φr = 0.85, 0.05 (clay). The compression model (thin black curves) closely overlays the data.
3. Model validation Several recent, high-quality consolidation tests conducted on marine sediments in the Gulf of Mexico (URSA) (Flemings et al., 2012) and the Nankai Trough (Saffer et al., 2010) have provided observations of void ratio vs. the log of the vertical effective stress (e − log10 (σ ) curves) for deep-water marine sediments (Fig. 1). These observations were necessarily made ex situ on samples recovered at depth in the sediment column. Fig. 1 shows only the portion of the compression test data (thick gray curves) that represents the virgin compression curve, i.e., the stress as it would exist in situ, as defined in Flemings et al. (2012) and Saffer et al. (2010). The narrow black curves in Fig. 1 were generated by fitting Eq. (4b) to the data using an annealing algorithm (Ingber, 1989) that minimizes the L1 norm of the error in void ratio. This fit yields values of e 0 , e r and σ0 , fully constraining Eq. (4b) for each compression test. A similar fit was performed using the best previous model (Long et al., 2011). Values of the fitting parameters and mean errors are listed in the Supplementary Material. For small stress ranges, the difference in curvature between the two models produces only slightly different results. But over two or more orders of stress magnitude, the difference between our model and that from Long et al. (2011) becomes more evident. A comparison of the efficacy of the two best-fit curves is shown in the residual plots in Fig. 2. The gray points are the difference between the observations and our model. The black points are the difference between the same observations and the model of Long et al. (2011). The tests chosen for inclusion in this figure are the three that cover the largest range of vertical effective stress out of the 50 compression tests analyzed for this paper. Note that not only is the overall error smaller for our model, but also that the form of the error is roughly constant over the range of effective stress observed (remarkable for such a simple empirical formulation). A model with the form of Long et al. (2011) yields an arc, suggesting that model form is less appropriate. Note also that any line of constant slope (constant C c in Eq. (1)) on the e − log10 (σ ) graph in Fig. 1a will at some point result in physically impossible negative values for the void ratio. Eq. (3c) not
23
Fig. 2. The three compression tests out of fifty in the two data sets (Flemings et al., 2012 and Saffer et al., 2010) exhibited virgin compression spanning ranges of vertical effective stress greater than a factor of ten: GoM CRS003, Nankai U88, and Nankai U68. The points are the void ratio residual of the best fit of this model (gray) and that of a power-law fit (Butterfield, 1979, black). Note primarily that the power-law form cannot fit the observations over a broad range of vertical effective stress.
only more accurately models the observations from two disparate environments, but also ensures that void ratio ranges only between e 0 and e r , as long as σ > σ0 , which is implicit in our definition of σ0 as the minimum possible stress. 3.1. Determining maximum and minimum void ratio; e 0 , e r Although intuitively related to porosity and stress at the time of deposition, as well as the minimum porosity achievable without crushing grains, an accepted quantitative method for determining e 0 and e r across all sediments does not yet exist. They are, strictly speaking, empirical compression parameters. Where compression test data exists, e 0 and e r can be determined by fitting Eqs. (3) or (4) to the data. For more general modeling purposes, one could assume that e 0 and e r depend uniquely on the grain shape and inter-grain forces, which in turn are characteristics of the sediment type (e.g. sand, silt, or clay). Void ratio extrema (e 0 , and e r ) could then be taken from the literature as reported for a given sediment type; e.g., clays and sands exhibit φ0 /φr of 0.85/0.05 (e 0 /e r = 5.67/0.05) and 0.40/0.25 (e 0 /e r = 0.67/0.33) respectively (Manger, 1963). A value of volume fraction of clay ( f cly ) for a specific parcel of sediment can then be used to determine a minimum and residual porosity (φ0 and φr respectively) for that parcel. The simplest estimation of φ0 and φr using f cly is a weighted average,
φ0 = f cly φ0 cly + (1 − f cly )φ0 snd φr = f cly φr cly + (1 − f cly )φr snd
(5)
where, for example,
φ0 cly = 0.85,
φ0 snd = 0.40,
φr cly = 0.05,
φr snd = 0.25.
However, a non-linear relationship between φ and f cly , such as that proposed by Revil et al. (2002) may also be used. In our numerical fits, the resulting range of e 0 values corresponds to the porosity range 0.65 < φ0 < 0.85, with the siltier Nankai trough samples consistently yielding lower maximum (depositional, or initial) porosities than the more clay-rich URSA mudstones. Similarly the range of e r values produced corresponds to the porosity range 0.3 > φr > 0.05, with the siltier Nankai trough samples consistently yielding higher minimum porosities than the
24
K. Martin, W. Wood / Earth and Planetary Science Letters 477 (2017) 21–26
more clay-rich URSA samples. The correspondence of the parameter values produced by the numerical fit to the values we would expect for initial and minimum porosity in sediments of the lithologies considered in this study, suggests that this method is appropriate in ocean sediments across an array of lithologies. 3.2. Determining parameter σ0 / g 0 The value of the vertical effective stress at the instant of deposition, σ0 (or g 0 ) is somewhat more enigmatic. We know that σ0 is extremely small (<104 kPa in our numerical fits), and yet finite. Interestingly, when we plot the values of e 0 and σ0 from the numerical fits of Gulf of Mexico and Nankai Trough data (Fig. 3), the results fall remarkably near a straight line in e − log10 (σ ) space: e 0 = u ∗ log10 (σ0 ) + v, where u = −0.580, v = 3.167. Such an inverse relationship between initial void ratio and initial vertical stress at the instant of deposition (the ‘critical’ point where grains first begin to exert a stress on a framework of grains that are in contact) makes intuitive sense. The higher the porosity, the less dense the sediment is, and the lower the stress exerted by that top layer of grains. The empirical nature of our model, and the ex-situ nature of the laboratory measurements prevent us from asserting that the specific log-linear relationship exhibited in Fig. 3 is an inherent property of the in-situ physical system. None-the-less, for the purposes of using our compression model in the absence of consolidation tests on physical samples, deriving σ0 from e 0 using this median fit equation is consistent with a wide range of observations. 4. Implications An accurate relation between vertical effective stress (σ ) and porosity (φ ), that uses only intrinsic properties of the sediment (φ0 , φr ), and no tuning parameters, allows for the calculation of one from the other over all observable ranges of either. For a given sediment type (specific values of φ0 , φr ) the vertical effective stress and the porosity do not vary independently. Therefore, estimates of porosity in the subsurface can be used to constrain estimates of the vertical effective stress. Because there are many common measurements strongly affected by porosity, e.g. seismic response via density and sound speed, electrical resistivity, neutron porosity logs, etc., these common measurements can also be used to constrain estimates of the stress state of the sediment column. We focus here on the seismic/acoustic response. Because porosity is the primary parameter controlling sound speed in sediments (e.g., Wyllie et al., 1956), seismic velocity models can be used to make quantitative estimates of porosity throughout a seismic section (e.g., Yuan et al., 1994). Areas where the porosity does not decrease with depth in accordance with normal consolidation (high porosity, low velocity) are inferred to have lower effective stresses and higher pore pressures. Porosity and pore pressure estimates are improved by adding additional constraints, such as the compression behavior of cored sediments (e.g., Prieux et al., 2013) to the inversion model. By using a compression model that better fits actual compression data, one could improve these inversions to not only get estimates of porosity but also subsurface stress. 4.1. Practical application An example application of our model to geologic inversions is the estimation of porosity with depth using sonic logs or P-wave interval velocity (V p ), such as shown in Fig. 4. In this inversion (iterative forward modeling) the sonic log measurements are the data space, and φ0 and φr are the model space, all as functions of depth. To the extent that φ0 and φr represent the lithology,
we are inverting for lithology. The overpressure ratio is held constant at 0.2. Our forward model, in which sediments are deposited, stresses accumulated, and porosities calculated, yields profiles of realistic porosity that can only vary from φ0 to φr . These porosities, combined with a model of sound speed such as that used by Dvorkin et al. (1999), put realistic constraints on the velocity over any given interval, which in turn put realistic constraints on the sonic velocity at each depth. The sediment physics model of Dvorkin et al. (1999) yields V p using φ , a critical (depositional) porosity (φ0 ), effective stress (σ ) and the elastic moduli (bulk, shear and density) of both the pore fluid and grains. In our formulation we assume that the vertical effective stress (σ ) is equal to the effective stress ( P eff ) required by the Dvorkin et al. (1999) model. This formulation assumes isotropy, but the assumption is not important for the near vertically traveling sound waves of sonic logging tools, or near-offset seismic P-waves. Therefore, σ = P eff = P lith − P por , where P eff is the effective pressure (stress), P por is the fluid pressure within the pores, and P lith is the lithostatic or undrained pressure, i.e. the pressure from the entire load of overlying grains and fluid. Following existing convention (e.g., Schneider et al., 2009), we define the overpressure ratio: λ∗ = ( P por − P hyd )/( P lith − P hyd ), where Phyd is the hydrostatic pressure (the pressure due to only the overlying fluids). The overpressure ratio is zero when the grains support all of the load, and is unity when the fluids support all of the load. Thus, if we know enough about the grain type to estimate the effective grain elastic moduli, and can estimate the overpressure ratio (λ∗ ), we can calculate all the stresses above any given depth and calculate V p from φ , or equivalently φ from V p . The results of our iterative forward modeling are shown in Fig. 4. In this example, an annealing algorithm (Ingber, 1989; Sen and Stoffa, 1995) is used to iteratively update the model parameters, run the forward model, and compute the objective function. The sonic log data in Fig. 4a (open triangles) are compared with the final inverted V p model (heavy line). We assume there was no gas or gas hydrate, and φ0 was allowed to range between 0.7 and 0.85, and φr between 0.15 and 0.25. The overpressure ratio was held constant at 0.2, so the vertical effective stress increases monotonically with hydrostatic and lithostatic stress. The output includes a posterior estimate of the V p , as well as an estimate of the porosity. The observed porosities were not used in the inversion, and are presented only as a check on the final modeled porosity. The fit was achieved by annealing φ0 and φr at 10 m increments (layers) down through the section. The value of σ0 was determined from e 0 using the relationship in Fig. 3. The minimized objective function was the absolute difference (in sound speed) between each sample and the modeled profile. The layers were inverted in groups of about 10–15 layers each, staggering such that the search space on the bottom 5 layers from the previously annealed group is opened up to the original a priori bounds (melted, by analogy, from the point of view of the inversion), and re-annealed with layers below. For this inversion the layers annealed as a group were layers 1–11, 6–21, 16–31, 26–41, 36–51 and 46–53. For each forward model of the inversion, the stresses (hydrostatic, lithostatic, and effective) were calculated from the seafloor downward by accumulating sediment density for all the layers above, even those already annealed, and calculating porosity from the relation described in this study. Calculations of density and sound speed included the effects of varying density and bulk modulus of water with increasing temperature and pressure at depth. The modeling resulted in profiles not only of V p but also porosity, φ0 , φr , etc. – all the inputs to the Dvorkin et al. (1999) model. The porosity profile throughout this earth model is shown in Fig. 4c (heavy line) along with porosity log data from the same
K. Martin, W. Wood / Earth and Planetary Science Letters 477 (2017) 21–26
Fig. 3. Fifty compression tests were performed on GoM (31, Flemings et al., 2012) and Nankai Trough (19, Saffer et al., 2010) cores, and these tests were modeled using the compression model presented in this study. The estimate of e 0 vs. σ0 (in kPa) for each test is shown here. A median linear fit in log10 -linear space is also shown. The y-intercept of the linear equation occurs at σ0 = 100 = 1 kPa.
25
site (grey diamonds, not used in the inversion). Correspondence between modeled and observed porosity is good, even in the top 100 mbsf, where no sonic log data exist. The value of the parameterization used in our compaction model and inversion is demonstrated by the accuracy of the porosity estimate, even where no V p data exist. The top-most group of layers included some layers for which there was no data. The model parameters (φ0 , φr ) were free to take on any value within their a priori bounds (0.35 to 0.90, and 0.05 to 0.25 respectively). The porosity estimate, verified by actual observations not used in the inversion, is based only on the limited values allowed for φ0 and φr , and their relation to effective stress via the model described in this study. Using the limited search space, the inversion was able to find the porosity structure of the upper layers required to explain the V p in the lower layers. This predictive ability is one of the benefits of a model such as ours, which only permits geologically reasonable porosities. An additional practical use of the compaction model developed here would be to use both the porosity and sonic log observations as the data, and invert for the overpressure ratio as well as φ0 and φr . Going further, one could modify the inversion algorithm described to use seismic stacking velocities in place of sonic log velocities to estimate porosity and stress state. This requires the forward calculation of stacking velocity from interval velocity, (not necessarily invoking the approximation that stacking velocity is the same as the rms average of interval velocity). The relationships described herein also suggest that if V p , porosity, and lithology are known well enough, we can quantitatively constrain the effective stress throughout the section. This has significant implications for remotely (via reflection seismic techniques) determining zones of weakness and incipient failure in sediment columns. A Bayesian approach, where inputs and outputs are described on input as probability density functions (PDFs), will yield a posterior PDF for the effective stress profile. 5. Conclusions The empirical model described here 1) establishes a direct relationship between void ratio and vertical effective stress using only fundamental sediment properties; 2) applies over all observable void ratios and stresses; 3) can never result in negative stress, porosity, void ratio, or any other nonphysical occurrence; 4) fits data from highly controlled compression tests on actual sediments more accurately than previous models. As such, this model is ideally suited not just for relating porosity and effective stress, but also as a key relation in using measurements of sound speed to remotely quantify estimates of effective stress, over-pressure and potential slope failure. Acknowledgements This work was funded through the Naval Research Laboratory Base Program, and is patent pending. We thank Editor Peter Shearer, Tim Minshull and an anonymous reviewer for their comments and questions which helped us improve this text. UTIG Contribution number 3130. Appendix A. Supplementary material
Fig. 4. Solid lines are modeled values. Point data are measured log values. a) Sonic log data from IODP site 1324B (triangles) were binned into 10 m layers (solid line) and used in an inversion for the most likely earth model. b) The inversion model space consisted of the depositional and residual porosity (black and gray lines, respectively). Everywhere in this model, the vertical effective stress was assumed to be 0.2 of the way from hydrostatic to lithostatic. c) The porosity profile predicted by the model (heavy line) compares well with actual porosity log readings (circles) from the same hole, which were not used in the inversion. The excellent fit of the porosity estimate in the upper 100 m, where no sonic data exist, highlights the value of constraining the model space via φ0 and φr .
Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.epsl.2017.08.008. References Archie, G.E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. AIME 146, 54–62. http://dx.doi.org/10.2118/ 942054-G.
26
K. Martin, W. Wood / Earth and Planetary Science Letters 477 (2017) 21–26
Butterfield, R., 1979. A natural compression law for soils (an advance on e − log p ). Geotechnique 29, 469–480. http://dx.doi.org/10.1680/geot.1979.29.4.469. Collett, T.S., Ladd, J., 2000. Detection of gas hydrate with downhole logs and assessment of gas hydrate concentrations (saturations) and gas volumes on the Blake Ridge with electrical resistivity data. Proc. Ocean Drill. Program Sci. Results 164, 179–191. http://dx.doi.org/10.2973/odp.proc.sr.164.219.2000. Dugan, B., 2012. Petrophysical and consolidation behavior of mass transport deposits from the northern Gulf of Mexico, IODP Expedition 308. Mar. Geol. 315–318, 98–107. http://dx.doi.org/10.1016/j.margeo.2012.05.001. Dvorkin, J., Prasad, M., Sakai, A., Lavoie, D., 1999. Elasticity of marine sediments: rock physics modeling. Geophys. Res. Lett. 26, 1781–1784. Flemings, P.B., John, C., Behrmann, J., Expedition 308 Scientists, 2012. Expedition 308 synthesis: overpressure, consolidation, and slope stability on the continental slope of the Gulf of Mexico. In: Proc. IODP, vol. 308. Goto, S., Matsubayashi, O., 2009. Relations between the thermal properties and porosity of sediments in the eastern flank of the Juan de Fuca Ridge. Earth Planets Space 61, 863–870. Ingber, L., 1989. Very fast simulated re-annealing. Math. Comput. Model. 12, 967–973. Long, H., Flemings, P.B., Germaine, J.T., Saffer, D.M., 2011. Consolidation and overpressure near the seafloor in the Ursa Basin, Deepwater Gulf of Mexico. Earth Planet. Sci. Lett. 305, 11–20. http://dx.doi.org/10.1016/j.epsl.2011.02.007. Manger, G.E., 1963. Porosity and Bulk Density of Sedimentary Rocks. United States Government Printing Office, Washington, DC. Prieux, V., Brossier, R., Operto, S., Virieux, J., 2013. Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field, part 2: imaging compressive-wave and shear-wave velocities. Geophys. J. Int. 194, 1665–1681. http://dx.doi.org/10.1093/gji/ggt178.
Revil, A., Cathles, L.M., 1999. Permeability of shaly sands. Water Resour. Res. 35, 651–662. Revil, A., Grauls, D., Brevart, O., 2002. Mechanical compaction of sand/clay mixtures. J. Geophys. Res. 107, 2293. http://dx.doi.org/10.1029/2001JB000318. Saffer, D.M., Guo, J., Underwood, M.B., Likos, W., Skarbek, R.M., Song, I., Gildow, M., 2010. Data report: consolidation, permeability and fabric of sediments from the Nankai continental slope, IODP Sites C0001, C0008, and C0004. In: Kinoshita, M., Tobin, H., Ashi, J., Kimura, G., Lallemant, S., Screaton, E.J., Curewitz, D., Masago, H., Moe, K.T., the Expedition 314/315/316 Scientists (Eds.), Proc. IODP. Integrated Ocean Drilling Program. Management International, Inc., Washington, DC. Schneider, J., Flemings, P.B., Dugan, B., Long, H., Germaine, J.T., 2009. Overpressure and consolidation near the seafloor of Brazos–Trinity Basin IV, northwest deepwater Gulf of Mexico. J. Geophys. Res., Solid Earth 114, 1–13. http://dx.doi.org/ 10.1029/2008JB005922. Sen, M.K., Stoffa, P.L., 1995. Global Optimization Methods in Geophysical Inversion. Elsevier Science Publications, Netherlands. Terzaghi, K., Peck, R.B., 1948. Soil Mechanics in Engineering Practice. J. Wiley, New York. Waite, W.F., Santamarina, J.C., Cortes, D.D., Dugan, B., Espinoza, D.N., Germaine, J., Jang, J., Jung, J.W., Kneafsey, T.J., Shin, H., Soga, K., Winters, W.J., 2009. Physical properties of hydrate-bearing sediments, 1–38. http://dx.doi.org/ 10.1029/2008RG000279.Table. Wyllie, M.R.J., Gregory, A.R., Gardner, L.W., 1956. Elastic wave velocities in heterogeneous and porous media. Geophysics 21, 41–70. http://dx.doi.org/10.1190/ 1.1438217. Yuan, T., Spence, G.D., Hyndman, R.D., 1994. Seismic velocities and inferred porosities in the accretionary wedge sediments at the Cascadia margin. J. Geophys. Res. 99, 4413–4427. http://dx.doi.org/10.1029/93JB03203.
Contents lists available at ScienceDirect
Earth and Planetary Science Letters www.elsevier.com/locate/epsl
A new model of marine sediment compression Kylara Martin a,∗ , Warren Wood b a b
University of Texas Institute for Geophysics, 10100 Burnet Rd (R2200), Austin, TX 78758-4445, United States Naval Research Laboratory, 1005 Balch Blvd, Stennis Space Center, MS 39529, United States
a r t i c l e
i n f o
Article history: Received 29 March 2017 Received in revised form 31 July 2017 Accepted 3 August 2017 Available online xxxx Editor: P. Shearer Keywords: compaction dewatering consolidation seismic inversion sediment properties
a b s t r a c t Marine sediments cover two-thirds of the earth, and porosity (or void ratio) is a major controlling parameter in virtually every model of seafloor properties, including strength, sound speed, hydrology, thermal conductivity, and electrical resistivity. Our new model of void ratio (e) is based on the proportional void ratio, [e p = (e − e r )/(e 0 − e r )], where e 0 is the depositional maximum at the sea floor, and e r is the minimum residual void ratio at depth. We assume the values of e 0 and e r are inherent characteristics of the sediment type. Our model further defines the compression index C c to be the square root of the proportional void ratio (C c (e ) = (e p )1/2 ). This new formulation establishes a direct relation between void ratio and effective stress: e = (e 0 − e r )−1 [log10 (σ0 /σ ) + 2(e 0 − e r )]2 /4 + e r and exhibits several advantages over previous models that we demonstrate with compression test data from the Gulf of Mexico and Nankai Trough. © 2017 Elsevier B.V. All rights reserved.
1. Introduction In order to understand sub-seafloor processes relevant to resource exploration, fluid cycling, slope stability and hazard analysis, scientists must first model the physical properties of the sediment column. Physical models of sub-seafloor properties, including density, sound speed, thermal conductivity, electrical resistivity and others, depend significantly on porosity (or equivalently, void ratio). Knowing the stress–strain relationship for each layer in a column of sediment allows estimation of the porosity from deposition at the seafloor to deep burial. From porosity, one can then estimate numerous sediment geophysical parameters and implement existing models, such as the sediment physics model of Dvorkin et al. (1999), the thermal property models of Goto and Matsubayashi (2009), and Waite et al. (2009), the permeability model of Revil and Cathles (1999) and the resistivity relationships of Archie (1942) and Collett and Ladd (2000), to create an extensive physical model of the sediment. We present here a new formulation for the reduction of sediment porosity (sediment compression) with increasing effective stress. The result is an equation relating void ratio (e) to the log10 of the vertical effective stress (σ ), which can be used to model void ratio as a function of depth. From the void ratio, and other inputs, we could use the models mentioned above to estimate
*
Corresponding author. E-mail addresses: [email protected] (K. Martin), [email protected] (W. Wood). http://dx.doi.org/10.1016/j.epsl.2017.08.008 0012-821X/© 2017 Elsevier B.V. All rights reserved.
pressure, temperature, density, sound speed, shear wave speed, and electrical resistivity. Because these models have roots in both marine geophysical and geotechnical literature, we use in this paper both porosity (φ ), and void ratio (e), to describe the volume fraction of void space, where e = φ/(1 − φ), and φ = e /(1 + e ). The expression for the compression of marine sediment during normal consolidation used in recent literature (e.g., Long et al., 2011; Dugan, 2012) dates back to Terzaghi and Peck (1948) and is based on the change in void ratio being proportional to the base 10 logarithm of the vertical effective stress:
e 2 = e 1 − C c log10
σ
(1)
where C c is the constant of proportionality, also called the compression index, and e 1 and e 2 are the initial and final void ratios. The vertical effective stress, σ (in kPa), is the load supported by the grains, equal to the difference between the total lithostatic pressure (σ ) and the pore pressure. In Eq. (1), C c is an empirical constant, equal to the slope of the e − log10 (σ ) curve for a given sediment type. Since C c is assumed to be constant, the virgin consolidation curve in e − log10 (σ ) space is modeled as a straight line. However, it has been observed for some time (e.g., Butterfield, 1979) that the stress–strain curve, or e − log10 (σ ) curve, is concave upward (as in Fig. 1), implying that C c is not constant, but actually a function of void ratio (Long et al., 2011). The variation in C c with void ratio is particularly notable in clay-rich marine sediments where void ratio varies over a larger range than it does in sand.
22
K. Martin, W. Wood / Earth and Planetary Science Letters 477 (2017) 21–26
Nomenclature Porosity: pore volume/total volume = e /(1 + e ) Void ratio: pore volume/solid volume = φ/(1 − φ) σ Total stress: assumed to be lithostatic pressure σ Vertical effective stress: load supported by the grains (difference between lithostatic and pore pressure) σ0 , e0 , φ0 Initial stress, void ratio, porosity: depositional values (at the seafloor) σr , er , φr Residual stress, void ratio, porosity: values at maximum compression (minimum pore volume) without crushing sediment grains
φ
φp
e
ep Cc v g a b
In order to build the best possible model of sediment parameters, the modeled e − log10 (σ ) curve should fit observed data as closely as possible at all stresses. Therefore, the modeled relationship must be concave upward, and produce only porosities which are physically possible. Butterfield (1979) achieved a concave upward e − log10 (σ ) curve by setting specific volume (v = 1 + e) proportional to change in effective stress, producing the equation v = v 0 (σ )C , where v 0 is the initial void ratio, and C is a fitting parameter. Long et al. (2011) showed that this equation fit consolidation tests on samples from IODP Sites U1322 and U1324 in the Gulf of Mexico better than the traditional geomechanical models (e.g., Eq. (1)) because it better mimicked the concave upward behavior of the e − log10 (σ ) curve. However, the Butterfield (1979) relationship still produces negative porosities at high effective stresses. Since porosity cannot reach negative values, the curvature of the true stress–strain curve must be fundamentally different from that represented by the Butterfield equation.
Our goal is to derive a better model for the virgin consolidation curve by allowing C c to vary with void ratio. Long et al. (2011) suggest a linear relationship between C c and void ratio (e). Such a relationship can theoretically produce negative porosities, which we wish to avoid. To determine a better functional form for C c , we first make some of the same assumptions as these previous models: that the dominant method of porosity reduction due to increased vertical effective stress is by the re-arrangement of grains through rotation, sliding, and bending with no significant mineral precipitation, crushing, or melting. We define the depositional porosity, φ0 , as the porosity at which the sediment falls out of suspension, which Dvorkin et al. (1999) refers to as the critical porosity. The vertical effective stress at φ0 is σ0 , which is not equal to zero, but is very small. The particular value of σ0 is important for determining the starting point of a stress–strain curve, and is discussed in Section 3.2. We define φr as the residual porosity, which is the porosity at which the grains have been re-arranged to their maximum packing efficiency. Beyond this limit, the dominant method of porosity reduction is by crushing, melting, or chemical alteration of individual grains. We then define the proportional porosity, φp , as the fraction of the way from φ0 to φr , or equivalently for proportional void ratio, e p , from e 0 to e r ;
1/2
de = −C c (e )dg = − (e − e r )/(e 0 − e r )
dg
Substituting a = 1.0/(e 0 − e r ) and b = −e r /(e 0 − e r ) produces
de = −(ae + b)1/2 dg . Accumulating terms of “e” on one side yields
dg = −(ae + b)−1/2 de .
1/2
C c (e ) = (e p )1/2 = (e − e r )/(e 0 − e r )
.
g = −2(ae + b)1/2 /a + C , where C is the constant of integration. To determine a value for C , we apply a boundary condition: at minimum void ratio e = e r , making ae r + b = 0. Therefore, C = g (e r ) = g r = log10 (σr ). In other words, the constant of integration, C , is the log10 of the vertical effective stress required to reach the minimum void ratio. Substituting g r for C ,
g = g r − 2(ae + b)1/2 /a or
log10
(2)
σ = log10 (σr ) − 2(ae + b)1/2 /a
(3a)
2 a log10 σr /σ /2 − b e =a 1 σr + er log210 e= 4(e 0 − e r ) σ −1
(3b) (3c)
The initial condition of the system is at maximum void ratio (e = e 0 ), where ae 0 + b = 1, g 0 = g r − 2/a, and therefore
e=
1 4(e 0 − e r )
log10
1 4(e 0 − e r )
log210
2
σ0 + 2(e 0 − e r ) σ
Expanding and simplifying yields
e=
e p = (e − e r )/(e 0 − e r ).
As our results will show, we have found empirically that the compression index is extremely well represented by the square root of the proportional void ratio;
We offer no physical justification for this functional form; it is entirely empirical. However, it offers a significant advantage over previous forms in that it is everywhere geologically reasonable. Neither negative nor infinite void ratios are ever encountered for any stress. (Equivalently, no porosity is ever less than 0 or greater than 100%.) Eq. (2) now gives us an expression for the slope at any point along the stress–strain curve. Letting g = log10 (σ ), and putting Eq. (1) in differential form yields
Integrating both sides yields
2. Model formulation
φp = (φ − φr )/(φ0 − φr ),
Proportional Porosity = (φ − φr )/(φ0 − φr ) Proportional Void Ratio = (e − e r )/(e 0 − e r ) Compression Index (e.g., Long et al., 2011) Specific Volume (Butterfield, 1979) = log10 (σ ) = 1.0/(e 0 − e r ) = −e r /(e 0 − e r )
+ er
(4a)
σ0 σ0 + log10 + e0 σ σ
(4b)
Eqs. (3) and (4) are closed form expressions for vertical effective stress as a function of void ratio and vice versa. Only three parameters are required to span the spaces described in these equations, namely e 0 , e r and either σ0 or σr . Relating e 0 to σ0 will reduce the necessary parameters to two.
K. Martin, W. Wood / Earth and Planetary Science Letters 477 (2017) 21–26
Fig. 1. Compression tests from the Gulf of Mexico (CRS, Flemings et al., 2012) and Nankai Margin (U, Saffer et al., 2010) are shown in gray. Tests shown were chosen to span the maximum range of vertical effective stress within the two datasets. Data fall within end members of φ0 , φr = 0.65, 0.31 (silty clay) and φ0 , φr = 0.85, 0.05 (clay). The compression model (thin black curves) closely overlays the data.
3. Model validation Several recent, high-quality consolidation tests conducted on marine sediments in the Gulf of Mexico (URSA) (Flemings et al., 2012) and the Nankai Trough (Saffer et al., 2010) have provided observations of void ratio vs. the log of the vertical effective stress (e − log10 (σ ) curves) for deep-water marine sediments (Fig. 1). These observations were necessarily made ex situ on samples recovered at depth in the sediment column. Fig. 1 shows only the portion of the compression test data (thick gray curves) that represents the virgin compression curve, i.e., the stress as it would exist in situ, as defined in Flemings et al. (2012) and Saffer et al. (2010). The narrow black curves in Fig. 1 were generated by fitting Eq. (4b) to the data using an annealing algorithm (Ingber, 1989) that minimizes the L1 norm of the error in void ratio. This fit yields values of e 0 , e r and σ0 , fully constraining Eq. (4b) for each compression test. A similar fit was performed using the best previous model (Long et al., 2011). Values of the fitting parameters and mean errors are listed in the Supplementary Material. For small stress ranges, the difference in curvature between the two models produces only slightly different results. But over two or more orders of stress magnitude, the difference between our model and that from Long et al. (2011) becomes more evident. A comparison of the efficacy of the two best-fit curves is shown in the residual plots in Fig. 2. The gray points are the difference between the observations and our model. The black points are the difference between the same observations and the model of Long et al. (2011). The tests chosen for inclusion in this figure are the three that cover the largest range of vertical effective stress out of the 50 compression tests analyzed for this paper. Note that not only is the overall error smaller for our model, but also that the form of the error is roughly constant over the range of effective stress observed (remarkable for such a simple empirical formulation). A model with the form of Long et al. (2011) yields an arc, suggesting that model form is less appropriate. Note also that any line of constant slope (constant C c in Eq. (1)) on the e − log10 (σ ) graph in Fig. 1a will at some point result in physically impossible negative values for the void ratio. Eq. (3c) not
23
Fig. 2. The three compression tests out of fifty in the two data sets (Flemings et al., 2012 and Saffer et al., 2010) exhibited virgin compression spanning ranges of vertical effective stress greater than a factor of ten: GoM CRS003, Nankai U88, and Nankai U68. The points are the void ratio residual of the best fit of this model (gray) and that of a power-law fit (Butterfield, 1979, black). Note primarily that the power-law form cannot fit the observations over a broad range of vertical effective stress.
only more accurately models the observations from two disparate environments, but also ensures that void ratio ranges only between e 0 and e r , as long as σ > σ0 , which is implicit in our definition of σ0 as the minimum possible stress. 3.1. Determining maximum and minimum void ratio; e 0 , e r Although intuitively related to porosity and stress at the time of deposition, as well as the minimum porosity achievable without crushing grains, an accepted quantitative method for determining e 0 and e r across all sediments does not yet exist. They are, strictly speaking, empirical compression parameters. Where compression test data exists, e 0 and e r can be determined by fitting Eqs. (3) or (4) to the data. For more general modeling purposes, one could assume that e 0 and e r depend uniquely on the grain shape and inter-grain forces, which in turn are characteristics of the sediment type (e.g. sand, silt, or clay). Void ratio extrema (e 0 , and e r ) could then be taken from the literature as reported for a given sediment type; e.g., clays and sands exhibit φ0 /φr of 0.85/0.05 (e 0 /e r = 5.67/0.05) and 0.40/0.25 (e 0 /e r = 0.67/0.33) respectively (Manger, 1963). A value of volume fraction of clay ( f cly ) for a specific parcel of sediment can then be used to determine a minimum and residual porosity (φ0 and φr respectively) for that parcel. The simplest estimation of φ0 and φr using f cly is a weighted average,
φ0 = f cly φ0 cly + (1 − f cly )φ0 snd φr = f cly φr cly + (1 − f cly )φr snd
(5)
where, for example,
φ0 cly = 0.85,
φ0 snd = 0.40,
φr cly = 0.05,
φr snd = 0.25.
However, a non-linear relationship between φ and f cly , such as that proposed by Revil et al. (2002) may also be used. In our numerical fits, the resulting range of e 0 values corresponds to the porosity range 0.65 < φ0 < 0.85, with the siltier Nankai trough samples consistently yielding lower maximum (depositional, or initial) porosities than the more clay-rich URSA mudstones. Similarly the range of e r values produced corresponds to the porosity range 0.3 > φr > 0.05, with the siltier Nankai trough samples consistently yielding higher minimum porosities than the
24
K. Martin, W. Wood / Earth and Planetary Science Letters 477 (2017) 21–26
more clay-rich URSA samples. The correspondence of the parameter values produced by the numerical fit to the values we would expect for initial and minimum porosity in sediments of the lithologies considered in this study, suggests that this method is appropriate in ocean sediments across an array of lithologies. 3.2. Determining parameter σ0 / g 0 The value of the vertical effective stress at the instant of deposition, σ0 (or g 0 ) is somewhat more enigmatic. We know that σ0 is extremely small (<104 kPa in our numerical fits), and yet finite. Interestingly, when we plot the values of e 0 and σ0 from the numerical fits of Gulf of Mexico and Nankai Trough data (Fig. 3), the results fall remarkably near a straight line in e − log10 (σ ) space: e 0 = u ∗ log10 (σ0 ) + v, where u = −0.580, v = 3.167. Such an inverse relationship between initial void ratio and initial vertical stress at the instant of deposition (the ‘critical’ point where grains first begin to exert a stress on a framework of grains that are in contact) makes intuitive sense. The higher the porosity, the less dense the sediment is, and the lower the stress exerted by that top layer of grains. The empirical nature of our model, and the ex-situ nature of the laboratory measurements prevent us from asserting that the specific log-linear relationship exhibited in Fig. 3 is an inherent property of the in-situ physical system. None-the-less, for the purposes of using our compression model in the absence of consolidation tests on physical samples, deriving σ0 from e 0 using this median fit equation is consistent with a wide range of observations. 4. Implications An accurate relation between vertical effective stress (σ ) and porosity (φ ), that uses only intrinsic properties of the sediment (φ0 , φr ), and no tuning parameters, allows for the calculation of one from the other over all observable ranges of either. For a given sediment type (specific values of φ0 , φr ) the vertical effective stress and the porosity do not vary independently. Therefore, estimates of porosity in the subsurface can be used to constrain estimates of the vertical effective stress. Because there are many common measurements strongly affected by porosity, e.g. seismic response via density and sound speed, electrical resistivity, neutron porosity logs, etc., these common measurements can also be used to constrain estimates of the stress state of the sediment column. We focus here on the seismic/acoustic response. Because porosity is the primary parameter controlling sound speed in sediments (e.g., Wyllie et al., 1956), seismic velocity models can be used to make quantitative estimates of porosity throughout a seismic section (e.g., Yuan et al., 1994). Areas where the porosity does not decrease with depth in accordance with normal consolidation (high porosity, low velocity) are inferred to have lower effective stresses and higher pore pressures. Porosity and pore pressure estimates are improved by adding additional constraints, such as the compression behavior of cored sediments (e.g., Prieux et al., 2013) to the inversion model. By using a compression model that better fits actual compression data, one could improve these inversions to not only get estimates of porosity but also subsurface stress. 4.1. Practical application An example application of our model to geologic inversions is the estimation of porosity with depth using sonic logs or P-wave interval velocity (V p ), such as shown in Fig. 4. In this inversion (iterative forward modeling) the sonic log measurements are the data space, and φ0 and φr are the model space, all as functions of depth. To the extent that φ0 and φr represent the lithology,
we are inverting for lithology. The overpressure ratio is held constant at 0.2. Our forward model, in which sediments are deposited, stresses accumulated, and porosities calculated, yields profiles of realistic porosity that can only vary from φ0 to φr . These porosities, combined with a model of sound speed such as that used by Dvorkin et al. (1999), put realistic constraints on the velocity over any given interval, which in turn put realistic constraints on the sonic velocity at each depth. The sediment physics model of Dvorkin et al. (1999) yields V p using φ , a critical (depositional) porosity (φ0 ), effective stress (σ ) and the elastic moduli (bulk, shear and density) of both the pore fluid and grains. In our formulation we assume that the vertical effective stress (σ ) is equal to the effective stress ( P eff ) required by the Dvorkin et al. (1999) model. This formulation assumes isotropy, but the assumption is not important for the near vertically traveling sound waves of sonic logging tools, or near-offset seismic P-waves. Therefore, σ = P eff = P lith − P por , where P eff is the effective pressure (stress), P por is the fluid pressure within the pores, and P lith is the lithostatic or undrained pressure, i.e. the pressure from the entire load of overlying grains and fluid. Following existing convention (e.g., Schneider et al., 2009), we define the overpressure ratio: λ∗ = ( P por − P hyd )/( P lith − P hyd ), where Phyd is the hydrostatic pressure (the pressure due to only the overlying fluids). The overpressure ratio is zero when the grains support all of the load, and is unity when the fluids support all of the load. Thus, if we know enough about the grain type to estimate the effective grain elastic moduli, and can estimate the overpressure ratio (λ∗ ), we can calculate all the stresses above any given depth and calculate V p from φ , or equivalently φ from V p . The results of our iterative forward modeling are shown in Fig. 4. In this example, an annealing algorithm (Ingber, 1989; Sen and Stoffa, 1995) is used to iteratively update the model parameters, run the forward model, and compute the objective function. The sonic log data in Fig. 4a (open triangles) are compared with the final inverted V p model (heavy line). We assume there was no gas or gas hydrate, and φ0 was allowed to range between 0.7 and 0.85, and φr between 0.15 and 0.25. The overpressure ratio was held constant at 0.2, so the vertical effective stress increases monotonically with hydrostatic and lithostatic stress. The output includes a posterior estimate of the V p , as well as an estimate of the porosity. The observed porosities were not used in the inversion, and are presented only as a check on the final modeled porosity. The fit was achieved by annealing φ0 and φr at 10 m increments (layers) down through the section. The value of σ0 was determined from e 0 using the relationship in Fig. 3. The minimized objective function was the absolute difference (in sound speed) between each sample and the modeled profile. The layers were inverted in groups of about 10–15 layers each, staggering such that the search space on the bottom 5 layers from the previously annealed group is opened up to the original a priori bounds (melted, by analogy, from the point of view of the inversion), and re-annealed with layers below. For this inversion the layers annealed as a group were layers 1–11, 6–21, 16–31, 26–41, 36–51 and 46–53. For each forward model of the inversion, the stresses (hydrostatic, lithostatic, and effective) were calculated from the seafloor downward by accumulating sediment density for all the layers above, even those already annealed, and calculating porosity from the relation described in this study. Calculations of density and sound speed included the effects of varying density and bulk modulus of water with increasing temperature and pressure at depth. The modeling resulted in profiles not only of V p but also porosity, φ0 , φr , etc. – all the inputs to the Dvorkin et al. (1999) model. The porosity profile throughout this earth model is shown in Fig. 4c (heavy line) along with porosity log data from the same
K. Martin, W. Wood / Earth and Planetary Science Letters 477 (2017) 21–26
Fig. 3. Fifty compression tests were performed on GoM (31, Flemings et al., 2012) and Nankai Trough (19, Saffer et al., 2010) cores, and these tests were modeled using the compression model presented in this study. The estimate of e 0 vs. σ0 (in kPa) for each test is shown here. A median linear fit in log10 -linear space is also shown. The y-intercept of the linear equation occurs at σ0 = 100 = 1 kPa.
25
site (grey diamonds, not used in the inversion). Correspondence between modeled and observed porosity is good, even in the top 100 mbsf, where no sonic log data exist. The value of the parameterization used in our compaction model and inversion is demonstrated by the accuracy of the porosity estimate, even where no V p data exist. The top-most group of layers included some layers for which there was no data. The model parameters (φ0 , φr ) were free to take on any value within their a priori bounds (0.35 to 0.90, and 0.05 to 0.25 respectively). The porosity estimate, verified by actual observations not used in the inversion, is based only on the limited values allowed for φ0 and φr , and their relation to effective stress via the model described in this study. Using the limited search space, the inversion was able to find the porosity structure of the upper layers required to explain the V p in the lower layers. This predictive ability is one of the benefits of a model such as ours, which only permits geologically reasonable porosities. An additional practical use of the compaction model developed here would be to use both the porosity and sonic log observations as the data, and invert for the overpressure ratio as well as φ0 and φr . Going further, one could modify the inversion algorithm described to use seismic stacking velocities in place of sonic log velocities to estimate porosity and stress state. This requires the forward calculation of stacking velocity from interval velocity, (not necessarily invoking the approximation that stacking velocity is the same as the rms average of interval velocity). The relationships described herein also suggest that if V p , porosity, and lithology are known well enough, we can quantitatively constrain the effective stress throughout the section. This has significant implications for remotely (via reflection seismic techniques) determining zones of weakness and incipient failure in sediment columns. A Bayesian approach, where inputs and outputs are described on input as probability density functions (PDFs), will yield a posterior PDF for the effective stress profile. 5. Conclusions The empirical model described here 1) establishes a direct relationship between void ratio and vertical effective stress using only fundamental sediment properties; 2) applies over all observable void ratios and stresses; 3) can never result in negative stress, porosity, void ratio, or any other nonphysical occurrence; 4) fits data from highly controlled compression tests on actual sediments more accurately than previous models. As such, this model is ideally suited not just for relating porosity and effective stress, but also as a key relation in using measurements of sound speed to remotely quantify estimates of effective stress, over-pressure and potential slope failure. Acknowledgements This work was funded through the Naval Research Laboratory Base Program, and is patent pending. We thank Editor Peter Shearer, Tim Minshull and an anonymous reviewer for their comments and questions which helped us improve this text. UTIG Contribution number 3130. Appendix A. Supplementary material
Fig. 4. Solid lines are modeled values. Point data are measured log values. a) Sonic log data from IODP site 1324B (triangles) were binned into 10 m layers (solid line) and used in an inversion for the most likely earth model. b) The inversion model space consisted of the depositional and residual porosity (black and gray lines, respectively). Everywhere in this model, the vertical effective stress was assumed to be 0.2 of the way from hydrostatic to lithostatic. c) The porosity profile predicted by the model (heavy line) compares well with actual porosity log readings (circles) from the same hole, which were not used in the inversion. The excellent fit of the porosity estimate in the upper 100 m, where no sonic data exist, highlights the value of constraining the model space via φ0 and φr .
Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.epsl.2017.08.008. References Archie, G.E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. AIME 146, 54–62. http://dx.doi.org/10.2118/ 942054-G.
26
K. Martin, W. Wood / Earth and Planetary Science Letters 477 (2017) 21–26
Butterfield, R., 1979. A natural compression law for soils (an advance on e − log p ). Geotechnique 29, 469–480. http://dx.doi.org/10.1680/geot.1979.29.4.469. Collett, T.S., Ladd, J., 2000. Detection of gas hydrate with downhole logs and assessment of gas hydrate concentrations (saturations) and gas volumes on the Blake Ridge with electrical resistivity data. Proc. Ocean Drill. Program Sci. Results 164, 179–191. http://dx.doi.org/10.2973/odp.proc.sr.164.219.2000. Dugan, B., 2012. Petrophysical and consolidation behavior of mass transport deposits from the northern Gulf of Mexico, IODP Expedition 308. Mar. Geol. 315–318, 98–107. http://dx.doi.org/10.1016/j.margeo.2012.05.001. Dvorkin, J., Prasad, M., Sakai, A., Lavoie, D., 1999. Elasticity of marine sediments: rock physics modeling. Geophys. Res. Lett. 26, 1781–1784. Flemings, P.B., John, C., Behrmann, J., Expedition 308 Scientists, 2012. Expedition 308 synthesis: overpressure, consolidation, and slope stability on the continental slope of the Gulf of Mexico. In: Proc. IODP, vol. 308. Goto, S., Matsubayashi, O., 2009. Relations between the thermal properties and porosity of sediments in the eastern flank of the Juan de Fuca Ridge. Earth Planets Space 61, 863–870. Ingber, L., 1989. Very fast simulated re-annealing. Math. Comput. Model. 12, 967–973. Long, H., Flemings, P.B., Germaine, J.T., Saffer, D.M., 2011. Consolidation and overpressure near the seafloor in the Ursa Basin, Deepwater Gulf of Mexico. Earth Planet. Sci. Lett. 305, 11–20. http://dx.doi.org/10.1016/j.epsl.2011.02.007. Manger, G.E., 1963. Porosity and Bulk Density of Sedimentary Rocks. United States Government Printing Office, Washington, DC. Prieux, V., Brossier, R., Operto, S., Virieux, J., 2013. Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field, part 2: imaging compressive-wave and shear-wave velocities. Geophys. J. Int. 194, 1665–1681. http://dx.doi.org/10.1093/gji/ggt178.
Revil, A., Cathles, L.M., 1999. Permeability of shaly sands. Water Resour. Res. 35, 651–662. Revil, A., Grauls, D., Brevart, O., 2002. Mechanical compaction of sand/clay mixtures. J. Geophys. Res. 107, 2293. http://dx.doi.org/10.1029/2001JB000318. Saffer, D.M., Guo, J., Underwood, M.B., Likos, W., Skarbek, R.M., Song, I., Gildow, M., 2010. Data report: consolidation, permeability and fabric of sediments from the Nankai continental slope, IODP Sites C0001, C0008, and C0004. In: Kinoshita, M., Tobin, H., Ashi, J., Kimura, G., Lallemant, S., Screaton, E.J., Curewitz, D., Masago, H., Moe, K.T., the Expedition 314/315/316 Scientists (Eds.), Proc. IODP. Integrated Ocean Drilling Program. Management International, Inc., Washington, DC. Schneider, J., Flemings, P.B., Dugan, B., Long, H., Germaine, J.T., 2009. Overpressure and consolidation near the seafloor of Brazos–Trinity Basin IV, northwest deepwater Gulf of Mexico. J. Geophys. Res., Solid Earth 114, 1–13. http://dx.doi.org/ 10.1029/2008JB005922. Sen, M.K., Stoffa, P.L., 1995. Global Optimization Methods in Geophysical Inversion. Elsevier Science Publications, Netherlands. Terzaghi, K., Peck, R.B., 1948. Soil Mechanics in Engineering Practice. J. Wiley, New York. Waite, W.F., Santamarina, J.C., Cortes, D.D., Dugan, B., Espinoza, D.N., Germaine, J., Jang, J., Jung, J.W., Kneafsey, T.J., Shin, H., Soga, K., Winters, W.J., 2009. Physical properties of hydrate-bearing sediments, 1–38. http://dx.doi.org/ 10.1029/2008RG000279.Table. Wyllie, M.R.J., Gregory, A.R., Gardner, L.W., 1956. Elastic wave velocities in heterogeneous and porous media. Geophysics 21, 41–70. http://dx.doi.org/10.1190/ 1.1438217. Yuan, T., Spence, G.D., Hyndman, R.D., 1994. Seismic velocities and inferred porosities in the accretionary wedge sediments at the Cascadia margin. J. Geophys. Res. 99, 4413–4427. http://dx.doi.org/10.1029/93JB03203.