Rate-and-state models are widely used physical models for the relation betweenchanges in pore pressure due to fluid injection or gas extraction and the induced seismic hazard in a field. We consider the modification where the pore pressure measurements are affected by noise and provide explicit expressions for the first and second moments of the state variable. We show that when the pressure increases, there is positive correlation. In the case of decreasing pressure, both positive and negative correlation is possible. Using the delta method, approximate first and second moments of the rate variable are derived and compared to empirical moments

Baki, Z.

van Lieshout, M.N.M.

English

University of Twente Research Information

arXiv:2403.13413v1  [cs.CE]  20 Mar 2024On the moments of Cox rate-and-state modelsZ.Bakia and M.N.M. van Lieshoutb,aa Department of Applied Mathematics, University of TwenteP.O. Box 217, NL-7500 AE, Enschede, The Netherlands.b CWI, P.O. Box 94079, NL-1090 GB Amsterdam, The NetherlandsABSTRACT: Rate-and-state models are widely used physical models for the relation betweenchanges in pore pressure due to fluid injection or gas extraction and the induced seismichazard in a field. We consider the modification where the pore pressure measurements areaffected by noise and provide explicit expressions for the first and second moments of thestate variable. We show that when the pressure increases, there is positive correlation. Inthe case of decreasing pressure, both positive and negative correlation is possible. Using thedelta method, approximate first and second moments of the rate variable are derived andcompared to empirical moments.Keywords: Cox process, induced seismicity, pore pressure, rate-and-state model.2020 MSC Classification: 60G55.1 IntroductionThe study of induced earthquakes caused by extraction or injection of fluids or gases is animportant research topic. In the Netherlands, the Groningen gas field discovered in the late1950s was viewed in a mostly positive light for decades until the occurrence of earthquakesfrom the 1990s onward caused a shift in public opinion that eventually led to the closing ofthe field. Thus, for planning purposes, it is important to be able to predict seismic hazardbased on field measurements, for instance of pore pressure or, equivalently, Coulomb stress.One of the most widely used methodologies to do so is the rate-and-state model [1, 4, 7]which is now considered to be the state-of-the-art technique [5].In the rate-and-state model, the earthquake intensity λ(t), the rate at time t, is assumedto be inversely proportional to a state variable Γ(t), that is,λ(t) =1Γ(t), t ≥ 0.The state variable Γ(t) is defined by the ordinary differential equationdΓ(t) = α [dt+ Γ(t)dX(t)] ,where X(t) ≥ 0 is the pore pressure at time t and α > 0 is a parameter. Multiplyingboth sides by exp(−αX(t)) and discretising in time steps of length ∆ > 0, we obtain thecorresponding Euler difference equationΓ((k + 1)∆) = (Γ(k∆) +∆α) exp [α {X((k + 1)∆)−X(k∆)}] . (1)Additionally, a spatial location may be added to the model. For a fuller discussion, we referto [5].In this paper, we consider stochastic rate-and-state models in which noise is added tothe pore pressure X. After the formal definition of this model in Section 2, we give explicitexpressions of the first and second moments of the random state variable. We then applythe delta method to approximate the first and second moments of the random rate variable.The paper finishes with an illustration of the accuracy of the approximations for various porepressure scenarios.22 Cox rate-and-state modelsLet Ψ be a Cox process (cf. [3, Section 5.2]) on k∆, k = 0, 1, 2, . . . , with driving randommeasureΛ(k∆) =1Γ(k∆), k ∈ N0. (2)In contrast to the classic deterministic case discussed in Section 1, we assume that the stateΓ is a stochastic process defined by equation (1) for a deterministic starting state γ0 > 0and random pore pressure X of the form X(k∆) = m(k∆) + E(k∆). Here, m ≥ 0 is adeterministic mean and E is Gaussian white noise. In other words, each E(k∆) is normallydistributed with mean zero and variance σ2, and for k 6= k′, E(k∆) and E(k′∆) are mutuallyindependent.The Euler difference equation (1) can be solved explicitly. Indeed, for each k ∈ N,Γ(k∆) = exp [αX(k∆)]{α∆k−1∑i=0exp [−αX(i∆)] + γ0 exp [−αX(0)]}(3)under the convention that Γ(0) = γ0. The Cox rate-and-state model description is completeupon requiring that N(k∆), the number of induced earthquakes during the time period[k∆, (k + 1)∆), is Poisson distributed with rate parameter Λ(k∆)∆, k = 0, 1, . . . .3 Moments of the state variableSince the randomness in the driving random measure (2) of our Cox process is induced bythe state process Γ, we investigate its first and second-moment properties first.Theorem 1. Let Γ be defined by (3). Set c = eα2σ2 and define, for i, j ∈ N0,fij = exp [α {m(i∆)−m(j∆)}] .3Then, for k ∈ N,EΓ(k∆) = c(α∆k−1∑i=0fki + γ0fk0),VarΓ(k∆) = α2∆2c2(c2 − 1)k−1∑i=0f2ki + α2∆2c2(c− 1)k−1∑i=0k−1∑i 6=j=0fkifkj+ 2α∆γ0c2f2k0(c2 − 1 + (c− 1)k−1∑i=1f0i)+ γ20f2k0c2(c2 − 1)and, for 0 < k < l,Cov(Γ(k∆),Γ(l∆)) = α2∆2c2k−1∑i=0[fkifli(c− 1)− fli(1− 1c)]+ (2α∆γ0 + γ20)c2fk0fl0(c− 1)− α∆γ0c2fl0(1− 1c). (4)Proof. Recall the moment-generating function of a normal distribution to conclude thatEetZ = exp(µt+ σ2t2/2)when Z is normally distributed with mean µ and variance σ2. The formula for the meanfollows immediately. Moreover, for 0 < k ≤ l,Cov(Γ(k∆),Γ(l∆)) = α2∆2k−1∑i=0l−1∑j=0Cov(eα[X(k∆)−X(i∆)], eα[X(l∆)−X(j∆)])+ α∆γ0k−1∑i=0Cov(eα[X(k∆)−X(i∆)], eα[X(l∆)−X(0)])+ α∆γ0l−1∑j=0Cov(eα[X(k∆)−X(0)], eα[X(l∆)−X(j∆)])+ γ20Cov(eα[X(k∆)−X(0)], eα[X(l∆)−X(0)]). (5)Let us work out the terms in the expression on the right-hand side of equation (5) one by one.4First, consider the double sum over i and j. Because of the independence of the componentsof the random vector E, summands for which k, l, i and j are all different do not contribute.Therefore, for 0 < k = l, the total contribution of the double sum isα2∆2c2k−1∑i=0f2ki(c2 − 1)+ α2∆2c2k−1∑i=0fkik−1∑j 6=i;j=0fkj (c− 1) .For 0 < k < l, the double sum contributes non-zero entries for i = j (which cannot be equalto k or l) and for i 6= j = k. Their total contribution isα2∆2c2k−1∑i=0fkifli (c− 1)− α2∆2c2k−1∑i=0fli(1− 1c).Next consider the two single sums on the right-hand side of (5). If 0 < k = l, they areidentical and each one is equal toα∆γ0c2[f2k0(c2 − 1) + (c− 1)k−1∑i=1fkifk0].In the case that 0 < k < l, the sum over i has a non-zero contribution only for i = 0, and thesum over j has non-vanishing contributions for j = 0 and j = k. Adding these up, we obtainα∆γ0c2fl0[2fk0(c− 1)−(1− 1c)].Finally, the last term in the expression on the right-hand side of equation (5) readsγ20f2k0c2(c2 − 1) for 0 < k = l and γ20fk0fl0c2(c − 1) when 0 < k < l. Expression (4) nowfollows from tallying up the various contributions.When fluid is injected into a field, the pore pressure typically increases. In this case, thestate variables are positively correlated.Corollary 1. Let Γ be as in Theorem 1. Then, if the function m(k∆) is increasing in k ∈ N0,Cov(Γ(k∆),Γ(l∆)) ≥ 0 for all k, l ∈ N0.5Proof. Obviously, the claim holds for k = l. For 0 ≤ i < k < l, α2∆2fli{fkic2(c− 1)− c(c −1)} ≥ 0 if σ2 = 0 or cfki ≥ 1. The latter condition is equivalent toα2σ2 + α{m(k∆) −m(i∆)} ≥ 0and is implied by the assumption that m is increasing. By the same argument, for increasingpore pressure, α∆γ0fl0{fk0c2(c − 1) − c(c − 1)} ≥ 0 and an appeal to (4) completes theproof.When the pore pressure decreases, for example due to gas extraction, the picture is morevaried.Corollary 2. Let Γ be as in Theorem 1 and assume that the function m(k∆) is decreasingin k ∈ N0.(i) If ασ2 > m(0), then Cov(Γ(k∆),Γ(l∆)) ≥ 0 for all k, l ∈ N0.(ii) Ifασ2 < mini∈N0{m(i∆)−m((i+ 1)∆)} ,the minimal drop in pressure in between observation epochs,Cov(Γ(k∆),Γ(l∆)) − (α∆γ0 + γ20)c2(c− 1)fk0fl0 ≤ 0for all 0 ≤ k < l.Proof. For 0 ≤ i < k < l, consider α2∆2flic(c − 1){cfki − 1}. The term in between curlybrackets is negative if and only ifα {m(k∆)−m(i∆)}+ α2σ2 < 0.6Years1995–2001 179.81 177.39 174.86 172.20 169.42 166.50 163.482002–2008 160.32 157.05 153.65 150.13 146.49 142.72 138.822009–2015 134.81 130.68 126.43 122.04 117.53 112.91 108.162016–2021 103.28 98.29 93.17 87.94 82.56 77.08Table 1: Estimated pore pressure in bara on January 1st in the years 1995–2021 near thetown of Slochteren, The Netherlands.Since m is decreasing, under assumption (ii),α2σ2 < m(i∆)−m((i+ 1)∆) ≤ m(i∆)−m(k∆).We conclude that α2∆2flic(c − 1){cfki − 1} ≤ 0. Since the same argument can be used toshow that α∆γ0fl0c(c− 1){cfk0 − 1} ≤ 0, (ii) is seen to hold. Under assumption (i),α2σ2 > m(0) ≥ m(i∆) ≥ m(i∆)−m(k∆)and therefore α2∆2flic(c − 1){cfki − 1} ≥ 0. Similarly, α∆γ0fl0c(c − 1){cfk0 − 1} ≥ 0. Theobservation that Γ(0) = γ0 is deterministic completes the proof of (i).Example 1. Table 1 lists estimated pore pressure values near the town of Slochteren inthe Groningen gas field in The Netherlands for January 1st, 1995–2021 [6]. The estimatedstandard deviation is σ = 7.17.Note that the pore pressure values are decreasing due to gas extraction. Since the intensityof induced earthquakes was very low in 1995, γ0 can be considered infinite. Equation (4) thenimplies that the covariance matrix of the random vector Γ(k∆)k has positive entries only.Example 2. Next, let us suppose that – in contrast to the previous example – the initialseismicity is very high, i.e. γ0 = 0. Assume a linearly decreasing sequence of pore pressuresm(0) = 3, m(1) = 2 and m(2) = 1 and set α = 1. Then the covariance matrix of therandom vector (Γ(1),Γ(2)) is readily calculated. Indeed, VarΓ(1) = e−2(e4σ2 − e2σ2) and7VarΓ(2) = (e−2 + e−4)(e4σ2 − e2σ2) + 2e−3(e3σ2 − e2σ2). As for the off-diagonal entry,Cov(Γ(1),Γ(2)) = Cov(e−1eE(1)−E(0), e−2eE(2)−E(0) + e−1eE(2)−E(1))= e2σ2{e−3(eσ2 − 1)− e−2(1− e−σ2)}is negative for σ2 < 1 and positive for σ2 > 1.4 Approximate moments of the rate variableRecall that in the rate-and-state model defined in Section 2, the rate of induced earthquakesis inversely proportional to the state. However, due to its form (3), the moments of 1/Γ(k∆)are intractable. Thus, in this section, we use the delta method to approximate them in termsof the tractable moments of Γ(k∆) derived in the previous section.Theorem 2. Let Γ be defined by (3). Then, for k ∈ N0,E[1Γ(k∆)]≈ 1EΓ(k∆)+VarΓ(k∆)(EΓ(k∆))3.Note that the approximation of the expectation of Γ(k∆)−1 is at least as large as its‘plug-in estimator’ 1/EΓ(k∆).Proof. Apply the delta method based on the Taylor expansion1x0 + h≈ 1x0− hx20+12!2h2x30around x0 equal to the expectation of Γ(k∆). Upon taking the expectation, one obtains thatE[1Γ(k∆)]≈ 1EΓ(k∆)− E[Γ(k∆)− EΓ(k∆)(EΓ(k∆))2]+ E[(Γ(k∆)− EΓ(k∆))2(EΓ(k∆))3].The middle term on the right-hand side is zero, and the claim follows.8Theorem 3. Let Γ be defined by (3). Then, for k, l ∈ N0,E[1Γ(k∆)Γ(l∆)]≈ 1E [Γ(k∆)] E [Γ(l∆)]+Cov(Γ(k∆),Γ(l∆))(EΓ(k∆))2(EΓ(l∆))2+VarΓ(k∆)EΓ(l∆)(EΓ(k∆))3+VarΓ(l∆)EΓ(k∆)(EΓ(l∆))3.Proof. The Taylor expansion of the function (x, y) 7→ 1/(xy) around the point (EΓ(k∆),EΓ(l∆))yields the approximation1Γ(k∆)1Γ(l∆)≈ 1EΓ(k∆)1EΓ(l∆)+(Γ(k∆)− EΓ(k∆)) (Γ(l∆)− EΓ(l∆))(EΓ(k∆))2(EΓ(l∆))2− Γ(k∆)− EΓ(k∆)(EΓ(k∆))2EΓ(l∆)− Γ(l∆)− EΓ(l∆)(EΓ(l∆))2EΓ(k∆)+(Γ(k∆)− EΓ(k∆))2(EΓ(k∆))3EΓ(l∆)+(Γ(l∆)− EΓ(l∆))2EΓ(k∆)(EΓ(l∆))3up to second-order moments. Taking expectations proves the claim.Combining Theorems 2 and 3, we obtain the following corollary.Corollary 3. Let Γ be defined by (3). Then, for k, l ∈ N0,Cov(1Γ(k∆),1Γ(l∆))≈ Cov(Γ(k∆),Γ(l∆))(EΓ(k∆))2(EΓ(l∆))2.5 Simulation examplesTo investigate the accuracy of the approximations in Section 4, we compare the formulaein Theorem 2 and Corollary 3 with population estimates. Recall that for an i.i.d. sampleX1, . . . ,Xn of some random variable X1 with mean µ and variance σ2, approximate confidenceintervals for µ and σ2 take the form(X̄n −Sn√nξ1−α/2, X̄n +Sn√nξ1−α/2)90 5 10 15 204.54.64.74.84.95.00 5 10 15 200.0060.0080.0100.0120.014Figure 1: 95% pointwise confidence intervals for mean (left) and variance (right) of Γ(k∆)−1as a function of k when m(∆k) = 6−1/(0.5k + 1), α = 0.01, ∆ = 0.1, γ0 = 0.2 and σ2 = 2.0.The dots correspond to the approximation values of Theorem 2 and Corollary 3.0 5 10 15 205.05.15.25.30 5 10 15 200.0060.0080.0100.0120.014Figure 2: 95% pointwise confidence intervals for the mean (left) and variance (right) ofΓ(k∆)−1 as a function of k when m(k∆) = 5+ 10/(2k + 1), α = 0.01, ∆ = 0.1, γ0 = 0.2 andσ2 = 2.0. The dots correspond to the approximation values of Theorem 2 and Corollary 3.10for µ, and (S2n1 + ξ1−α/2√( 2ζn−1),S2n1− ξ1−α/2√( 2ζn−1))for σ2. Here ξ1−α/2 is the (1−α/2)-quantile of the standard normal distribution. Furthermore,ζ = 12(γ4 − 1) for γ4 the ratio of the fourth central moment and the squared variance. Sinceγ4 is unknown, we estimate it by γ̂4 = n∑(Xi −m)4/((n − 1)2S4n) where m is the samplemedian [2]. In our case, Xi = 1/Γi(k∆). A similar approach may be taken for the covariance,although we do not pursue it here.We consider two cases: increasing and decreasing pore pressure. For Figure 1, define thepore pressure by the increasing functionm(∆k) = 6− 11 + k/2.With parameter values α = 0.01, ∆ = 0.1, γ0 = 0.2 and σ2 = 2.0, the 95% pointwiseconfidence intervals for the mean and variance are given in the left- and right-most panels.In both cases, the sample size was n = 500. It can be seen that the approximations are quiteadequate.For Figure 2, take the decreasing functionm(∆k) = 5 +102k + 1.With the same parameter values and sample size, the 95% pointwise confidence intervals forthe mean and variance are given in the left- and right-most panels. Again the approximationsare satisfactory.6 ConclusionIn this paper, we modified the rate-and-state model proposed in the earth sciences to describethe relation between changes in pore pressure and induced seismic hazard [1, 4, 7] in such a11way that measurement error in the pore pressure can be taken into account. For the timeseries of state variables Γ in the resulting Cox process [3], closed-form expressions for the firstand second moments were derived. Additionally, we provided sufficient conditions for positiveand negative correlation. On top of that, we gave approximate expressions for the first twomoments of the rate variable and verified in simulated examples that the approximations aresatisfactory.AcknowledgementsThis research was funded by the Dutch Research Council NWO through their DEEP.NLprogramme (grant number DEEP.NL.2018.033).References[1] Candela, T. et al. (2019). Depletion-induced seismicity at the Groningen gas field:Coulomb rate-and-state models including differential compaction effect. Journal of Geo-physical Research: Solid Earth 124:7081–7104.[2] Curto, J.D. (2023). Confidence intervals for means and variances of non-normal distri-butions. Communications in Statistics – Simulation and Computation 52:4414–4430.[3] Chiu, S.N., Stoyan, D., Kendall, W.S., Mecke, J. (2013). Stochastic Geometry and itsApplications. Wiley, 3rd ed.[4] Dempsey, D., Suckale, J. (2017). Physics-based forecasting of induced seismicity atGroningen gas field, the Netherlands. Geophysical Research Letters 22:7773–7782.[5] Kűhn, D., Hainzl, S., Dahm, T., Richter, G., Vera Rodriguez, I. (2022). A review ofsource models to further the understanding of the seismicity of the Groningen field.Netherlands Journal of Geosciences, 101:e11.12[6] Lieshout, M.N.M. van, Baki, Z. (2023). Exploring seismic hazard in the Groningen gasfield using adaptive kernel smoothing. Mathematical Geosciences, published online.[7] Richter, G., Hainzl, S., Dahm, T., Zőller, G. (2020). Stress-based statistical modelingof the induced seismicity at the Groningen gas field, The Netherlands. EnvironmentalEarth Sciences 79, 252.13

On the moments of Cox rate-and-state models

https://ris.utwente.nl/ws/files/359411676/2024-XBakiLieshout_2403.13413.pdf

On the moments of Cox rate-and-state models

Abstract

Similar works

Full text

Available Versions

University of Twente Research Information