# Memory Effects on Nonlinear Temperature and Pressure Wave Propagation in the Boundary between Two Fluid-Saturated Porous Rocks

Authors: R.  Garra; E. Salusti; R. Droghei
Journal: Advances in Mathematical Physics
 (2015)

---

## Abstract

The evolution of strong transients of temperature and pressure in two adjacent fluid-saturated porous rocks is described by a Burgers equation in an early model of Natale and Salusti (1996). We here consider the effect of a realistic intermediate region between the two media and infer how transient processes can also happen, such as chemical reactions, diffusion of fine particles, and filter cake formations. This suggests enlarging our analysis and taking into account not only punctual quantities but also “time averaged” quantities. These boundary effects are here analyzed by using a “memory formalism”; that is, we replace the ordinary punctual time-derivatives with Caputo fractional time-derivatives. We therefore obtain a nonlinear fractional model, whose explicit solution is shown, and finally discuss its geological importance.


---

## Body

## 1. Introduction

Modeling of nonlinear wave transition in fluid-saturated porous rocks is very important for a large number of applications. A realistic model of these phenomena is however rather hard to obtain. As an example, if we consider a transient wave crossing a two-layer medium, the boundary between the two rocks can have a particular importance but it is often schematized as just a jump in the initial/boundary conditions. It is important to notice that a quickly moving wave, coming from a hot pressurized rock and entering into a similar colder medium, has to cross their realistic boundary layer which certainly is not an abrupt mathematical plane. It is indeed more realistic to describe it as a thin region of transition, characterized by continuous evolution. All this is schematized here both in general grounds and as the model of a detailed boundary that describes a continuous linear transition.The sudden arrival of a transient wave can moreover cause perturbations, as those due to fine particle migration, chemical reaction, filter cake formation, and so forth [1–3]. In order to gain a more realistic knowledge of these time dependent phenomena, one can focus the attention onto time averages of pressure P and temperature T. In such context, it is of interest that Caputo [4] discussed generalized Darcy’s law for fluid fluxes in porous media, by using memory formalism. It mathematically rests on replacing the ordinary time-derivatives with fractional derivatives, which correspond to a time average of pore pressure P and fluid-rock temperature T. Notice how recently these fractional models were used to describe transient waves in inhomogeneous porous media [5] while the issue of fractional thermoelasticity has been the object of the recent monography by Povstenko [6].In such context, we here generalize a model proposed by Rice and Cleary [7], McTigue [8], and Bonafede [9], among many others, of thermomechanical transients of P and T propagating in two adjacent layers of fluid-saturated porous rocks, in a simple 1D case. This model has found relevant geological applications in the study of bradyseismic crisis at the Campi Flegrei [9] and at the Izu Peninsula [10]. In addition, Natale and Salusti [11] studied the full nonlinear properties of such model, showing how it can be reduced to a single relation, the classical Burgers equation. Its solutions are strong quick transients of P and T. These early models consider an initial sharp jump of pressure and temperature between the two layers, at the space coordinate z
=
0. More recently, Garra [12] found an analytic solution of the fractional generalization of this model in the inviscid case.The aim of this paper is to study the effect of a smooth transition in the boundary layer between different media, in the propagation of nonlinear thermoelastic waves. In particular, we discuss the physical meaning related to the novel boundary layer condition, as discussed in Section2 of the present paper. This finally gives additive terms in the model equation for temperature and pressure fields. We then also consider the role played by memory effect in the wave propagation by means of a time-fractional approach. Fractional calculus indeed provides useful mathematical tools to deal with physical systems where memory effects can not be neglected and a power law behavior of the memory kernel is present. This is the reason why the time-fractional derivative is by now a popular tool in mathematical models for wave propagation in porous media and in many other different fields of applied sciences. On the other hand, nonlinear fractional models, even if physically motivated, are still poorly studied, because few exact results are known and the main part of the literature is concentrated on numerical and semianalytical methods (see, e.g., [13] and references therein). In this paper, we are able to present a new exact solution of the nonlinear fractional model obtained from physical considerations.In this study, we first recall the early models and the equation describing the effect of the boundary between two rocks (Section2). Then, we study the fractional generalization of the Natale-Salusti model, taking into account the complex structure of the thin region at the boundary of the two rocks, that is, assuming realistic initial profiles of P and T to simulate the transition zone between the two rocks (Section 3). Some realistic fractional models are moreover described in Section 4. A discussion of the obtained results is in Section 5.


## 2. Thermoelastic Response to Stresses of Fluid-Saturated Porous Rocks: The Early Models

We here analyze a homogeneous, isotropic, deep horizon, identified as a “source” of temperature or pressure change, covered by a fluid-saturated porous-permeable medium with lower temperature and pressure. These two media are first isolated but at a given moment, say at timet
=
0, they interact and give origin to transient phenomena. Their nature can be diffusive or strongly nonlinear where advection is considered, depending on the geophysical properties of the rocks. The control parameter is a Reynolds number as we are going to discuss. In many cases, this sharp transition between two media with strongly varying temperature and pressure can be modelized by means of a step function.To describe such interaction, first we recall [8] the fluid mass balance (1)
∂
ς
∂
t
+
∇
·
ϕ
U
=
0
,where ς is the fluid-rock content, ϕ the porosity, and the fluid velocity vector is U. Darcy’s law gives the relation between U and P inside the porous medium:(2)
U
=
-
K
f
ϕ
μ
∇
P
.Here K
f is the medium permeability and μ is the fluid viscosity. Taking also into account the linear stress and strain relation, McTigue [8] shows that in a thermoporoelastic medium pressure P is ruled by the following evolution equation, depending on the temperature field:(3)
∂
∂
t
-
k
∗
∇
2
t
r
σ
i
j
+
3
P
B
=
4
G
k
∗
1
+
ν
∗
α
s
1
-
ν
∗
∇
2
T
+
2
G
B
1
+
ν
∗
1
+
ν
u
∗
3
ν
u
∗
-
ν
∗
ϕ
α
f
-
α
S
∂
∂
t
T
=
3
α
∗
B
∇
2
T
+
3
α
B
∂
∂
t
T
,where k
∗
=
(
2
K
f
G
/
μ
)
B
2
(
1
-
ν
∗
)
(
1
+
ν
u
∗
)
2
/
9
(
1
-
ν
u
∗
)
(
ν
u
∗
-
ν
∗
) plays the role of the hydraulic diffusion coefficient, B is Skempton parameter, G is the shear modulus, σ
i
j is the stress tensor, ν
∗ is the drained Poisson ratio and ν
u
∗ is the undrained Poisson ratio, and α
s is the volumetric thermal expansion coefficient for the solid (fluid). In (3), α and α
∗ are moreover implicitly defined as(4)
α
∗
=
4
G
B
k
∗
1
+
ν
∗
α
S
3
1
-
ν
∗
,
α
=
2
G
B
2
1
+
ν
∗
1
+
ν
u
∗
9
ν
u
∗
-
ν
∗
ϕ
α
f
-
α
S
.Equation (3) plays a fundamental role in such problems since it strictly relates of P and T. It shows that the evolution of P is due to the gradients of T, and the opposite as well, namely, the resulting transients, are strictly interconnected as transients of P and T.McTigue [8] proved also that, for a simple 1D geometry with space coordinate z, the stress is constant and therefore (3) reduces to(5)
∂
∂
t
-
k
∗
∂
2
∂
z
2
P
-
α
∗
∂
2
∂
z
2
T
-
α
∂
∂
t
T
=
0
.As a second equation, Bonafede [9] shows that the energy balance for the whole system, assuming a local thermal equilibrium between the fluid and the solid, is(6)
∂
T
∂
t
-
K
T
φ
ρ
f
c
f
+
1
-
φ
ρ
m
c
m
∇
2
T
=
-
ρ
f
c
f
U
·
∇
T
φ
ρ
f
c
f
+
1
-
φ
ρ
m
c
m
,where K
T is the average thermal conductivity, ρ
m is the matrix density, c
f is an empirical fluid heat capacity, and c
m is that of the matrix. Equation (6) is the classical equation of heat transfer ([14], but see also [15]) where convection is also considered. Notice how these empirical quantities are often difficult to estimate but McTigue [8] was able to obtain a similar relation theoretically.In a simple 1D model, the governing equations (2), (5), and (6) can be arranged in(7a)
∂
T
∂
t
-
k
∂
2
T
∂
z
2
+
β
∂
T
∂
z
∂
P
∂
z
=
0
,
(7b)
∂
∂
t
-
k
∗
∂
2
∂
z
2
P
-
α
∗
∂
2
∂
z
2
T
-
α
∂
∂
t
T
=
0
,where k
=
K
T
/
ϕ
ρ
f
c
f
+
(
1
-
ϕ
)
ρ
m
c
m is the average thermal diffusivity due to diffusion and the average thermal diffusivity due to convection is β
=
-
K
f
ρ
f
c
f
/
μ
ϕ
[
ϕ
ρ
f
c
f
+
(
1
-
ϕ
)
ρ
m
c
m
].To describe the above-mentioned flows from a first rock into the other matrix, the boundary conditions of the previous works of Natale and Salusti [11] for P and T are a sharp jump at z
=
0, namely,(8a)
at
z
⟶
∞
P
=
P
0
+
P
1
T
=
T
0
+
T
1
at
z
=
0
P
=
P
0
T
=
T
0
,and the initial conditions at t
=
0 are(8b)
at
z
<
0
P
=
P
0
T
=
T
0
at
z
>
0
P
=
P
0
+
P
1
T
=
T
0
+
T
1
,where P
0, P
1, T
0, and T
1 are constant. This means that such transient waves are generated around t
~
0 at z
~
0 and T and P are constants for large z.The above equations depend on many empirical coefficients and we have neglected the possible temperature or pressure variability of these geological parameters. Moreover, the knowledge of all these coefficients is a relevant research issue in order to know the cases where nonlinear phenomena play an important role. We refer, for example, to Bonafede [9] and Natale and Salusti [11] for an estimate of these parameters for clay rocks, Berea sandstones, and Ruhr sandstones.We observe that these are solutions for an infinite domain: rather regular patches of rocks represent a more realistic idealization, assuming space-scales larger than the evolution region but certainly with a finite extension.In order to treat analytically the model equations (7a)-(7b) can be simplified as follows (see, e.g., [9]):(9a)
∂
P
∂
t
-
k
∗
∂
2
P
∂
z
2
=
α
∂
T
∂
t
,where the fluid diffusivity k
∗ and the thermal expansion coefficient α are defined in (3). In the early treatment of NS, it is moreover assumed that k
∗
~
0, that is, the diffusive term, is neglected. We refer to Merlani et al. [3] for a discussion about characteristic rocks where diffusive terms can be neglected.In some cases, one can obtain a less stringent result. Indeed we can introduce here the assumption of Merlani et al. [3] related to a mathematical analysis of the group invariants of these equations, which allow solutions of rigid wave translation f
z
-
v
t or of autosimilarity f
z
/
t. This suggests investigating solutions such that(9b)
k
∗
∂
2
P
∂
z
2
+
α
∗
∂
2
T
∂
z
2
=
d
d
t
f
tis anansatz that clearly must be checked a posteriori. Both approaches (9a) and (9b) imply that (10)
P
z
,
t
-
α
T
z
,
t
≈
g
z
+
f
t
,where g
(
z
)
=
P
(
z
,
0
)
-
α
T
(
z
,
0
) is determined by the initial conditions of both T and P. Notice that (2) and (10) also give for the Darcy velocity that (11)
U
z
,
t
=
-
K
f
ϕ
μ
d
d
z
α
T
z
,
t
+
g
z
.IfP
(
z
,
0
)
-
α
T
(
z
,
0
) are constants as in (8a) and (8b), we have g
=
c
o
n
s
t. Inserting (10) in (7a), Natale and Salusti [11] obtain a classical nonlinear Burgers-type equation: (12)
∂
T
∂
t
-
k
∂
2
T
∂
z
2
+
α
β
∂
T
∂
z
2
=
0that can be linearized by the Hopf-Cole transformation [16]. Its solution is ruled by the fundamental Reynolds number (13)
R
=
α
β
kand in its general form (up to an additive constant term) is(14)
T
z
,
t
=
k
t
e
R
-
1
exp
⁡
-
z
2
/
4
k
t
π
+
e
R
-
1
∫
z
/
4
k
t
∞
exp
⁡
-
ζ
2
d
ζbut for a small k, that is, R
>
10, one has a simple solution and front for z
=
z
B [16], namely,(15)
T
z
,
t
≈
T
0
+
z
2
4
α
β
t
,
0
<
z
<
z
B
=
4
α
β
T
1
t
.It is important to stress that, for a large R (see (13)), the front z
B
≈
4
α
β
T
1
t is characterized by T
(
z
B
)
=
T
I. The corresponding Darcy’s velocity is(16)
U
=
-
K
f
ϕ
μ
∂
P
∂
z
≅
-
z
2
β
t
.On general grounds, since in (12) there are only derivatives of T and the space derivatives are of even order, if T
(
z
) is a solution also T
(
-
z
) and T
(
z
)
+
c
o
n
s
t are still solutions. The roles of T
1 and P
1 are therefore essential, and different dynamics are strongly related to their sign: in this work, we focus on T
I
>
0 and P
I
>
0.To conclude this section, we observe that, assuming more realistic continuous transition between the two rocks (parameterized by means of the functiong
(
z
) that schematizes a steady large scale trend of P and/or T) and using (10) in (7a), we arrive to the more general model equation: (17)
∂
T
∂
t
-
k
∂
2
T
∂
z
2
+
α
β
∂
T
∂
z
2
+
β
∂
g
∂
z
∂
T
∂
z
=
0
.


## 3. Fractional Thermoelasticity: A Short Survey

We have seen that in these problems one has to deal with a great number of empiric geological parameters. Moreover, some of them, as permeability, thermal diffusivity, and so on are not constants but can be altered by the transient flow and become time-dependent. In the framework of the theory of dynamical coefficients in fluid-saturated porous media [17], on the basis of experimental studies on sound waves propagation, Fellah et al. [5] have shown that some coefficients are frequency-dependent. In addition, many analyses [2, 3] show that transient flows can bring fine and large particles from one rock to the adjacent media and this can affect their permeability, porosity, and so on which finally give a time delay in the front evolution. So, to take into account these time-dependences of the coefficients, it is useful to introduce a time averaged approach by means of fractional operators [4]. Fractional models in linear viscoelasticity are widely used, starting from the works by Caputo and Mainardi [18] and Bagley and Torvik [19] that discuss the general theory and implications of the fractional approach in viscoelasticity problems (see the recent book by Mainardi [20]).Following a similar heuristic viewpoint, we now introduce the fractional generalization of the early Natale-Salusti model. We recall that there are many attempts in literature to treat the governing equations of thermoelasticity with an integrodifferential approach ([21], and references therein), but few exact results about nonlinear models are present.The definition of the Caputo fractional derivative of orderν
∈
(
0,1
) is [22](18)
∂
t
ν
a
h
t
=
1
Γ
1
-
ν
∫
a
t
∂
t
h
τ
t
-
τ
ν
d
τ
=
h
′
τ
∗
t
-
τ
-
ν
Γ
1
-
ν
,where this ∂
t
ν
a is an integrodifferential operator that satisfies the following properties for ν
∈
(
0,1
) and β
>
-
1
/
0:(19)
∂
t
ν
a
J
ν
h
t
=
h
t
∂
t
ν
a
h
t
=
J
1
-
ν
d
d
t
h
t
∂
t
ν
a
t
β
=
Γ
β
+
1
Γ
β
-
ν
+
1
t
β
-
ν
J
ν
∂
t
ν
a
h
t
=
h
t
-
h
0
,where(20)
J
ν
h
t
=
1
Γ
ν
∫
0
t
h
τ
t
-
τ
ν
-
1
d
τis the Riemann-Liouville fractional integral. This memory formalism can also be seen as a time averaging method. Following these ideas, one can consider a time-fractional version of (5), namely,(21a)
∂
ν
∂
t
ν
-
k
f
∂
2
∂
z
2
P
-
α
∗
∂
2
T
∂
z
2
-
α
∂
ν
∂
t
ν
T
=
0
,and the nonlinear energy balance equation, under the assumption (10)(21b)
∂
ν
T
∂
t
ν
-
k
∂
2
T
∂
z
2
+
α
β
∂
T
∂
z
2
+
β
∂
g
∂
z
∂
T
∂
z
=
0
.In (21a) and (21b), we consider a unique time-scale evolution for both P and T and therefore the real order ν is the same in (21a) and (21b), with 1
/
2
<
ν
<
1 for technical reasons that will be explained in the next section. This physical assumption of our model is due to the fact that the time-evolution of pressure and temperature fields is strictly related by means of the constitutive equations (3) and (5) discussed in the previous section. On the other hand, from the mathematical point of view, it is not difficult to generalize these techniques in order to obtain more general results also if the fractional order appearing in (21a) and (21b) is different.


## 4. Fractional Generalizations of the Natale and Salusti Model

We have shown a fractional generalization of the model equations (7b) and (17), to deal with time averaged quantities. We consider in particular the following initial condition on the pressure field:(22)
P
z
,
0
=
g
z
=
P
I
1
-
z
ψ
,
d
g
d
z
=
-
P
I
ψ
,where P
I
>
0 and ψ is the thickness of the boundary layer between the two rocks, while the initial condition on the temperature field T
(
z
,
0
) is taken as in (8b). Since in this boundary interval only P
(
z
,
0
) varies, equation (21b) becomes(23)
∂
ν
T
∂
t
ν
-
k
∂
2
T
∂
z
2
+
α
β
∂
T
∂
z
2
+
β
∂
g
∂
z
∂
T
∂
z
=
∂
ν
T
∂
t
ν
-
k
∂
2
T
∂
z
2
+
α
β
∂
T
∂
z
2
-
β
∂
T
∂
z
P
I
ψ
=
0
.This equation schematizes the time averaged effect of the border between the two rocks for pressure wave arrival. In such thin layer, much more complex phenomena can happen, as fine particles can leave the source rocks and filtrate towards the downstream rock until the eventual formation of filter cakes, chemical reactions, or also local rock fracturing [3]. As suggested by the previous cases, we tentatively assume that (24)
T
z
,
t
=
m
t
z
2
+
n
t
z
+
s
t
.Thisansatz can be rigorously justified by means of the invariant subspace method introduced by Galaktionov [23] to find explicit solutions of some cases of nonlinear partial differential equations.Inserting suchT in (23), we obtain(25)
d
ν
m
d
t
ν
+
4
α
β
m
2
=
0
,
d
ν
n
d
t
ν
+
4
α
β
n
-
2
β
P
1
ψ
m
=
0
,
d
ν
s
d
t
ν
-
2
k
m
+
α
β
n
2
-
β
P
1
ψ
n
=
0
.The first is a time-fractional Riccati-type equation whose solution is(26)
m
t
=
-
1
4
α
β
Γ
1
-
ν
Γ
1
-
2
ν
t
ν
.We stress how, because of the Gamma function at the numerator of (26), these solutions have different signs, negative in the case ν
<
1
/
2 and positive in the case ν
>
1
/
2 while, for ν
=
1
/
2, (26) diverges. Substituting (26) in the second equation of (25), we find that (27)
n
t
=
P
1
2
ψ
α
+
n
0
t
-
ν
.Finally, we obtain the function s
(
t
) by substituting (26) and (27) in the third equation of (25) and by fractional integration:(28)
s
t
=
-
k
2
α
β
Γ
1
-
ν
Γ
1
-
2
ν
f
t
-
α
β
Γ
1
-
2
ν
Γ
1
-
ν
n
0
2
t
ν
+
β
P
1
2
t
ν
4
α
ψ
2
Γ
1
+
ν
,where we denoted with f
(
t
) a locally integrable function such that t
-
ν is its Caputo fractional derivative of order ν. In this way, we finally obtain (29)
T
z
,
t
=
T
I
-
1
4
α
β
Γ
1
-
ν
z
2
Γ
1
-
2
ν
t
ν
+
P
1
2
ψ
α
+
n
0
t
ν
z
+
s
twhich for ν
>
1
/
2 gives a solution somehow similar to the previous solution (15). This is the reason why we have considered the constraint 1
/
2
<
ν
<
1 in order to obtain a physical meaningful generalization of the solution of the original Natale-Salusti model.In more detail, such novel solution has a structure similar to the punctual solutions (15) but with extra terms linear in z.One has to remember that the group theoretical analysis of (7a)-(7b) due to Merlani et al. [10] suggests also the importance of solutions as f
(
z
-
V
t
), namely, rigidly moving waves that however have different initial/boundary conditions and less practical interest.


## 5. Discussion

The detailed analysis of the thin boundary layer between two different fluid-saturated rocks and its effects on transients of fluid pressure or/and temperature waves is an important but rather complex problem. We here focus attention on one of such effects. We assume that a homogeneous “source” matrix is warmer and more pressurized than another adjacent matrix and thus we investigate the effects in the thin barrier separating these two systems. Such barrier is seen as a boundary region with smooth continuous trends of pressure and/or temperature between the two matrices. Then, at a timet
=
0 the barrier is somehow broken and consequently a flux of heat and pressure start moving from the “source,” first in the boundary layer and then in the adjacent rock.It is evident that all this may happen in many different ways, as a pure pressure flux or pure heat flux or a mixture of the two, that the “source” can have higher temperature and pressure, temperature only, pressure only, and so forth. We find a general equation describing these cases. Moreover, we have discussed in detail one of these cases, chosen because of its interest and mathematical simplicity, that is, a coupled temperature and pressure transient generated by a pure pressure transition from the “source” rock and then flowing in the rock boundary and its final evolution as aT
/
P wave in the adjacent rock. We see that all this gives rather unexpected acceleration of such transient. On the other hand, we show that a thermal increase in the boundary is at the origin of a transient delay, which is most often seen in experiments.All this concerns punctual characteristics of the system, but introducing time averaged quantities, by using fractional time-derivatives, one has a different situation. If a linear pressure variation in the boundary is considered, then we have a situation that reminds us of the previous punctual cases. If in the boundary region the initial pressure is however sharper, mathematical difficulties increase and only a different model solution is found, a kind of rigid wave. All this will be discussed in a more mathematical viewpoint in a work in progress.We finally outline some interesting issues related to this research paper for possible future work.A first relevant problem for applied geophysical problems is the estimate of the fractional order of the time-derivative appearing in the evolution equations governing the physical phenomena by inverse problem. Different studies about fractional inverse problems have recently appeared in literature; we refer, for example, to the review of Jin and Rundell [24].Another interesting issue regards the generalization of the considered model in different more realistic directions: multidimensional cases and different orders of fractionality in the governing equations. A further relevant effect that cannot be neglected in many cases is the dependence of porosity by the temperature or pressure fields that introduce higher order nonlinearities in the field equations. Therefore, our research is a first step in this rich field of research.



---

*Source: 532150-2015-11-05.xml*