Integrating rock physics laboratory data and modelling to improve uplift characterization methodology

Modelling changes in elastic wave velocities in sediments and rocks as they undergo burial and uplift can reveal aspects of their burial history. This has implications in the petroleum risking procedures, as aspects such as hydrocarbon maturation and reservoir quality are key ingredients. Previous modelling for such a purpose assumes constant velocity and porosity during uplift after the cessation of cementation. Laboratory data of a synthetic sandstone formed under stress, and subjected to loading and unloading cycles are used to test this assumption. The experimental P‐wave velocities show a clear increase in stress dependence during the uplift simulation. Besides, the P‐wave anisotropy transforms from negative to positive. These observations are combined to make an updated modelling workflow to mimic the experimental results better. The modelling is conducted in different phases. Loading before cementation is modelled using a triaxial granular medium theory. For the laboratory data, a novel anisotropic formulation of the patchy cement model captures the observed trends. A modified crack model incorporates the stress sensitivity seen during the unloading that simulates uplift. The results of the integrated rock physics workflow mimic the laboratory data very well. A conceptual field example using the rock physics workflow illustrates how the exclusion of the stress dependence during uplift can lead to an underestimation of the exhumation magnitude.


I N T RO D U C T I O N
The in situ characteristics of rocks at a given depth in the subsurface are a consequence of the burial history. This has implications for economic risking procedures, as aspects such as hydrocarbon maturation and reservoir fracture development causing leakage and fluid flow are defined through the relationships of time, temperature and stress. Furthermore, velocities and density change as the sediments experience different temperature and stress conditions through geological time. Bringing rocks from greater depths to their current depth may, therefore, cause them to possess different properties compared or overburden removal by tectonic forces that causes a stress release, and a decrease in temperature.
One method to characterize the burial history is to relate the measured velocity to velocity development through time.
A framework for such a modelling sequence is described in Avseth and Lehocki (2016), where three phases are described: mechanical compaction, cementation, and uplift. In the workflow by Avseth and Lehocki (2016), neither stress dependence during unloading nor anisotropy was incorporated. In the integrated rock physics modelling suggested in this paper, the framework of modelling in different stages is adopted. Alteration of the models in each phase from Avseth and Lehocki (2016) is necessary to account for the experimental observations. Holt et al. (2014) studied the effects of loading and unloading of a synthetic sandstone formed under stress. Taking the loading and unloading to simulate burial and uplift, their experiments show that the assumption of constant porosity during uplift seems to be reasonable. In contrast, the ultrasonic velocities decrease sharply and show a definite increase in stress dependence during simulated uplift, which causes the methodology that assumes no change in the velocity to underpredict the exhumation magnitude. Bredesen (2017) incorporated stress dependence of the rock during uplift. He implements a differential elastic medium (DEM) method from Avseth et al. (2014) to account for diagenetic sandstone modelling producing erroneous velocities in the Kobbe sandstone formation in the Norwegian Sea. Crack-like pore geometries are used in the DEM methodology, and he found that the inclusion of microcracks improved the velocity predictions.
In addition to the increased stress dependence, the experiments with synthetic sandstones also showed a sign change in the P-wave anisotropy during unloading. This reversal is of further interest in an interpretation setting but requires that the models in an integrated rock physics modelling sequence can handle anisotropic stress and strain fields.
The first stage in the modelling sequence is mechanical compaction. As the burial depth increases after deposition, the effective stress increases, and the porosity decreases. These factors act to stiffen the grain assemblage, and several models exist to quantify these effects on the stiffness. The models are subdivided into theoretical, heuristic and empirical (Avseth et al., 2010).
One set of theoretical models assuming various contact boundary conditions between grains are those of Digby (1981) and Walton (1987). The basis of these models is the elastic behaviour of the contact between two elastic spheres, found in Mindlin (1949). The advantage in this context of the contactbased models is that it is easier to derive expressions that relate stiffness to stress. That is not to say that the contact models are without inherent assumptions; in particular, the assumption of spherical grains cannot be expected to be met in sediments in nature.
Another suite of theoretical models is based on inclusion theory and incorporate cavities representing the pore space, which reduces the stiffness of the rock as these inclusions are more compliant than the surrounding solid (Avseth et al., 2010). Examples of inclusion-based models are described in (but are not limited to) Eshelby (1957), Walsh (1965), Kuster and Toksöz (1974), Budiansky and O'Connell (1976), Cheng (1978Cheng ( , 1993 and Hudson (1980Hudson ( , 1981. These models utilize assumptions about the pore geometries that are not necessarily representative. This means that parameters often used as input, such as the pore aspect ratio or inclusion density, are difficult to relate to observations of rock textures viewed in, for example, microscope (Avseth et al., 2010).
In this work, the model presented in Walton (1987) is used to model the effect of increased stress on the granular package, as this model allows for the use of an anisotropic stress field. The modelling of the effects of porosity loss is not contained in this model. Instead, a heuristic interpolation along a Hashin-Shtrikman lower bound from the 36% porosity point to the mineral point (0% porosity) is used to model the effects of porosity loss through compaction (Dvorkin and Nur, 1996). This interpolation is necessary, as the theory in Walton (1987) is for random dense packing, of porosity 36%. Another pragmatic alteration made is to let the final stiffness be a combination of the no-slip and slip limits in Walton (1987), referred to as a binary mix in Bachrach and Avseth (2008). In this work, a novel way to define the binary mixing factor which incorporates both anisotropy and stress dependence is included. Lander and Walderhaug (1999) present a model that gives the intergranular volume as a function of increased effective stress. This model is used to model the porosity loss in the field case.
A rise in temperature with depth accompanies an increasing burial. Quartz cementation starts once the temperature reaches 70-80 • C (Bjørlykke and Jahren, 2010). Walderhaug (1996) developed a model to quantify the cement volume as a function of temperature and time. If the amount of cement is sufficiently high to render the stiffness of the rock isotropic and stress independent, the cement models based on the work in Dvorkin and Nur (1996) can be implemented directly. Experimental results (Holt et al., 2014), however, suggest that sparsely cemented sandstones are stress-dependent and retain any anisotropy imposed by the stress field. To incorporate these observations into the modelling results, the patchy cement model from  is implemented. As the cement volume increases, the results from this model approach those found in Dvorkin and Nur (1996).
An anisotropic effective medium model incorporating cracks and pores (Fjaer, 2006) models the stiffness during uplift. This model is consistent with the physical interpretation of the processes causing the stress dependence, and further directly relates the stiffnesses to anisotropic stress and strain states. The model implementation is by inputting small changes in the strain that are estimated based on small changes in the stress using Hooke's law. This process is iterated until the rock reaches its current stress level.
Inputting a given burial history with corresponding stress and temperature development into the modelling outlined above yields an estimated stiffness for the rock at its present depth and stress level. An advantage of such an integrated modelling sequence is that the building blocks can be substituted with other models when desired. We compare our modelled velocities and P-wave anisotropy with corresponding laboratory tests performed on synthetic sandstones. Implementation of the methodology on a conceptual burial history shows how ignoring the stress dependence of the velocities during unloading causes an under-prediction of the estimated exhumation magnitude.

Method
The experimental data originate from various studies conducted by SINTEF (Holt et al., 1996;Holt et al., 2000;Holt et al., 2014), where the main aim was to quantify how stress release during coring affects the mechanical properties of a rock core and to compare this with the behaviour of the 'virgin' (un-cored) in situ rock. Although the motivation was not directly linked to uplift, the experimental results may still be interpreted with that perspective.
A typical timeline of one such test with synthetic sandstone is illustrated in Figure 1. Quartz grains were mixed with a sodium silicate solution and loaded to predefined axial and radial stresses, mimicking rock stresses in the subsurface. At the pre-defined stress level (here 15 MPa axial and 7.5 MPa radial stress, reached by following a stress path where σ z = 2σ r ), the sample was flushed with carbon dioxide, which causes amorphous silica to precipitate as cement at the grain contacts. The samples were left to settle for an hour or two before being loaded (simulating compaction of a reservoir, or burial of a cemented rock) and unloaded (mimicking re-pressurization of a depleted reservoir, or uplift). The postcementation loading and unloading stages were performed along uniaxial strain paths. For the test shown in Figure 1, unloading was stopped when the axial and radial stress became equal, close to the initial radial stress. Note that all tests were done with unsaturated samples, and the axial and radial stresses should thus be thought of as effective vertical and horizontal stresses, respectively.
In the following, the results of one such test with a weakly cemented sandstone subjected to loading and unloading after cementation will be shown and used in the modelling of the experimental data.
In a parallel sequence of tests, identically manufactured samples were unloaded directly after cementation, in order to simulate coring. In those cases, both the axial and radial stress was reduced towards zero along different stress paths.
The stiffness (and strength) of the synthetic rock can be modified by changing the amount of sodium silicate present prior to loading, and the experiments were repeated with varying amounts of cement. Minor variations in the stress paths do exist between samples of different stiffness, but this does not affect the overall observed trends. The use of silicate cement limits the strength and stiffness of synthetic rocks. In order to prepare a very stiff sample, epoxy was used instead of synthetic quartz cement. A brief description of the main results from these tests will be given below to see how variations in the degree of cementation affect the results.
Samples were approximately 75 mm long cylinders with a diameter of 38 mm. During the experiments, the axial strain was measured, while the radial strain was obtained from strain gauges attached to a cantilever system. Uniaxial strain conditions were obtained by increasing the confining pressure as the axial stress was increasing, maintaining zero radial strain. Ultrasonic transducers were embedded in the loading pistons, permitting axial P-and S-wave velocities to be obtained by pulse transmissions. In addition, radial P-wave velocities were obtained by transmission across a diameter of the sample. Only P-wave data have been analyzed in the current study.

M O D E L L I N G WO R K F L OW
The modelling of the elastic parameters as a function of stress and temperature history is split into three phases: mechanical compaction, cementation, and uplift. A flowchart illustrating the workflow is given in Figure 2. The axial and radial stress paths over the experimental time interval utilized in this work. During pre-cementation loading and the cementation interval it can be observed that σ z = 2σ r . During loading and unloading after cementation, the stress path is controlled by uniaxial strain conditions.

Mechanical compaction
Following deposition, unconsolidated sediments are subjected to increased stress as they are buried deeper. To model the effect of this increased stress on the stiffness of the granular assemblage at critical porosity, the model from Walton (1987) is used. This model is a contact-based theory, where the average strain in the medium is related to the average stress, which allows for the calculation of the effective moduli. The advantage of this model is that the general equation is valid also in anisotropic stress fields. Several assumptions enter into Walton's model. The grains are assumed to be identical, spherical, homogeneous and elastically isotropic. Furthermore, the strain is considered uniform and interpreted as the average strain, and the number of grain contacts is fixed. The general expression for stiffness of no-slip contacts is in Walton (1987) given as B (− pq I p I q ) 1/2 I j I k δ il + (− pq I p I q ) 1/2 I i I k δ jl + (− pq I p I q ) 1/2 I j I l δ ik (1) where n is the coordination number, φ rd is the porosity of random dense packing and i j are the strain values. The I i, j,k,l represent the direction cosines of the grain contacts (Bandyopadhyay, 2009). δ i j is the Kronecker delta, defined as The parameters B and C are values given by the properties of the grains, given in Walton (1987) as Flowchart for modelling sequence. The first step is to input a burial history or stress history. Calculation of the stiffnesses at the critical porosity is done using a set of triaxial expressions derived from the general equation in Walton (1987). Porosity is estimated either based on measured strain, or from the model in Lander and Walderhaug (1999). An interpolation along a Hashin-Shtrikman lower bound allows for the estimation of stiffness at other porosities than that of random dense packing. In the laboratory, the suite of available data provides an estimate for the cement volume. The model from Walderhaug (1996) gives the cement volume estimation in the field case. Direct implementation of the contact cement model with its modifications is implemented if the cement volume exceeds the stress independent cement limit. If the rock is sparsely cemented and thus stress-dependent, the patchy cement model , modified to tackle anisotropy, is used. The crack model from Fjaer (2006) in iteration with Hooke's law (see for example Fjaer et al. (2008)) models the effects of stress release. limited by an assumption of small anisotropy of the strain. The strain anisotropy is defined by Bandyopadhyay (2009)  . In Torset (2018), a solution motivated by the work in Bandyopadhyay (2009), but without the limitation of small strain anisotropy, is presented. It is demonstrated in Torset (2018) that the assumption in Bandyopadhyay (2009) is satisfactory down to a strain anisotropy value of −0.15. The equations found in Torset (2018) are used to model the effects of the increased stress in the integrated rock physics modelling sequence.
The work in Walton (1987) is limited to no-slip and slip, relating to the tangential stiffness of the two-grain contact. Mindlin (1949) shows how the tangential stiffness is related to the normal and tangential forces on the grain interface, which again are a factor of the externally applied stress. Duffaut et al. (2010) suggest from this a methodology where the effective moduli are given as an average of the no-slip and slip contacts. This approach accounts for the fact that all grains probably exist somewhere between the limits of no-slip and slip. A methodology that gives the same result is the binary mixing factor described in Bachrach and Avseth (2008). The factor from Bachrach and Avseth (2008) treats the grain assembly as a mixture of no-slip and slip contacts, instead of having all contacts at an intermediary state. The limits of Walton (1987) are in this work combined using a stress-dependent, anisotropic binary mixing factor. Mindlin (1949) gives a cuberoot dependence of the tangential stiffness on the force term. Inspired by this (but also from the nice fit it provides to the experimental data, see below), the dimensionless binary mixing factors are given a cube-root dependence on the applied axial stress.
The subscript indicates which stiffness the factor is relevant for. In this work, 1 is used for the radial direction, and 3 is used to denote the axial direction. C 11 is therefore referred to as the radial P-wave modulus, and C 33 is referred to as the axial P-wave modulus. a, b, c and d are values that are chosen to match the measured velocities. The stress levels used to define the functions are given in MPa. The way that the binary mixing factors are implemented is to make the overall stiffness dependent on the no-slip and slip values: The physical processes that cause the binary mixing factor to be stress-dependent are believed to change after the sample is cemented. So, it seems reasonable that in a sparsely cemented sample, the nature of the stress dependence of the binary mixing factor changes after cementation. This is equivalent to defining new values for a, b, c and d, which is relevant in scenarios where the sample retains stress sensitivity after cementation.
The equations derived by Walton (1987) are assuming the porosity of a random dense packing, 36%. Increased stress as the burial depth increases reduces the porosity. When considering a real burial history in the workflow, a model from Lander and Walderhaug (1999) is used. This model gives the intergranular volume of a clean, uncemented sand as In this equation, IGV stands for intergranular volume, which is the part of a rock that is between the grains. In a clean, uncemented sandstone, this is equivalent to the porosity. φ 0d is the starting (depositional) porosity, and IGV f is the stable packing configuration, which is the smallest porosity possible by mechanical compaction alone. IGV f is selected to be 28%, which is the value given for rigid grains in Lander and Walderhaug (1999). The β factor describes the stress dependence of the porosity and is, in this work, chosen as 0.06, which is the same as was given for rigid grains in Lander and Walderhaug (1999). σ represents the effective stress. If the stress field is anisotropic, the implementation should be with the principal effective stress. If the laboratory branch is followed in Figure 2, the porosity is estimated from the measured strain, by an equation given in Fjaer (2006) In this equation, φ is the porosity, φ 0s is the starting porosity and vol is the volumetric strain.
To account for different porosities than that of the random packing, an interpolation is made along a Hashin-Shtrikman lower bound from the random dense packing porosity to the mineral point (Dvorkin and Nur, 1996). There are thus two causal factors increasing the framework stiffness: The increased stress on the grain contacts, and the impact of the porosity reduction itself on the stiffness.
Hashin-Shtrikman (Hashin and Shtrikman, 1963) bounds themselves are robust bounds, and give the narrowest range of elastic moduli of a mixture of elastic materials, without more detail about the constituents' geometries (Avseth et al., 2010). A physical interpretation of choosing to interpolate along the lower bound is that the causal effect of the mechanical packing on the stiffness is the softest way to mix the mineral and the critical porosity limit.
The Hashin-Shtrikman bounds were in their original formulation limited to isotropic constituents. Parnell and Calvo-Jurado (2015) provide an anisotropic version of the Hashin-Shtrikman bounds. Their anisotropic lower bound is used with the anisotropic granular media model.

Cementation
As the burial of sediments continues, the temperature increases. Once the sediments reach a temperature of 70-80 • C, quartz cementation is expected to start (Bjørlykke and Jahren, 2010). The cementation significantly stiffens the sediments, and reduce their sensitivity to stress (Avseth et al., 2010). In a field scenario, the integrated rock physics workflow estimates the amount of cement from the model in Walderhaug (1996). The model gives the cement volume at a given time step as where a c , b c are constants that should essentially be the same for all quartzose sandstones, with values given in Walderhaug (1996) as a c = 1.98 · 10 −22 moles cm 2 s and b c = 0.022 1 deg C . T represents the temperature in Celsius at the desired time step, and M = 60.09 g mole is the molar mass of quartz. The initial porosity is given by φ oc . In this scenario, the initial porosity refers to the porosity at the onset of cement. The mechanical compaction determines this porosity. Porosity loss before cementation thus, to some extent, affects the rate of cementation. The density of quartz is ρ = 2650 kg m 3 . A 0 is the initial area available for cementation, and c c is the heating rate, which is the gradient of temperature as a function of time (burial rate multiplied by the geothermal gradient).
Along the laboratory branch in Figure 2, the cement volume is estimated based on the available data and observation. A cement volume of 1.6% is estimated based on the knowledge that the samples are sparsely cemented.
Following the estimation of the cement volume, the effect of the cement on the stiffnesses can be calculated. A common way to calculate this effect is the contact cement model found in Dvorkin and Nur (1996). The contact cement model is a high-porosity model, so for low-intermediate porosity sands losing porosity as a result of further cementation, a heuristic modification, known as the increasing cement model, can be made (Avseth et al., 2010). The increasing cement model involves joining the contact cement model at some cement volume (2-4%) to the mineral point with a Hashin-Shtrikman upper bound (Avseth et al., 2010). Another heuristic modification is to attach a Hashin-Shtrikman lower bound to the contact cement model, which is done to diagnose sands that have low-intermediate porosities, but not necessarily high cement volumes. This is known as the constant cement model (Avseth et al., 2010).
Common for all the models mentioned above is that at the onset of cementation, the rock becomes completely stress insensitive as well as isotropic. Experimental data, on the other hand, suggest that sparsely cemented sandstones can be stress-sensitive and anisotropic. The patchy cement model is used to incorporate these observations into the model sequence . This is another heuristic model, employing Hashin-Shtrikman bounds. In the patchy cement model, the starting points are stiffness estimations at 0% cement, using a model such as that in Walton (1987) and a stiffness estimation at some critical cement volume, based on the contact/increasing cement models. The rock is stress insensitive and isotropic at this critical cement volume. These points are then joined by a Hashin-Shtrikman bound.  suggest that both an upper and lower bound can be used, interpreted to represent different physical scenarios. Connected patches of cement are represented by an upper bound, and disconnected patches of cement are considered by a lower bound. An interpolation is made along the chosen bound to the estimated cement volume. At this point, another Hashin-Shtrikman bound is employed, joining the estimated stiffness at the actual cement volume and the mineral point. Another interpolation, to the actual porosity, is made along this lower bound. As the actual cement volume approaches the chosen critical cement limit, the patchy cement model approaches the contact/increasing cement model (depending on which one is used to define the critical cement point).
In the modelling of a general field case, the increasing cement model is used directly instead of the patchy cement model. The cement volume where the Hashin-Shtrikman upper bound is connected to the contact cement model is 4%. The patchy cement model is used for the modelling of the lab data, with a critical cement limit of 6%. As pointed out in , the choice of the critical cement limit is subjective.

Uplift
At some point in the life cycle of the rock, tectonic forces and surface erosion might remove overlying layers, causing a stress reduction. When considering the development of the elastic properties of a rock through time, it is essential to consider the effects such a stress release might have. Avseth and Lehocki (2016) model under the assumption that after the cessation of cementation, the porosity and velocities remain constant. Although laboratory data support that porosity changes are small, the velocities are not constant, and so it is desirable to find a model capable of modelling these laboratory observations.
The crack model from Fjaer (2006) is used to model the observed stress and anisotropy effects during the simulated uplift interval. This model relates stiffness changes to stress and strain changes, incorporating both crack formation and crack closure.
In Fjaer (2006), the full set of stiffness moduli are given, but only the ones relevant for this work are provided here. The axial and radial P-wave moduli are given as where φ is the porosity, and the Q's are constants accounting for the interaction between the cracks and the surrounding medium that depend on the Poisson's ratio, ν. The Q's are given as (Fjaer et al., 2008) The D-parameter describes the ability of fluid to flow between a crack and the surrounding pore space during the passage of a wave. D = 1 for dry rocks, while it becomes smaller than 1 for saturated cracks. D approaches 0 for cracks of very low aspect ratio at high frequencies, when pore fluid will not be able to flow so that cracks appear sealed (Fjaer et al., 2008). The crack densities are represented by ζ x and ζ z . For discshaped crack geometries, these densities are defined as the number of cracks per unit volume, times the average of the crack radius cubed (Fjaer, 2006). The subscript represents the directions of the normal of the cracks. The stress sensitive crack densities are suggested in Fjaer (2006) as In these equations, i represents strain in different directions. In some reference stress state, σ 0 i , the strain is set to 0 and the crack densities at this point are denoted as ζ 0 i . The stress sensitivity related to the action of normal stresses is denoted N, and the stress sensitivity related to shear deformation is denoted β. In addition, a term related to the maximum shear strain, is included by Fjaer (2006), which in transversely isotropic media equals z − r . The definition of η given in Fjaer (2006) is simply as a third factor controlling the stress dependence of the cracks, the others being N and β. T 0 is a parameter accounting for tensile strength, in order to prevent diverging crack densities and hence zero velocities at zero stress. In the case where the horizontal stresses and strains are assumed to be equal, it can be seen that ζ x = ζ y , given that ζ 0 x = ζ 0 y . Cracks are implemented in the model by the values of ζ x and ζ z at the stress level where the unloading starts. This is done by taking the modelled stiffnesses from the cement/granular media models. It was shown in Torset (2018) that equations (10) and (11) can be solved for ζ 0 x and ζ 0 z in the reference state as (under the assumption of equal horizontal stresses and strains up to the point of unloading such that ζ 0 x = ζ 0 y ). and where H 0 is the P-wave modulus of the isotropic quartz grains given as H 0 = K s + 4Gs 3 , where K s and G s are the bulk and shear modulus of the solid material, respectively, given in Table A1. The Q s are the same as defined previously. φ is the porosity. The C s are the stiffnesses at the starting point of unloading, obtained from granular media and cementation models.
Microcrack formation is thought to happen both during the unloading process in the laboratory (Holt et al., 1996), and during uplift in the field (Bredesen, 2017). To make the model more applicable to a setting where the strain might not be known, strain must be estimated. Strain estimation in the modelling sequence is from the generalized Hooke's law (Fjaer et al., 2008), the strain is estimated for small changes in stress. The stiffness is updated based on the estimated strain recursively over the unloading interval. This iteration methodology led to the need to compensate for creep observed in the experiment. The addition of delayed strain accomplished this. The delayed strain is a pragmatic solution, adding some delayed strain to the strain predicted by Hooke's law at each stiffness interval, such that the cumulative strain profile looks like that observed in the experiment.

Parameter selection
Due to the large number of models implemented throughout the integrated rock physics modelling, this also leaves a large number of parameters to be defined. Some of these are more robust than others, and there is room for disagreement regarding the choice of parameters.
For the modelling of the experimental data, the full range of input parameters is given in Table A1, and for the conceptual field case the parameters are given in Table A2. An overview of the symbols used throughout this paper is given in Table A3. The grain parameters are chosen to be the same in both cases. The value of the coordination number is subject to some debate, see, for example, Mavko et al. (2009). Different coordination numbers were tested in the modelling sequences. For the laboratory data, the result of using a coordination number of 9, in line with the results from Smith et al. (1929), is shown. The results of using a value of 7, more in line with García and Medina (2006), is displayed for the conceptual field example. The effect of changing the coordination number is to move the no-slip and slip limits up or down, while the trends remain the same.
For the laboratory data, the binary mixing has been made stress-dependent and anisotropic to account for the observed velocities. After the cementation, the factors are redefined to account for the fact that the value of the binary mixing factor is believed to be less stress-dependent after cementation. This is done by reducing the value of a and c by about 0.1 (based on observation of the velocities after cementation) and adjusting b and d such that the binary mixing values are continuous as the sample is cemented (at 15 MPa).
In the conceptual field case, a constant value is used (represented by a = c = 0) and it is chosen to be isotropic (represented by b = d). The appropriate values in a real field case  Figure 2 with the input parameters given in Table A1.
are likely to vary depending on the local lithologies and stress history, so calibration as done in the lab case requires that enough data are available. Note that in the laboratory data, the binary mixing factor for the axial P-wave velocity falls below 0 at around 1 MPa (as d ≈ −c). The consequence of this is that if the data are extrapolated to lower stress levels than the range of measurements, the axial P-wave velocity is predicted to fall below the slip limit. Authors such as Makse et al. (2004) show that at low confining stresses (0.1 MPa) simulated moduli do indeed fall below the slip limit, and attributes this to non-affine grain relaxation. The choice of parameters that are input into the porosity estimation model from Lander and Walderhaug (1999) is the same as used in their paper. The same is true for the parameters used in the cement estimation from Walderhaug (1996). In the laboratory, the porosity is estimated based on the measured strain. The laboratory cement volume and critical cement volume are estimated based on the available observations. For the conceptual field example, burial histories must be input into the modelling sequence. The two scenarios that will be compared is a rock buried at a rate of 50 m/Myr to a maximum depth of 2800 m before being uplifted at a rate of 100 m/Myr, where the assumption of no stress effects during uplift is used. The second case will use the same burial rate, but in this case the maximum burial depth is 3200 m. Following this, the crack model is implemented during uplift with a rate of 100 m/Myr to reach the same depth of 600 m. The effective stress for the conceptual example is assumed to increase linearly at a rate of 12 MPa/km. Factors such as the burial and uplift rate as well as the effective stress gradient are inputs that are determined by the local geological history. Figure 3 shows the development of axial P-wave velocity throughout the experiment for one silica-cemented synthetic sandstone. The loading phase from around 3 MPa up to 15 MPa represents the stage of mechanical compaction. At 15 MPa, the cement is formed, and this increases the velocity sharply, while the stress remains constant. If the sample had been stress insensitive, the velocity would remain constant during the stage of loading after cementation. The patchy cement model captures this stress sensitivity after cementation. After the sample reaches 40 MPa, unloading starts. This unloading is the stage taken to represent the uplift phase in the life cycle of a rock. The velocity decrease during the unloading is significant and is captured by the crack model. The porosity only changes 0.3 porosity units, so the assumption that the porosity is unchanged by a stress release seems to be justified. Figure 4 shows the measured P-wave anisotropy, together with the modelled P-wave anisotropy. The P-wave anisotropy  Figure 2 with the input parameters given in Table A1.

Modelling versus laboratory data
where V p r is the radial P-wave velocity, and V p z is the axial P-wave velocity. During the loading phase prior to cementation, the anisotropy becomes slightly more negative. The model can capture this decrease because of the anisotropic and stress-dependent binary mixing factors. Without this modification of the binary mixing factor, the estimated anisotropy would be close to constant (but not zero) during the loading interval, as the relationship between the axial and radial stresses is constant. In such a scenario, the only anisotropy variations in the model will arise from the Hashin-Shtrikman interpolation affecting the radial and axial P-wave velocities differently. Cementation of the sample reduces the magnitude of the anisotropy. This reduction is expected and captured in the patchy cement model. During unloading, the P-wave anisotropy is reversed. The crack model replicates the measured P-wave anisotropy to a satisfactory degree.

Results from modelling a conceptual field example
The purpose of modelling the development of elastic parameters through time is to diagnose the burial history of a rock based on its present-day properties. Assume now that a formation at 600 m depth has an estimated dry P-wave velocity of 2852 m/s. Figure 5 shows two different conceptual P-wave velocity developments as a function of burial history that yield the desired modelled velocity. The green dashed line represents such a diagnostic where the effects of unloading after cementation are ignored. The other curves represent the modelling sequence that includes the crack model. The black curve extends to a depth of 0, to illustrate how the stress dependence continues to develop as the stress is further reduced. Exclusion of the stress effects during unloading would cause a 400-m underestimation of the exhumation magnitude. This underestimation would be determined by comparing the maximum depths of the different scenarios.

Varying the cement volume in the experimental data
As mentioned in the experimental background, synthetic sandstones have been manufactured with different contents and types of cement, as well as at different forming stress levels. The ultrasonic velocities increase with the degree of cementation and with increasing forming stress. Figure 6 illustrates the variations in stress dependence as a function of material stiffness expressed through the axial P-wave velocity after cementation at the forming stress level. The stress dependence is quantified by the absolute change in P-wave velocity per 10 MPa unloading. In terms of effective vertical stress change, 10 MPa corresponds to approximately 1 km depth change. For the samples that were loaded after the cementation, the stress  Table A2. The green dashed curve is obtained by following the field branches in Figure 2, but without considering stress release effects. In the full modelling sequence, including stress release effects, the different stages are colour-coded. The mechanical compaction stage combines the models of Walton (1987), Dvorkin and Nur (1996) and Lander and Walderhaug (1999). The increasing cement model (Avseth et al., 2010) is used directly to model the effects of cementation, corresponding to following the high cement volume branch in the flow-chart in Figure 2, and this gives the dark-blue curve. As uplift starts, the effects of stress release modelled by the crack model from Fjaer (2006) and cementation are superposed to give the cyan curve. Once the rock exits the cementation domain, the only physical process in the model is the effects of stress release, displayed by the black curve. Both sets of diagnostic curves reach the desired velocity at the desired depth, but the difference in exhumation magnitude is 400 m.
dependence during the post-cementation loading is given by the burial after cementation data points. The stress sensitivity during the subsequent unloading of these samples is given by the unloading after burial and cementation data points. The arrows indicate the data points for the sample used to make Figures 3 and 4.

D I S C U S S I O N
The rock physics modelling sequence illustrated in Figure 2 was able to fit experimental data to a satisfactory degree. In particular, the crack model captured the increased stress de-pendence and decreased velocity seen during unloading. The increased stress dependence in the laboratory was attributed to micro-crack formation. In terms of the strain and stress anisotropy, the observed anisotropy reversal is then not surprising. In a physical sense, as the stress is unloaded more rapidly in the axial direction, and extension is limited to this direction, cracks are expected to form at a greater rate with normals parallel to the axial direction. This means that the axially propagating P-wave is affected to a larger extent than the radially propagating P-wave. This brings the P-wave anisotropy towards isotropy and eventually it is reversed, even though the axial stress remains larger than the radial stress. It is reasonable to expect that the effect of these microcracks on a given rock is dependent on the maximum burial, degree of cementation and the magnitude of the stress release. The maximum burial dependence can be seen in the samples unloaded directly after cementation which showed a greater stress dependence than those samples that were loaded after forming. It is important, however, to bear in mind that in the laboratory the cementation took place only at a given state of stress, whereas in the field the rock could be buried further, and the stress increase, as the cementation is ongoing.
The amount of cement depends on the temperature history. In the laboratory, samples were prepared with similar 'maximum depths' but different amounts of cement. Loading after cementation (simulated burial) causes a velocity increase, and it is intuitive that the stress sensitivity is higher the softer the sandstone is: In the lower limit, the sand is uncemented, and velocities increase due to increasing stiffness between grain contacts. In the ultimate upper limit, the sand grains would be welded to each other, and the rock exhibits no stress sensitivity unless cemented bonds break and form microcracks. When the samples were unloaded after reaching their peak stress, mimicking uplift after burial, the velocity reductions are in general larger than the corresponding velocity increases during loading, except for the endpoints. The reduction of velocity during unloading can be largely attributed to damage in the form of microcracks, but for an uncemented sand or a strongly cemented sandstone, such microcracks are not likely to form. For intermediate cementation, the difference between velocity increase during burial and velocity reduction by uplift is larger, as for the synthetic sandstone modelled in Figure 3. Note that the data for all samples confirm a significant change in P-wave velocity anisotropy during simulated uplift, in agreement with Figure 4. The volumetric strain data also confirm the hypothesis of negligible porosity changes during simulated uplift (<0.2% units change per 10 MPa axial stress reduction), as well as small changes in porosity during burial after cementation (<1% units per 10 MPa axial stress reduction).
Note that stress reduction was performed under uniaxial strain conditions until the axial and radial stress became equal, as seen in Figure 1. The observed changes in velocities cover a wide range from 20 m/s to 200 m/s per 10 MPa axial stress change. In Figure 3, one sees that the axial P-wave velocity drops non-linearly and at an increased rate as the stresses are reduced. It is hence quite clear that a more substantial uplift towards lower stresses would lead to larger velocity changes. Considering the data obtained during coring simulations, where the stresses were reduced towards zero, this is indeed confirmed. Because of increased attenuation at low stresses, caused by significant crack formation, axial velocity measurements were usually not possible below 5 MPa axial stress. Using the available data, however, the reduction in velocities caused by 'uplift' ranged from 20 m/s to 600 m/s per 10 MPa axial stress reduction for all synthetic rocks used in the analysis. Thus, an uplifted reservoir rock at shallow depth (below 1000 m) may be expected to have experienced much larger velocity changes than a similar rock at a larger depth.
The observations and interpretations from the laboratory data were used to modify existing models in order to more accurately characterize the velocity development as a function of burial history. This is extended to the field by assuming that the grain-scale microcracks also develop in the field, as suggested in Bredesen (2017). In the field, macroscopic fractures and faults might also be formed, further adding to the variability in elastic wave velocities. Cracks may also be closed during creep or as a result of sealing processes over geological time.
The samples in the laboratory are dry. Investigating the effect of fluids could be done by, for example, Gassmann fluid substitution (Gassmann, 1951) and reveals that fluid in the pore space increases the velocities and decreases the magnitude of the anisotropy. This means that the absolute change in anisotropy over the unloading is reduced, but the reversal is still observed. The ability to model the effects of simulated uplift on P-wave anisotropy provided in this updated technique offers an interesting possibility in terms of uplift characterization. The observation that the anisotropy is significantly altered is of interest in the interpretation of seismic data, as input in AVO analysis.
To model the data as well as accomplished for the experimental data requires selection of input parameters that are commonly used as 'fudge factors' (for example, the coordination number). Making the no-slip fraction anisotropic and stress-dependent also grants more freedom as to the range of results that can be replicated. This freedom is useful in the case where the stress history is known, and data points are plentiful, such that the parameters are used to create potential solutions; however, it puts a degree of ambiguity on the forward modelling process. If these parameters are not known to any degree of certainty, any predicted uplift is naturally subject to errors. Some of the input parameters required for tailoring the fit seen in Figures 3 and 4 are dependent on environmental factors (temperature, stress state, etc.), so it is not expected that the parameters used for the laboratory data (especially the values of the anisotropic, stress-dependent no-slip fractions) will be universally transferable to a field case, which is why looking at the slip and no-slip limits might be more relevant. The slip and no-slip limits are dependent on parameters such as coordination number and grain shear modulus, whose values are also subject to some debate, but they are physical quantities. The quantitative nature of any predicted uplift is wholly dependent on the parameters chosen in the input model. This is, however, not a caveat unique to the work in this paper. Acknowledging the limitations of the model is an important aspect of facilitating the proper implementation of models. We want to remind the reader that most rock physics papers claiming to 'match' measured data with a model, whether it is seismic, well or laboratory, will have made a choice regarding which parameters they use to obtain these results, this paper being no exception. If quantitative assessments are made, the robustness of these choices should be challenged, especially if the authors do not address the choice of parameters. Although 'fudging' might be necessary, it should be done with a conscious awareness of the underlying uncertainty.
To calibrate the parameters, it is relevant to consider a 'reference' basin, in which rocks and sediments considered analogous to those in a potentially uplifted basin have only undergone pure subsidence. In such a scenario, the modelling down to the onset of the uplift can be used to characterize different exhumation scenarios. This might provide an indicator of which parameters could be employed in further modelling.
It is, however, important to remember the non-uniqueness of some of these parameter combinations. The prime example is probably the coordination number and grain shear modulus, where one can be increased and the other decreased, to produce identical results.

C O N C L U S I O N
A new modelling scheme to characterize velocity development as a function of burial history was developed and compared with rock physics laboratory experiments of synthetic sandstones. The proposed modelling scheme includes a crack model that captures the observed stress sensitivity in cemented samples that are unloaded. Where existing methodologies that diagnose the velocity development considering its burial history assumes that the velocity is constant during stress unloading, the novel modelling scheme predicts a velocity decrease during the uplift stage, attributed to the formation of microcracks. Excluding the stress dependence during the uplift causes an underestimation of the exhumation magnitude.
Another benefit of the updated modelling workflow is that the P-wave anisotropy can be modelled throughout the life-cycle of a rock. Improved insight into the development of the P-wave anisotropy as a function of geological time can lead to more robust interpretation strategies in heavily uplifted areas.
The modular nature of the modelling workflow means that models at a given stage in the life cycle can be replaced by other models if desired. Further work into determining the best models to use at each stage, as well as a more consistent definition of the broad range of input parameters will improve the diagnostic of the burial history further.

AC K N OW L E D G E M E N T S
Two of the authors (Kenneth Duffaut and Rune M. Holt) wish to thank the Norwegian Research Council and the industry partners of the GAMES consortium at NTNU (Grant No. 294404) for partial financial support to this work. The authors wish to thank SINTEF's Petroleum department for sharing laboratory data.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data subject to third party restrictions.