Multi-exponential model to describe pressure-dependent P- and S-wave velocities and its use to estimate the crack aspect ratio

Mihály Dobróka, Norbert Péter Szabó, Tünde Edit Dobróka,Mátyás Krisztián Baraza

a Department of Geophysics, University of Miskolc, Miskolc-Egyetemváros, 3515, Hungary

b MTA-ME Geoengineering Research Group, University of Miskolc, Miskolc-Egyetemváros, 3515, Hungary

c Research Institute of Applied Earth Sciences, University of Miskolc, Miskolc-Egyetemváros, 3515, Hungary

Keywords:Multi-exponential rock physical model Spectral inversion method (SIM)Crack aspect ratio Characteristic pressures

ABSTRACT We present new quantitative model describing the pressure dependence of acoustic P- and S-wave velocities.Assuming that a variety of individual mechanisms or defects(such as cracks,pore collapse and grain crushing) can contribute to the pressure-dependent change of the wave velocity, we order a characteristic pressure to all of them and allow a series of exponential terms in the description of the(Pand S-waves) velocity-pressure function.We estimate the parameters of the multi-exponential rock physical model in inversion procedures using laboratory measured P- and S-wave velocity data.As is known, the conventional damped least squares method gives acceptable results only when one or two individual mechanisms are assumed.Increasing the number of exponential terms leads to highly nonlinear ill-posed inverse problem.Due to this reason,we develop the spectral inversion method(SIM)in which the velocity amplitudes (the spectral lines in the characteristic pressure spectrum) are only considered as unknowns.The characteristic pressures (belonging to the velocity amplitudes) are excluded from the set of inversion unknowns, instead, they are defined in a set of fixed positions equidistantly distributed in the actual interval of the independent variable(pressure).Through this novel linear inversion method,we estimate the parameters of the multi-exponential rock physical model using laboratory measured P- and S-wave velocity data.The characteristic pressures are related to the closing pressures of cracks which are described by well-known rock mechanical relationships depending on the aspect ratio of elliptical cracks.This gives the possibility to estimate the aspect ratios in terms of the characteristic pressures.

1.Introduction

The pressure dependence of seismic/acoustic wave velocity in different types of rocks has been in the focus of research for several decades(Wyllie et al.,1958;Nur and Simmons,1969;Stacey,1976;Yu et al.,1993;Darot and Reuschlé,2000;He and Schmitt,2006;Ji et al., 2007; Sengun et al., 2011).It is well-known that the wave velocity changes more quickly at the beginning of loading phase,and later it is more slowly tending to a limiting value(Birch,1960;Yu et al.,1993;Best,1997;Singh et al.,2006).Two main ideas have been published to explain this process: the closure of microcracks(Walsh and Brace,1964; Yu et al.,1993; Best,1997; Sengun et al.,2011; Shen et al., 2020) and the closure of pores (Birch, 1960;Scholz and Kranz, 1974; Jones and Wang, 1981) in rocks placed under pressure.Experiments demonstrate that besides other factors, the type of pore fluid (Toks?z et al.,1979; Khazanehdari and McCann, 2005), the porosity and grain size (Prasad and Meissner,1992; Prasad, 2002; Yu et al., 2016) influence the scale of pressure dependence.To describe the nonlinear velocity vs.pressure relationship, exponential functions are most commonly applied(Wepfer and Christensen, 1991; Yu et al., 1993; Best, 1997; Wang et al., 2005; Singh et al., 2006).In these empirical equations, the coefficients of a regression curve fitted to the measured data are given without any explanation of the physical background.A semiempirical model was published by Ji et al.(2003, 2007, 2013), in which the pressure-dependent velocity was constructed usingsome analogy with natural phenomena such as radioactive decay,cooling, and vibration attenuation.Based on the Biot-Rayleigh theory, an exact numerical modeling method was published by Zhang et al.(2019a)to analyze the effect of crack characteristics on wave propagation.

Based on the closing mechanism of microcracks, Dobróka and Somogyi-Molnár (2012) presented a simple rock physical model to describe the pressure dependence of the P-and S-wave velocities in various rock samples (Somogyiné Molnár et al., 2015).The parameters of the model were determined through inversion of laboratory-measured acoustic data.In this single exponential model (SEM), the pressure closing almost all the cracks can easily be estimated.It is straightforward to assume that two or more kinds of defects (such as cracks, pore collapse and grain crushing)can simultaneously exist in the rock samples influencing the velocity-pressure relation.Due to this reason, the SEM should be generalized.A double exponential approximation was proposed by Saul and Lumley (2013) to describe the pressure-dependent porosity, elastic properties and propagation velocity in uncemented sediments.A rock physical model assuming two mechanisms influencing the pressure dependence of the propagation velocity of P- and S-waves was introduced by Somogyiné Molnár et al.(2019).In the present paper, the double exponential pressure dependence is extended to a multi-exponential one.The parameters of the model (characteristic pressures and velocity amplitudes) are determined through inversion of laboratorymeasured datasets.

The inversion of data depending on the combination of two or more exponential functions is frequent in natural sciences.Exponential decays are common in physics, medicine, biology,geophysics, engineering, etc.In many phenomena, several exponential processes are present simultaneously.Inversion or regression procedures are directed toward the determination of the amplitudes of the decaying processes and also their decay rates.The problem is broadly discussed in the literature (Householder,1950; Gardner et al.,1959; Prony,1795).An excellent overview of the problem is given in Istratov and Vyvenko (1999).

In the case when the data show multi-exponential nature, the inversion procedure is highly nonlinear and ill-posed.To reduce these difficulties, we introduce the linear spectral inversion method (SIM) in which the decay rates (characteristic pressures)are excluded from the set of inversion unknowns.Instead,they are defined in a set of fixed positions equidistantly distributed in the actual interval of the independent variable (the pressure).Discussing stress-dependent wave-velocities, the exponential factors contain stress and characteristic pressure instead of time and decay rate, respectively.The stability and accuracy of the new inversion procedure are analyzed in this study.The parameters of the multiexponential model (MEM) are determined employing laboratorymeasured P- and S-wave velocities.In the knowledge of the characteristic pressure, the closing pressure of the microcracks can be estimated.

In a rock mechanical consideration by Walsh(1965),the closing pressure was given as a function of the aspect ratio of spheroidal cracks.The aspect ratio distribution was derived from the dry rock compressibility curve by Morlier (1971).Cheng and Toks?z (1979)applied a joint inversion technique to determining the discrete pore aspect ratio spectrum through dry and saturated sandstone velocity data.Zimmerman (1991) extended the discussion of the problem in the framework of the effective medium theory.The method was further developed by David and Zimmerman(2012).It is shown in our discussion that the aspect ratios can also be estimated using the characteristic pressures found by the inversion based on the multi-exponential rock physical model.This gives the possibility to estimate the aspect ratios in terms of the characteristic pressures.

2.Theoretical background

A semi-empirical model containing a single exponential term in describing the pressure dependence of the acoustic wave velocity was published by Ji et al.(2003).A double exponential approximation was proposed by Saul and Lumley (2013) to describe the pressure-dependent elastic properties and wave velocity in unconsolidated sediments.A double exponential rock physical model based on two different mechanisms influencing the pressure dependence of the wave velocity of P-and S-waves was introduced by Somogyiné Molnár et al.(2019).In the following,we extend the model to a multi-exponential one.

2.1.The multi-exponential rock physical model

As a starting point, we summarize the single exponential rock physical model based on the physical phenomenon of the closure of microcracks.It is assumed that due to a dσ pressure increase,the dN change in the number of the open microcracks is proportional with the N number of open microcracks(as extensive quantity)and the pressure increase:

where the proportionality constant λ is a positive material characteristic (at increasing pressure, the number of open microcracks is decreasing).The solution gives N = N0exp( -λσ), where N0is the number of the open microcracks (in a unit volume) at σ = 0.This exponential decay law is accepted for the crack density in the literature(Zhang et al., 2019b).

Another basic assumption is that the dv infinitesimal change in the wave velocity is also proportional to the change in the number of microcracks:

where the material parameter α is also positive (at increasing pressure, dN ≤0 and dv ≥0).The solution of Eqs.(1) and (2) is

where v0is the propagation velocity at zero pressure and

According to Eq.(3) at increasing pressure, the acoustic wave velocity increases from v0and asymptotically tends to vm= v0+Δv0(at high pressure, when approximately all microcracks are closed).With these notations, the velocity-pressure function can also be written as

where the material parameter λ = -dln[Δv(σ)]/dσ can be considered as the logarithmic pressure sensitivity of the velocity change Δv(σ) = v0-v(σ) belonging to pressure σ (Dobróka and Somogyi Molnár, 2012), in which

A similar (dimensionless) material characteristic called elastic piezo sensitivity was introduced by Shapiro (2003) as θ = λ/C,where C is the compressibility of the hypothetic medium in which all the microcracks are closed.

Eq.(6) is analogous to the relaxation phenomena extensively studied in physics and chemistry following the expression φ(t) =φ0exp( -t/tc),where t is the time and tcis the characteristic time constant(relaxation time).The latter gives the duration of decrease in quantity φ from its initial value φ0to its fraction φ0/e,where e denotes the basis of natural logarithm.In other words, it gives the order of magnitude of time, along which the relaxation process runs its course.Employing this similarity, we re-write Eq.(6) as

where σc= 1/λ is the characteristic pressure representing the pressure interval in which the closing of the microcracks essentially runs its course.More formally,σcmeans the pressure at which the velocity increase is Δv(σc) = Δv0/e.Based on Eq.(7), the velocity change belonging to σ = 0 can be considered as the Δv0amplitude of the stress-dependent phase velocity change.This material parameter depends on N0(the crack density in the stress-free state).The pressure-dependent velocity in Eq.(3) can also be written as

The above single exponential rock physical model has three material parameters: phase velocity in the stress-free state v0,maximal phase velocity increment Δv0(related to the crack number N0),and characteristic pressure σc.Instead of v0,the maximal phase velocity vm= v0+Δv0can also be used as a material characteristic.

In a continuum mechanical consideration,Walsh(1965)pointed out that the closing pressure of a crack with an aspect ratio α is approximated by

where E is the Young’s modulus of the rock;υ is the Poisson’s ratio of the rock; α = b/a is the aspect ratio, in which a and b are the length and width of the cracks, respectively; and f(υ) denotes the expression defined differently for the penny-shaped and the elliptical cracks in plain strain or plain stress.This means that flatter cracks with aspect ratios of α ?1 are closed at smaller pressure compared to more spherical voids with α≈1.Through Eq.(8), we can see that by increasing the pressure,the velocity asymptotically tends to vmand reaches the limit with an accuracy of around 0.5%at σ = 5 σc.Thus, we can accept the closing pressure Pc≈5 σcas a good approximation giving the following relation:

Similar considerations were presented by Ji et al.(2003) in the relation of Pcand λ.The above equation shows that the characteristic pressure σccan be related to the aspect ratio of the defects(microcracks or pores) in the rock sample.

We call the above rock physical model containing only one defect as SEM.The model parameters (vm, Δv0, σc) can be determined in the framework of an inverse problem using velocity data measured at various pressures.Due to the terminology of inverse problem theory, the pressure-dependent velocity in Eq.(8) represents the solution of the forward problem.Due to the exponential term,this inverse problem is significantly nonlinear.The model can easily be generalized to the case when two rock physical mechanisms(or two families of randomly oriented cracks in the rock,each having a constant aspect ratio significantly different from that of the other family)exist in the rock sample.This model contains two exponential terms and is called double exponential model (DEM),while in case of three rock physical mechanisms we have the tripleexponential model (TEM), defining a highly nonlinear inverse problem.

The inversion of data depending on the combination of two or more exponential functions is frequent in physics, chemistry,biology, engineering sciences, etc.Inversion or regression procedures are directed toward the determination of the amplitudes (as Δv0i)of the decaying processes and also their decay rates(σci).The problem is broadly discussed in the international literature from the very beginning of Prony (1795) to a number of more recent works such as Householder (1950), Gardner et al.(1959), Láncos(1959), Istratov and Vyvenko (1999) and others.It is experienced that the inversion procedure using forward models containing two or more exponential terms often lead to an ill-posed problem with the non-uniqueness and instability of the solution.To avoid these difficulties, we developed a special inversion method (SIM) in which we exclude the σcicharacteristic pressures from the group of inversion unknowns retaining only the amplitudes resulting in a linear inverse problem.To find an acceptable estimate for the characteristic pressures, we define M equidistantly spaced points along with the relevant interval of characteristic pressures and order an unknown amplitude to all of them (M is limited by the number of the measurement data in defining an overdetermined inverse problem).This procedure requires defining the MEM.

In defining the MEM, we assume that we have M number of hypothetical intrinsic effects simultaneously influencing the pressure dependence of the wave propagation velocity.To each of them belongs an extensive quantity with infinitesimal changes taking the following form when the pressure has an increase of dσ:

Solving the equation,we find Ψi= Ψ0iexp( -λiσ),where Ψ0iis the value of the i-th extensive quantity in a stress-free state(σ =0).Extending the validity of Eq.(2) to all of the M hypothetical effects, we find

Eq.(12) results in the total infinitesimal velocity changes as

Assuming v0as the wave propagation velocity in the stress-free state, after integration one finds

where Δvi= αiΨ0i.With increasing pressure, the propagation velocity increases from v0to the limiting value of vm= v0+Δvi.Using this relation, Eq.(14) gives

where the characteristic pressure σci= 1/λiwas introduced.Here,vmis an upper limit to which the propagation velocity asymptotically converges.It is assumed in this model that during the uploading phase,new defects are not generated.Fortin et al.(2007)showed an investigation extended to a broad pressure interval(up to 300 MPa) that above a critical pressure (P*) large mechanical change in porosity and also in wave-velocity is observable due to pore collapse and grain displacements.Our investigations arerestricted to pressures much below P*, where the appearance of new cracks or any kind of new defects are excluded.

Eq.(15)is the required formula of the MEM giving the pressuredependent wave propagation velocity.This rock physical model has 2M + 1 material parameters: vm, Δv0i, σci(i = 1,???,M).In the framework of the SIM, the number of unknowns is J = M + 1, and the amplitudes belonging to physically irrelevant characteristic pressures(their contributions given in Eq.(15)are not fitting to the measured data)will tend to zero(the hypothetical intrinsic effects behind them diminish).

In the knowledge of propagation times measured at different pressures,the unknown model parameters can be determined using an inversion procedure.Due to the terminology of inverse problem theory, Eq.(15) represents a highly nonlinear solution of the forward problem.Special cases of the model can be defined as

2.2.Inversion algorithms

The Damped Least Squares(DLSQ)algorithm is used throughout the paper (Marquardt,1959).Depending on the forward modeling formulae two basic inversion procedures are defined.

2.2.1.Conventional DLSQ algorithm

The DLSQ algorithm is applied to determining the model parameters of the SEM, DEM and TEM.In the case of the SEM, the forward problem given in Eq.(16) contains the parameters of the model vector m = (Δv01,σc1,vm).We defined the DEM in Eq.(17)where the model parameter vector of the DLSQ inversion is m =(Δv01,Δv02,σc1,σc2,vm).A similar definition of the TEM is given in Eq.(18) with the unknown model of the DLSQ procedure m =(Δv01,Δv02,Δv03,σc1,σc2,σc3,vm).In this approach, the number of unknowns is J=2 M+1(M=1,2 or 3)and the inverse problem is significantly nonlinear because of the exponential nature of the model (Eqs.(16)-(18)).More details are given in Appendix A.

2.2.2.SIM

To avoid the nonlinearity in our inversion, we exclude the characteristic pressures from the unknowns.Instead, they are defined in M fixed positions equidistantly in the interval of[0, MΔσc]as σi= Δσc(i +1/2)(i = 0,???,M -1),where Δσcis the distance between two positions.In this approach, the modelparameter vector of the inversion is m = (Δv01,???,Δv0M,vm)resulting in J = M + 1 unknowns being linearly related to the calculated data (predictions) given in Eq.(15) for the MEM.The sequence of the (Δv01,???,Δv0M) velocity amplitudes can be considered as spectrum lines belonging to the characteristic pressures (σc1,???,σcM).Based on this analogy, we call the inversion procedure for the determination of the M velocity amplitudes SIM(more details are given in Appendix A).

There is no reason to expect that the characteristic pressures in a real inversion procedure (processing the acoustic velocity data measured on a certain rock sample) coincide with the pre-defined equidistantly sampled values of σci.Instead, when using SIM,the actual characteristic pressure σccan range between two neighbouring positions σci<σc<σc(i+1)and the inversion procedure gives spectrum lines at σciand σc(i+1).The amplitude Δv = Δvi+Δvi+1is shared between the two positions depending on the distance of σcfrom the fixed neighbours (σci,σc(i+1)).We assume the two spectral amplitudes in Eq.(19)giving the equivalent characteristic pressure as amplitude weighted mean(Eq.(20))and the equivalent velocity amplitude (Eq.(21)):

In the interpretation of neighbouring spectrum lines given by SIM,this note should be taken into account.

3.Laboratory measurement data

To test the proposed inversion algorithms,laboratory-measured acoustic velocity data were used.The first dataset contains acoustic longitudinal and transverse wave velocity data measured on a sandstone sample under uniaxial load.A cylindrical sample of 92 mm in length and 35 mm in diameter was prepared from a core sample taken from the depth interval of 2943 m and 2960 m of a hydrocarbon exploratory borehole drilled in southwest Hungary.The fine-grained sandstone sample with 2.31% porosity has a density of 2565 kg/m3.In uniaxial loading, a digitally controlled load frame (with maximum load of 300 kN) was used.For stress control, the software DION7 was used, which is an advanced software for planning and executing material tests as well as for controlling the elements of the system.The load cells include piezoelectric crystals (with eigen frequency of 1 MHz)for acoustic measurements.The wave velocities were measured at N = 36 equidistantly defined pressures in the[0,91]MPa interval using the pulse transmission technique (Toks?z et al., 1979).In the knowledge of the length of the rock sample and the arrival times,the Pand S-wave velocities were determined (dataset I).

In uniaxial stress conditions, the uniformly random orientation of the microcracks can be influenced by the mono-direction load resulting in a stress-induced anisotropy, which is not taken into account in our rock physical model.This effect is avoided under hydrostatic stress conditions.Because of this reason, we test the proposed inversion procedures by involving the P- and S-wave acoustic velocity dataset measured under hydrostatic pressure conditions by Fortin et al.(2007)on Vosges sandstone(dataset II).The ultrasonic velocities (Vp, Vs) were determined using 1 MHz transducers in the confining pressure range of [0, 110] MPa (for details,see Fortin et al.(2007)and David and Zimmerman(2012)).

4.Results

We tested the MEM combined with the SIM using at first the laboratory-measured P- and S-wave velocities of the dataset I.In our inversion examples, we apply the SIM to give an estimate for the model parameters(Δv01,???,Δv0M,vm)of the MEM as well as the conventional DLSQ to determine the model parameters of the SEM(Δv01,σc1,vm), DEM (Δv01,Δv02,σc1,σc2,vm) and TEM (Δv01,Δv02,Δv03,σc1,σc2,σc3,vm), respectively.

If we have no predefined number of the intrinsic mechanisms influencing the velocity-pressure dependence, the SIM is used.To perform it, we have to define a set of characteristic pressures that are fixed during inversion.The determination of the M number of the fixed points (each of them represents a hypothetic intrinsic mechanism) requires a balance between resolution and inversion stability, i.e.increasing M results in better resolution of thecharacteristic pressures up to the N number of the measurement data (while the over-determination is decreasing.For the sake of simplicity,we assumed an equidistantly determined set lying in the[0, 90] MPa range of the applied pressure in measurements.To define an overdetermined inverse problem, we assumed M = 30(hypothetic)mechanisms(with the same number of unknown Δv0iphase velocity amplitudes) and divided the range of [0, 90] MPa into M equal parts.The spectral lines belonging to σci(i = 1,???,M)were positioned in the middle of the Δσ = 3 MPa intervals.Taking vminto account, the J = M + 1 = 31 number of unknowns and the N = 36 number of data define a marginally overdetermined inverse problem.In the linear SIM, we used also the linear DLSQ inversion algorithm to determine the model parameters of the vector m = (Δv01,???,Δv0M,vm).We selected the starting value of the damping factor as ε2= 0?01(Eq.(A35))and reduced it in each iteration by a factor of 0.9.As a starting model, all the spectral amplitudes had the same value of 0.25.After 50 iterations (the damping factor is reduced to ε2= 5?15×10-5),the characteristic pressure spectrum of P-wave velocity data is shown in Fig.1.One can see that only four spectral lines could be given (amplitudes smaller than 0.0001 were neglected).

The fit between the measured data and those calculated on the estimated model (predictions) can be seen in Fig.2.The limiting value to which the velocity function tends is vm= 4?5807 km/s.The (fixed) characteristic pressures and the Δv0ispectral amplitudes (with their estimation errors) are given in Table 1.The equivalent characteristic pressures and amplitudes are calculated using Eqs.(20) and (21), respectively.In both the data and model space, we found accurate results: the relative data distance(introduced in Eq.(A37)of Appendix A)is D=0.051%,and the mean relative estimation error is MRE= 3.74%.

Fig.2.The fit of measurements and predictions in SIM inversion of P-wave velocity data (0 ≤σc ≤90).

The starting model is quite far from the solution,and the relative distance between the measured data and those calculated on the starting model is D0= 77.26%, which is more than 1500 times greater than the relative data distance calculated on the estimated model (D = 0.051%).To test the starting model dependence of the estimation results,we repeated the inversion with the new starting model, to which all the spectral amplitudes had the same value of 0.75(D0=234.21%).We found the same results as shown in Table 1.We performed the test using several starting models successfully with the same estimated model parameters, proving the starting model independence of the SIM.

Fig.1.Characteristic pressure spectrum in SIM inversion of P-wave velocity data(0 ≤σc ≤90) (M = 30).

To demonstrate the stability of the inversion algorithm,we repeat the inversion with two parametrizations.In the first one, the characteristic pressures are assumed to be positioned in the interval of[0,75] MPa.To define a moderately overdetermined inverse problem,we assume M = 25 internal mechanisms (this includes M = 25 spectral lines,equidistantly positioned in the middle of each Δσ = 3 MPa intervals with the J = M + 1 = 26 number of unknowns and N=36 number of data).The estimated parameters of the inversion results are given in Table 2.Introducing the relative model distance between the m1and m2vectors (Eq.(22)), we find D(m)= 0?171%concerning the model parameters found in Tables 1 and 2.

In the second test of discretization, we define a significantly overdetermined problem,with M=20 spectral lines positioned in the interval of[0,60]MPa(J=21,N=36.)The relative distance of the estimated model to that given in Table 1 is D(m)= 0?435%.The same distance is found between the equivalent model parameters(D(m)= 0.436%).One can see (inside the estimation error interval)essentially the same results in both cases.

The three SIM inversion procedures (running with different parametrizations and different levels of over-determination) lead to the same results demonstrating the stability and accuracy of the new inversion procedure.In both examples, the starting parameters were the same:the initial spectral amplitude values were 0.25 and the starting value ε2= 0?01 of the damping factor was reduced in each iteration by a factor of 0.9.Repeating these tests with various starting values of model parameters, the SIM inversion gives the same set of the estimated parameters, thus the starting model independence found in the marginally overdetermined problem of M=30 is still valid at the moderately(M=25),and also in the significantly overdetermined (M = 20) inverse problem.

In the next numerical experiments,we used the previous dataset and applied the conventional(linearized)DLSQ inversion method to determine the model parameters of the DEM and TEM.Assuming the DEM in forward modeling,we found a relative data distance(defined in Eq.(A10) of Appendix A) of D = 0.054%, and the mean relative estimation error is MRE=1.497%.The estimated values of the model parameter m = (Δv01,Δv02,σc1,σc2,vm) together with their estimation error are shown in Table 3.The distance between theparameters of this model and the equivalent model given by the SIM inversion(inTable1)isD(m)=1.882%,which showsagoodagreement.

Table 1SIM inversion results for P-wave velocity data (M = 30).

Table 2SIM inversion results of P-wave velocity data (M = 25).

Assuming the TEM in forward modeling, the linearized DLSQ inversion resulted in the relative data distance D=0.041%.On the other hand, compared to DEM results, the estimation accuracy is sufficiently decreased: 107.9% mean relative estimation error was found.Because of these reasons, the linearized DLSQ inversion based on TEM cannot be considered a stable procedure and accurate approach.By our experiences, this is also true for DLSQ inversion using a greater number of internal mechanisms.

On the contrary, we found the SIM approach based on MEM applicable involving even some tens of unknowns.Thus,because of its higher stability arising from its linear nature,SIM is favorable for general use (in interpreting pressure-dependent wave velocity data) even in the case when the number of different intrinsic mechanisms is not known.

The S-wave velocities were also measured on the rock sample under investigation at the same set of pressures.The SIM inversion results are shown in Fig.3 and Table 4.We assumed M=30 spectral lines equidistantly positioned in the middle of each 3 MPa interval inside the total characteristic pressure range of [0, 90] MPa.Both the data distance D = 0.031% and the mean relative estimation error of 1.37% show accurate parameter estimation.The limiting value to which the S-wave velocity function tends is vm= 2?8 km/s.

Reliability is also proven by Fig.4 showing the good fit between the measured and calculated data.The other parametrization with M=25 and 20 spectrum lines(discussed also in the case of P-wave data) gives essentially the same result, respectively.

In our measurements,the pressure varied within the[0,91]MPa interval.As a reasonable choice, we assumed the characteristic pressure lying in the[0,90]MPa interval.It is straightforward that in the applied measurement range,there can be information in the P- and S-wave velocity data about higher (but not too high) characteristic pressures.Due to this reason, we repeated the SIM inversion of the S-wave dataset assuming [0,180] MPa interval of the possible characteristic pressures.The result found for M=30 isshown in Table 5.The data distance D = 0.036% and the mean relative estimation error of 2.28% show accurate parameter estimation.It can be seen that the equivalent characteristic pressures and the spectral amplitudes are acceptably close to those found by assuming [0, 90] MPa interval of the possible characteristic pressures (Table 4), and the distance between the two equivalent models is D(m)= 7.387%.

Table 3DLSQ inversion results for P-wave velocity data using DEM.

Fig.3.Characteristic pressure spectrum in SIM inversion of S-wave velocity data(M = 30).

In our next test,we used SIM inversion by employing a cell-wise constant set of basis functions for discretization.We applied the forward modeling formula of the MEM with the Jacobi matrix given in Eq.(A26).In this experiment, the characteristic pressures were assumed again to lie in the total pressure interval of the measurement of[0,90]MPa.To define an overdetermined problem,the cell-wise constant set of basis functions was defined in M = 30 intervals of equal length.As a starting model, all the inversion unknowns had the same value of 0.25.After 50 SIM iterations, the characteristic pressure spectrum of P-wave velocity data shown in Fig.5 was found.We give the details of parameter estimation in Table 6.

Both the data distance D = 0.054% and the mean relative estimation error of 3.82% show accurate parameter estimation.In the cell-wise constant functions-based SIM inversion of P-velocity data,we found again four intervals with non-zero spectral densitybelonging to the same set of characteristic pressures included in Table 1.The first two spectral amplitudes given in Tables 1 and 6 are nearly the same.However, there is a difference in the last two amplitudes, and the sum of the two values is 0.6792 km/s in Tables 6 and 0.6849 km/s in Table 1,which shows the consistency of the two methods.Applying the cell-wise constant function based SIM inversion to the S-wave velocity of the dataset I,we also found stable and accurate results.

Table 4SIM inversion results for S-wave data with [0, 90] MPa interval of the possible characteristic pressures (M = 30).

Fig.4.The fit of measurements and predictions in SIM inversion of S-wave velocities of the dataset I.

The above tests were performed using the dataset I measured under uniaxial stress conditions.Unfortunately, the uniaxial load can cause to some extent a stress-induced anisotropy,not involved in our rock physical model.Because this effect is excluded under hydrostatic stress conditions, we test the SIM by involving dataset II.The dataset contains N = 23 data measured in the [0,110] MPa interval.To define an overdetermined inverse problem, we assumed the M = 15 number of unknown phase velocity amplitudes and divided the pressure range [0, 45] MPa into M equal parts.The spectral lines were positioned in the middle of 3 MPa intervals.The starting values of model parameters and the damping factor are similar to those used in the above tests.The SIM inversion of P-wave data found after 30 iterations is shown in Fig.6 and Table 7.

Fig.5.Characteristic pressure spectrum obtained by SIM inversion of P-wave velocity data using cell-wise constant basis functions.

The data distance is D=0.183%which shows a good fit between the measured and predicted data.The mean relative estimation error is 11.4%.The S-wave data of dataset II was processed in an independent SIM inversion procedure.The results are shown in Fig.7 and Table 8.The data distance is D = 0.161%, and the mean relative estimation error is 5.5%.The fit between the measured and predicted data can be seen in Fig.8 for P-waves,and in Fig.9 for Swaves.A similar result was found with [0, 90] MPa interval of the possible characteristic pressures.

Table 5SIM inversion results of S-wave data with [0,180] MPa interval of the possible characteristic pressures (M = 30).

Table 6SIM results for P-wave data using discretization by cell-wise constant functions (M = 30).

As a result of the two independent SIM procedures,the characteristic pressures differ(especiallyat higher pressures)for the P-and S-wave inversions.This is a consequence of the fact that in our phenomenological MEM, there are no tools to prescribe any interaction between theP-andS-wave data influencing thetwoinversion results.In several papers discussing the velocity-pressure relationship,the characteristic pressures as regressionparameters are found different for shear and compressional waves(Eberhart-Phillips et al.,1989;Jones,1995;Khakshar et al.,1999).In the more detailed(quasimicroscopic)models of David and Zimmerman(2011,2012)as well as MacBeth (2004), the stress dependence of the shear and bulk moduli(consequently the P-and S-wave velocities)are described by different characteristic pressures.

Fig.6.Characteristic pressure spectrum estimated by SIM inversion of P-wave velocities of the dataset II with 0 ≤σc ≤45 and M = 15.

On the other hand, analyzing the pressure dependence of the elastic properties of porous and fractured rocks, Shapiro (2003)concluded a (first-order) equality of the characteristic pressures belonging to P-and S-waves.In a mathematically sound discussion,Zaitsev et al.(2017a, b) formulated a model in which the effect of cracks with arbitrary ratios of normal and shear compliances can be described.The model was found to be capable of recovering the types of cracks from pressure-dependent wave velocities.It was shown that the characteristic pressures different for P-and S-waves are the consequence of the use of a penny-shaped crack model.This difference disappears in the case of cracks with properly chosen effects.

In agreement with Shapiro’s results and the notes of Zaitsev et al.(2017a), as another alternative, it is reasonable to analyze dataset II by requiring the equality of the P- and S-wave characteristic pressures.To do this,it is necessary to integrate the P-and S-wave velocities into a joint inversion procedure.As is known,joint inversion of two or more datasets is available only if there are common model parameters coupling them together(Vozoff,1975;Dobróka et al.,1991; Sz″ucs and Chivan,1996).Taking into account the two independent SIM results, in forward modeling, we allow the existence of two (lower and upper) significant characteristic pressures but require them to be the same for P- and S-waves(coupling).The result and its comparison with the two independent inversions (in Tables 7 and 8) are shown in Table 9.It can be seen that the value of the characteristic pressure in joint inversion(common for P- and S-waves) is between the values found in independent P-and S inversion(both in the cases of σc1and σc2).The spectral amplitudes (Δv1and Δv2) are increased because they belong to an increased number of input data.The relative data distance for the joint inversion is D = 0.18 which shows similar accuracy of data fitting.

Table 7SIM inversion results of P-wave velocities of the dataset II(measured under hydrostatic pressures)with[0,45]MPa interval of the possible characteristic pressures(M=15).

Table 8SIM inversion results for S-wave velocities of the dataset II(measured under hydrostatic pressures)with[0,45]MPa interval of the possible characteristic pressures(M=15).

Fig.7.Characteristic pressure spectrum given by SIM inversion of S-wave velocities of the dataset II with 0 ≤σc ≤45 and M = 15.

The closing mechanism of cracks is well discussed in the case of hydrostatic compression (Morlier, 1971; Sayers, 1999; David and Zimmerman, 2012; Zhang et al., 2019a).The physical meaning is clearer.If the orientation of cracks is isotropic in a stress-free state,it remains isotropic under hydrostatic stress conditions.This is not the case in uniaxial compression,the increasing stress can induce a certain kind of anisotropy of the crack distribution which is not included in the rock physical model used in our investigations.On the other hand, the SIM inversion of the datasets I and II shows a quite similar structure of the velocity spectrum.This can make the application of SIM to P- and S-wave velocity data measured in uniaxial load acceptable.

Fig.8.The fit of measurements and predictions in SIM inversion of P-wave velocities of the dataset II.

Fig.9.The fit of measurements and predictions in SIM inversion of S-wave velocities of the dataset II.

5.Discussion

We performed the inversion of P- and S-wave velocities of the dataset I using the conventional DLSQ method (containing the velocity amplitudes and also the characteristic pressures as unknown model parameters) and the linear SIM (with only the velocity amplitudes as unknowns), separately.

It was experienced that both velocity amplitudes given by DLSQ inversion combined with DEM as forward modeling are“resolved”into two spectral lines in the SIM approach.Due to the fact that in the SIM, the possible characteristic pressures σci(i = 1,???,M) are fixed in the middle of subsequent pressure intervals of 3 MPa in length, a spectral line of the rock sample lying between two fixed ones is always displayed in two neighbouring intervals.Its amplitude is shared between the two neighbours depending on the(pressure) distance measured from them.The true (or efficient)characteristic pressure of the intermediate line can be estimated as the amplitude-weighted mean of the two fixed values given in Eq.(20), while the true (efficient) amplitude is the sum of the two amplitudes given by Eq.(21).

Comparing the result of the conventional DLSQ method (using DEM in forward modeling) to that given by the SIM, we can consider them consistent.The characteristic pressure σc1= 6?26 MPa given in Table 3 is close to the(amplitude-weighted)mean of the first two characteristic pressures(6.43 MPa)included in Table 1.The sum of the first two spectral amplitudes(0.707 km/s)given in Table 1 is close to Δv1= 0?7 km/s in Table 3.Similarly, the σc2= 48?37 MPa is close to the(amplitude-weighted)mean of the last two characteristic pressures (47.51 MPa) given in Table 1.The sum of the last two spectral amplitudes in Table 1 is 0.6849 km/s,which is close to Δv2= 0?698 km/s given in Table 3.Thus,the result in Table 1 should be interpreted as the presence of a spectral line at 6.43 MPa with an amplitude of 0.707 km/s and another at 47.51 MPa with an amplitude of 0.6849 km/s,respectively.Taking into account the estimation errors, it can be stated that the results of the SIM inversion are in agreement with those given by the common DLSQ inversion with DEM in forward modeling.

The results of the above inversion of P-wave velocity data can be interpreted as the simultaneous presence of two internal mechanisms or two essentially different parameters of the same mechanism (M = 2) influencing the pressure dependence of the propagation velocity.In the latter case, we can assume that there are two families of randomly oriented cracks in the rock, each having a constant aspect ratio significantly different from that of the other family.Our SIM inversion results in two appreciably different characteristic pressures,thus,taking Eq.(10)into account,the rate of the characteristic pressures can be connected to the rate of aspect ratios:α2/α1= σc2/σc1.In the case of the dataset I, this ratio for P-waves is α2/α1= 7?39.The higher aspect ratio α2= 0?015 was determined utilizing the method of David and Zimmerman (2012).The mineralogical components of the rock sample are presented in Table 10.The effective compression and shear moduli of the grains were calculated as the Voight-Reuss-Hill average of the components moduli(Mavko et al.,2009)resulting in K0= 54.5 GPa and G0= 36.8 GPa, respectively.Thus, the smaller aspect ratio influencing the pressure dependence of the P-waves is α1= α2/7?39 = 0?002.

The SIM inversion results for S-wave velocities of the dataset I shown in Table 5 can again be interpreted as two(pressure)spectral lines.The first line has a characteristic pressure of 7.95 MPa(calculated using Eq.(20)) and a velocity amplitude of 0.195 km/s(given by Eq.(21)).Using again Eqs.(20) and (21), the second spectral line has a characteristic pressure of 40.75 MPa and a velocity amplitude of 0.347 km/s.The rate of aspect ratios is α2/α1=σc2/σc1= 5?12.Using the aspect ratio α2= 0?015 arising from the method of David and Zimmerman (2012), we find= α2/5?12 = 0?0029.The values found for P-and S-waves differ around 30%, and the arithmetic mean α1= 0?0025 can be an acceptable estimate.

Table 9Results of joint inversion of the P- and S-wave velocities of the dataset II and its comparison to the results of independent SIM inversions.

Table 10Mineral composition of the rock sample in dataset II.

In the case of the dataset II,the individual SIM inversion of the Pand S-wave data gives different characteristic pressures.The above ratio for P-wave is α2/α1= σc2/σc1= 3?92.The higher aspect ratio α2= 0?13 was given in David and Zimmerman (2012).Thus,the smaller aspect ratio influencing the pressure dependence of the P-wave is= 0?0332.For S-wave data, the characteristic pressure ratio is σc2/σc1= 4?39 resulting in= 0?0296.The two data agree within 12%, thus the mean value α1= 0?031 can be an acceptable estimate.The reason for the differences can be related to the polarization of the two waves.

The joint inversion of the dataset II gave the same characteristic pressures for P-and S-waves.The common value of the aspect ratio is α2/α1= σc2/σc1= 4?69 resulting in the lower aspect ratio α1=0?0277.This value agrees within 6.4%and 16.6%found for the lower aspect ratio in the case of the individual SIM inversion of S-and Pwave velocities,respectively.

6.Conclusions

To avoid the nonlinearity of the inverse problem in processing pressure-dependent acoustic velocities, we developed a special inversion procedure(SIM),in which we exclude the σcicharacteristic pressures from the group of inversion unknowns.Instead,we define M equidistantly spaced points along with the relevant interval of characteristic pressures and order an unknown amplitude to all of them.This procedure requires defining the MEM for describing the pressure dependence of acoustic P-and S-wave phase velocities.The proposed analytical model is based on the idea that to all of the M characteristic pressure points belong a contribution to the pressure dependence of seismic/acoustic velocities of the form Δv0iexp(-σ/σci) containing the Δv0iunknown amplitudes.Depending on the number of individual terms, we introduced the SEM, DEM, or even MEM and gave the forward modeling formulae governing the pressure dependence of the wave propagationvelocity.

We tested the rock physical models on laboratory-measured datasets using P- and S-wave velocities (dataset I) measured on a fine-grained sandstone sample under uniaxial stress conditions.To perform our investigations on acoustic velocity data measured under hydrostatic pressure conditions,dataset II from the literature was also used.We processed the measured datasets by inversion techniques.To characterize the accuracy of the inversion procedures,we used the data distance and the parameter estimation error.The conventional DLSQ was applied to determining the model parameters (including also the characteristic pressures σci) in the case of SEM and DEM models.We found that the inversion procedure based on DEM is stable and accurate enough.On the other hand,by increasing the number of the exponential terms with the unknown characteristic pressures,the nonlinearity of the problem quickly reduces the accuracy of the inversion procedure.Due to this reason, we suggest the use of the SIM, an innovative seriesexpansion based inversion algorithm, in which the unknown velocity amplitudes in the forward modeling formulae for MEM was considered as spectral lines of the characteristic pressure spectrum.The number of spectral amplitudes (playing the role of expansion coefficients) was chosen to ensure the overdetermined nature of the inversion procedure.

Using a laboratory-measured dataset, we found that the SIM gave accurate and stable results even in the case of some tens of unknowns(individual mechanisms considered).We compared this result to those given by DLSQ(with DEM in the forward modeling)and found them consequent.The SIM was used for the inversion of both P- and S-wave velocity datasets.Using the results of Walsh(1965) and the method of David and Zimmerman (2012), we related the characteristic pressures to the aspect ratios of the defects and gave an estimate for their values.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The research was supported by the European Union,co-financed by the European Social Fund and the GINOP-2.315-2016-00010“Development of enhanced engineering methods with the aim at utilization of subterranean energy resources” project in the framework of the Széchenyi 2020 Plan, funded by the European Union, co-financed by the European Structural and Investment Funds.Particular thanks are expressed to F.Kristály and F.Mádai for the X-ray diffraction analysis, as well as J.Somogyi-Molnár and A.Kiss for their contribution to laboratory measurements.

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2021.08.015.