C J Van

Duijn Darcy

Andro Mikelić
email: mikelic@univ-lyon1.fr

Mary F Wheeler

Keywords: Two-scale expansions, Porous media, Linear thermoelasticity, Fluid-structure interaction, Heat transfer, Upscaling

In this paper we derive a macroscopic model for thermoporoelasticity from the pore scale linearized fluid-structure and energy equations. We consider the continuum mechanics thermodynamically compatible pore scale equations corresponding to realistic rock mechanics parameters. They are upscaled using two-scale asymptotic expansions. For the upscaled equations a Lyapunov functional (a generalization of Biot's free energy) is constructed and the well-posedness

Introduction

There are serious energy and environmental problems that require advanced mathematical and high fidelity numerical modeling of thermoporoelasticity processes in geomaterials. These include geological disposal of heat-emitting nuclear fuel (nuclear waste management), geological sequestration of green house gases in supercritical form, extractions of geothermal energy, mechanics of frozen soils, and enhanced oil recovery using steam injection and hot water flooding. Scientific issues that need to be addressed involve the mechanics of porous skeletons due to elastic, plastic and creep behaviour, fluid transport characteristic of the intact skeleton and alternations in the skeleton due to mechanical, fluid flow and heat transfer processes. Additional complexities arise due to anisotropy, heterogeneity of the porous media, geochemistry effects, and multiphase fluids. Efforts in the engineering design and management of these thermoporoelastic examples are difficult and costly. Predictive computational simulation may be the only means to account for the lack of complete characterization of the subsurface environment, the multiple scales of the various interacting processes, the large areal extent of saline aquifers in the case of carbon sequestration , and the need for long time predictions. Moreover, since constitutive functionals and parameters must be determined experimentally there is a need for coupling modeling at both the pore laboratory and field scales.

To further motivate the importance of thermoporoelastic studies, we refer to Rummelhoff [30] who stated that without carbon capture and storage, it is not realistic to meet the global climate target as defined in the Paris Agreement. A massive scale up of the number of carbon capture and storage (CCS) projects is needed, and collaboration and sharing of knowledge are essential to accelerating the development.

There are numerous isothermal CCS studies, including those carried out one of the authors of this paper. However, very limited modeling that includes thermal effects has been done to date. Our objective in this first paper is to focus on fundamental modeling aspects of thermoporoelasticity, using homogenization and formal two-scale expansions. Later work will focus on rigorous mathematical results, benchmark problems and on numerical simulation of nuclear waste and CCS examples.

The equations describing the coupling between mechanical effects and flow in porous media were first introduced in the pioneering papers of Biot: see for instance the collection of Biot's papers in Tolstoy [START_REF] Tolstoy | Acoustics, elasticity, and thermodynamics of porous media. Twenty-one papers[END_REF]. In these works one roughly observes two regimes: consolidation processes, described by the quasi-static form of the equations, and wave propagation, described by their dynamic form.

To understand the origin of Biot's equations one distinguishes two approaches. One is thermodynamically based as explained in the work of Coussy et al ([11], [START_REF] Coussy | From mixture theory to Biot's approach for porous media[END_REF], [13]) or Lewis and Schrefler [START_REF] Lewis | The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media[END_REF].

The other approach -the so-called homogenization approach -considers Biot's equation as the effective (upscaled) form of the pore level (microscopic) equations.

Starting point in the homogenization approach are the continuum mechanics equations at the pore scale. Using the dimensionless pore size as the small parameter ε, the effective (upscaled) equations are found in the form of a power series in ε, with the coefficients being functions of the macroscopic variable x and the microscopic variable y = x/ε. Well-known references using such two-scale expansions are Auriault and Sanchez-Palencia [START_REF] Auriault | Etude du comportement macroscopique d'un milieu poreux saturé déformable[END_REF] and Burridge and Keller [START_REF] Burridge | Poroelasticity equations derived from microstructure[END_REF]. For a complete overview from an engineering perspective we refer to the monograph by Mei and Vernescu [START_REF] Mei | Homogenization methods for multiscale mechanics[END_REF].

Rigorous mathematical results concerning the dynamic case are found in the monograph by Sanchez-Palencia [START_REF] Sanchez-Palencia | Non-Homogeneous Media and Vibration Theory[END_REF]. Here the convergence of the Laplace transform in time of the displacements in the homogenization limit ε → 0 was established. The complexity of such multi-scale problems motivated Nguetseng ([26], [START_REF] Nguetseng | Asymptotic analysis for a stiff variational problem arising in mechanics[END_REF]) and Allaire [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] to introduce two-scale convergence techniques.

Rigorous results for the quasi-static regime, including interface conditions with an elastic body are given by Mikelić and Wheeler [START_REF] Mikelić | On the interface law between a deformable porous medium containing a viscous fluid and an elastic body[END_REF].

In this paper we are concerned with the problem of quasi-static thermal poroelasticity: coupling mechanics, flow and temperature in the regime where the acceleration terms are small.

The corresponding macroscopic equations being used in the engineering literature are often based on direct thermodynamic modeling. This approach extends Biot's pioneering work on the thermodynamics of deformable porous media (see for instance the collection [START_REF] Tolstoy | Acoustics, elasticity, and thermodynamics of porous media. Twenty-one papers[END_REF]) and is described in the work of Coussy et al ([11], [START_REF] Coussy | From mixture theory to Biot's approach for porous media[END_REF], [13]). It yields a global macroscopic model based on mixture theory, where the equations are posed separately for the fluid and solid phases. The theory also deals with nonlinear porous continua: a recapitulation of the governing equations is given in ( [13], p. 58). Macroscopic constitutive equations for thermoelastic processes in a fluid-saturated porous medium, with small heat advection and small thermoelastic coupling, are presented in [START_REF] Suvorov | Macroscopic constitutive equations of thermo-poroelasticity derived using eigenstraineigenstress approaches[END_REF] .

A formal homogenization approach in the thermoelastic case is considered by Lee and Mei ([19], [START_REF] Lee | Thermal consolidation in porous media by homogenization theory-II. Calculation of effective coefficients[END_REF]). Using the equations at the pore level they consider the perturbation of an equilibrium state. This state is characterize by a hydrostatic equilibrium for the fluid and a static equilibrium for the initial temperature. The perturbation satisfies the Boussinesq approximation for the viscous thermal flow in the pores and the linear elasticity equations coupled to heat conduction in the solid skeleton. The resulting system is written in the Eulerian description. In this description the fluidsolid interface moves in time. Small displacements of order ε , being the dimensional pore size, are assumed. This allows for an expansion of the interface conditions at the original static position. The imposed limitation on the microscopic displacement leads to a small upscaled displacement.

Comparing the homogenized constitutive laws in Lee and Mei to the results of Coussy et al, one observes that the poroelastic coefficients coincide but that the coefficients corresponding to the thermal effects deviate.

The goal of this paper is to revisit the derivation of the homogenized model, dropping the smallness conditions on the displacements and pressure (see equation [START_REF] Mctigue | Thermoelastic Response of Fluid-Saturated Porous Rock[END_REF] on page 129 of [START_REF] Lee | Thermal consolidation in porous media by homogenization theory-I. Derivation of macroscale equations[END_REF]). Moreover, we take into the account the effect of heat generation by deformation.

The pore scale equations are based on the assumptions:

i) The quasi-static deformation and heat flow in the elastic skeleton are described by the Navier linear elasticity system and the heat equation, while the stress-temperature tensor contribution is kept in the entropy. The equations are considered in the Lagrangian description.

ii) The fluid in the pores satisfies the Boussinesq approximation and inertia is disregarded.

Our aim is to derive an effective thermoporoelastic model in which the upscaled displacement can take values up to 1 -5% of the domain size and in which filtration rates have realistic values.

In addition we keep the thermoelastic coupling term in the energy balance, which allows us to investigate its consequences. For instance, its presence is crucial in a derivation of a Lyapunov functional of the upscaled system. We circumvent the difficulty caused by the moving fluid-solid interface by using the Lagrangian description for the equations of the elastic skeleton and the ALE (Arbitrary Lagrangian Eulerian) ( [START_REF] Donea | An arbitrary Lagrangian-Eulerian finite element method for transient dynamic fluid-structure interactions[END_REF], [START_REF] Hughes | Lagrangian-Eulerian finite element formulation for incompressible viscous flows[END_REF]) description for the fluid equations. The ALE description is such that the initial configuration of the fluid-solid interface remains unchanged for all times (e.g., [START_REF] Richter | Finite elements for fluid-structure interaction in ALE and fully Eulerian coordinates[END_REF], [START_REF] Wick | Solving monolithic fluid-structure interaction problems in arbitrary Lagrangian Eulerian coordinates with the deal. II library[END_REF]).

The pore scale equations are introduced in Section 2 in detail and realistic data are used to write them in a dimensionless form. In Section 3 we present the upscaling of the dimensionless micro-scale model, as introduced in Section 2. The standard steps in the derivation of Biot's equations by the two-scale expansion are alternating with the expansion of the temperature equations. The summary of the upscaled equations is presented in Section 4. They take a form similar to the equations in ( [13], page 75) and in [START_REF] Lee | Thermal consolidation in porous media by homogenization theory-I. Derivation of macroscale equations[END_REF]. However there are noticeable and non-trivial differences. In Coussy et al the equations are presented in a very general form, including large deviations. therefore it is difficult to compare terms in the equations and to quantify the (matricial) coefficients. In our approach one is able to compute all coefficients and matrices in terms of the microstructure. In Lee and Mei the term describing heat generation by the solid deformation is missing and some smallness assumptions are imposed. As a consequence, a quantitative comparison is virtually impossible.

In Section 4 we also give a Lyapunov functional for the upscaled equations. Its form is a generalization of the Biot free energy from the isothermal case. It has similarities with the free energy from [START_REF] Coussy | From mixture theory to Biot's approach for porous media[END_REF], but it does not coincide. This stability property can be used in numerical algorithms when solving the equations for large time intervals. This is crucial in simulating problems of practical interest.

Section 5 is dedicated to the calculation and discussion of the properties of the homogenized coefficients and Section 6 contains the conclusions.

Basic model

In this section we derive the equations on the microscopic level as a preparation for the homogenization procedure. First we present in Section 2.1, for the ease of the reader, a short description of the thermoelastic equations as they appear in the literature (for instance in the classical monographs by Ziegler [START_REF] Ziegler | An introduction to thermodynamics[END_REF] or Gurtin et al [START_REF] Gurtin | The mechanics and thermodynamics of continua[END_REF]).

Then in Section 2.2 we set up the equations and interface conditions at the microscale. In Section 2.3 we discuss the scaling of the microscale problem, yielding the dimensionless formulation and the homogenization parameter ε.

Thermo-elastic equations

Let e = e(u) = 1 2 (∇u + ∇ τ u)
be the symmetric (infinitesimal) strain tensor. Throughout this paper we use the notation

A : B = i,j A ij B ij = tr(AB τ )
for any two matrices A and B.

The fundamental equations describing the displacement (u) and the temperature (T ) of an elastic body are derived from the following considerations.

Momentum balance equation:

-

divσ = b, (1)
where σ denotes the (Cauchy) stress tensor and b a given body force. Often this is the gravitational force b = ρg with ρ denoting the mass density of the material and g the gravity vector. In this balance equation we have neglected inertia effects (ρ ∂ 2 u ∂t 2 ). Energy balance equation:

∂E ∂t = σ : e( ∂ ∂t u) -divq + h, (2)
q = -κ∇T.

Here E, q, κ and h denote, respectively, the internal energy of the material, heat flux, thermal conductivity tensor and external heat source. The dissipation term σ : e( ∂ ∂t u) will be removed from equations ( 2) and (3) when working in terms of the entropy (S).

Next we introduce the free energy

ψ = E -ST, (4)
where S denotes the entropy.

In the theory of thermo-elasticity a general assumption is (see for instance [START_REF] Ziegler | An introduction to thermodynamics[END_REF]) that

ψ = ψ(e, T ) (5)
only, with

σ = ∂ψ ∂e (σ ij = ∂ψ ∂e ij ) (6)
and

S = - ∂ψ ∂T . (7)
With this in mind, we first rewrite the energy balance [START_REF] Allaire | One-Phase Newtonian Flow[END_REF]. Differentiating [START_REF] Auriault | Etude du comportement macroscopique d'un milieu poreux saturé déformable[END_REF]

Differentiating now (4) with respect to time gives

∂ψ ∂t = ∂E ∂t -S ∂T ∂t -T ∂S ∂t . (9)
Equations ( 8) and ( 9) imply

T ∂S ∂t = ∂E ∂t -σ : e( ∂ ∂t u). (10)
Using this in (2) gives the rewritten energy balance

T ∂S ∂t = -divq + h. (11)
The free energy is expanded at a reference state which is strain free (e = 0) and has a reference temperature T 0 . Up to second order this gives

ψ(e, T ) = ψ 0 + ∂ψ ∂e | 0 : e + ∂ψ ∂T | 0 (T -T 0 ) + 1 2 e : ∂ 2 ψ ∂e 2 | 0 e + (T -T 0 ) ∂ 2 ψ ∂e∂T | 0 : e + 1 2
∂ 2 ψ ∂T 2 | 0 (T -T 0 ) 2 . ( 12
)
Assuming that the reference state is stress free and has zero entropy (i.e. assuming σ(0, T 0 ) and S(0, T 0 ) = 0, see ( 6)-( 7)), and neglecting the constant ψ 0 , expansion [START_REF] Coussy | From mixture theory to Biot's approach for porous media[END_REF] reduces to

ψ(e, T ) = 1 2 e : ∂ 2 ψ ∂e 2 | 0 e + (T -T 0 ) ∂ 2 ψ ∂e∂T | 0 : e + 1 2
∂ 2 ψ ∂T 2 | 0 (T -T 0 ) 2 . ( 13
)
In classical linear thermo-elasticity one takes for (13) the specific expression

ψ(e, T ) = 1 2 e : Ae + (T -T 0 )M : e - ρc 2T 0 (T -T 0 ) 2 , ( 14
)
where A is the elasticity tensor (generalized Hooke's law) satisfying

A ijkl = A jikl = A ijlk = A klij and Ae : e > 0 ( 15
)
and where M is the symmetric stress-temperature tensor. Finally, ρ denotes the density of the material and c the specific heat.

In the case of an isotropic material one has

A ijkl = λδ ij δ kl + µ(δ ik δ jl + δ il δ jk ) (16)
implying

σ ij = λδ ij k e kk + 2µe ij . (17)
Here λ and µ are the Lamé moduli. They are given by

λ = Λν (1 + ν)(1 -2ν) and µ = Λ 2(1 + ν) , ( 18
)
where Λ is Young's modulus and ν is Poisson's ratio. In the isotropic case, the stress-temperature tensor takes the form

M ij = -(3λ + 2µ)βδ ij , ( 19
)
where β is the coefficient of thermal expansion.

With expression [START_REF] Donea | An arbitrary Lagrangian-Eulerian finite element method for transient dynamic fluid-structure interactions[END_REF] one finds for the stress and entropy the expressions (the constitutive equations)

σ = ∂ψ ∂e = Ae + (T -T 0 )M (20)
and

S = - ∂ψ ∂T = -M : e + ρc T 0 (T -T 0 ). ( 21
)
This implies for the momentum balance

-div Ae(u) + (T -T 0 )M = b, (22)
and for the energy balance, replacing T ∂S ∂t by T 0 ∂S ∂t in [START_REF] Coussy | Mechanics of porous continua[END_REF],

∂ ∂t ρc(T -T 0 ) -T 0 M : e = div(κ∇T ) + h. (23)
In this work we confine ourselves to isotropic solid materials. This makes the derivation of the upscaling procedure and its interpretation more transparent. For the non-isotropic case we refer to Remark 8. The isotropic form of equations ( 22) and ( 23) reads

-div Ae(u) -(3λ + 2µ)β(T -T 0 )I = b (24)
and

∂ ∂t ρc(T -T 0 ) + (3λ + 2µ)βT 0 divu = div(κ∇T ) + h. ( 25
)
Introducing the dimensionless temperature

ϑ = T -T 0 ∆T , ∆T = characteristic temperature difference, (26)
into equations ( 24) and ( 25) we arrive at

-div Ae(u) -(3λ + 2µ)β∆T ϑI = b (27)
and

∂ ∂t ρcϑ + (3λ + 2µ) βT 0 ∆T divu = div(κ∇ϑ) + h ∆T . (28)
Equations ( 27)-( 28) are used to describe the displacement and temperature in the solid part of the porous structure.

Equations at the micro scale

Let us consider a geometrical structure as in Figure 1. This structure is repeated in all three directions to form a pavement Υ of the whole space. The characteristic size of the structure is denoted by (characteristic micro-scale length or pore size). The structure is decomposed as Υ = Υ f ∪ Υ s ∪ Γ Υ , where Υ f (the pore space) and Υ s (the solid skeleton) are open and connected sets and where Γ Υ is the interface between them.

Let Ω be a smooth domain in R 3 with characteristic length L. Then we write Ω = Ω p f ∪ Ω p s ∪ Γ p , with Ω p s = Ω ∩ Υ s : Solid part (skeleton), connected domain, behavior is described by the thermoelastic equations ( 27)- [START_REF] Rajagopal | On the Oberbeck-Boussinesq approximation[END_REF].

Ω p f = Ω ∩ Υ f : Fluid part of Ω,
Γ p = Ω ∩ Γ Υ : Fluid-solid interface with sufficient smoothness.

When describing the equations in Ω we use the Lagrangian approach in the solid part Ω p s and the so-called ALE (Arbitrary Lagrange-Eulerian) approach ([14], [START_REF] Hughes | Lagrangian-Eulerian finite element formulation for incompressible viscous flows[END_REF]) in the fluid part Ω p f . This means that the domains Ω p f and Ω p s (and the interface Γ p ) are fixed in time and that the fluid velocity v f in Ω p f is written with respect to the velocity of the solid structure. To this end, we extend the displacement vector u = u(x, t), with x ∈ Ω p s , of the solid structure to the whole of Ω, e.g. see Richter and Wick [START_REF] Richter | Finite elements for fluid-structure interaction in ALE and fully Eulerian coordinates[END_REF], and write the fluid velocity as

v f = w + ∂ ∂t u in Ω p f . (29)
If the interface Γ p acts as a no-slip boundary for the fluid we impose v f = ∂ ∂t u on Γ p or, equivalently,

w| Γ p = 0. ( 30
)
Note that in the ALE formulation one uses the fluid velocity v f in the flow equations, but the relative velocity w in the energy transport equation, see for instance ( [START_REF] Allaire | Ion transport through deformable porous media: derivation of the macroscopic equations using upscaling[END_REF], [START_REF] Richter | Finite elements for fluid-structure interaction in ALE and fully Eulerian coordinates[END_REF] and [START_REF] Wick | Fluid-Structure Interactions using Different Mesh Motion Techniques[END_REF]).

Equations in the fluid domain Ω p f

In Ω p f we describe the fluid movement by Stokes' equations and the Boussinesq approximation, see ( [START_REF] Farina | On the equations governing the flow of mechanically incompressible, but thermally expansible, viscous fluids[END_REF] , [START_REF] Gallavotti | Foundations of Fluid Dynamics[END_REF] , [START_REF] Rajagopal | On the Oberbeck-Boussinesq approximation[END_REF], [START_REF] Kh | Modélisation Asymptotique en Mécanique des Fluides Newtoniens[END_REF]). Let the fluid stress σ f be given by

σ f = -pI + 2η(ϑ)e(v f ), ( 31
)
where p is the fluid pressure and η(ϑ) is the fluid viscosity in terms of the dimensionless temperature [START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF]. Furthermore, let the fluid density ρ f be given by

ρ f = ρ 0 f (1 -β f (T -T 0 )) = ρ 0 f (1 -β f ∆T ϑ) (32)
where ρ 0 f is the reference mass density of the fluid and β f the coefficient of thermal expansion. Then from the fluid momentum balance

-divσ f = ρ f g = -ρ f ge 3 , ( 33
)
where g is the gravity constant and e 3 the unit vector pointing upwards in positive

x 3 direction, we obtain ∇p -div(2η(ϑ)e(v f )) = -ρ 0 f ge 3 + ρ 0 f β f ∆T gϑe 3 . (34)
Absorbing the hydrostatic pressure into p by redefining

p := p + ρ 0 f gx 3 , (35)
gives for the pressure and velocity of the fluid

∇p -div(2η(ϑ)e(v f )) = ρ 0 f β f ∆T gϑe 3 . (36)
Incompressibility requires

divv f = div(w + ∂ ∂t u) = 0. ( 37
)
Finally the fluid temperature satisfies (same notation for fluid and solid temperature)

ρ 0 f c f ( ∂T ∂t + w • ∇T ) = div(κ f ∇T ) + h f . ( 38
)
Here c f is the specific heat of the fluid, κ f the thermal conductivity tensor and h f the external heat source. In terms of the dimensionless temperature ( 26) we obtain

ρ 0 f c f ( ∂ϑ ∂t + w • ∇ϑ) = div(κ f ∇ϑ) + h f ∆T . (39)
Remark 1. In the case of large solid deformations, the flow equations would contain coefficients depending on the displacement u (e.g., see [START_REF] Richter | Finite elements for fluid-structure interaction in ALE and fully Eulerian coordinates[END_REF] and [START_REF] Wick | Fluid-Structure Interactions using Different Mesh Motion Techniques[END_REF]). In the present ALE description this is disregarded.

2.2.2

Equations in the solid domain Ω p s They are given by equations ( 27) and [START_REF] Rajagopal | On the Oberbeck-Boussinesq approximation[END_REF]. For completeness and because of adjusted notation (ρ → ρ s , c → c s , β → β s , etc) they are repeated here:

-divσ s = -div(Ae(u) -(3λ + 2µ)β s ∆T ϑI) = b (40)
and

∂ ∂t (ρ s c s ϑ + (3λ + 2µ) β s T 0 ∆T divu) = div(κ s ∇ϑ) + h s ∆T . (41)

Interface conditions at Γ p

At Γ p we impose (30) for the relative velocity w. The balance of momentum gives

σ f n = σ s n at Γ p , (42)
where the fluid pressure is replaced by p -ρ 0 f gx 3 . For the temperature we impose

[ϑ] = 0 and κ f ∇ϑ • n = κ s ∇ϑ • n at Γ p . ( 43
)
In conditions ( 42) and (43), n denotes the unit vector at Γ p pointing in the direction of the fluid.

In the remainder of this paper, this will always be the case.

Scaling of equations and conditions

To discuss the appropriate scaling we first introduce a number of characteristic quantities. Some are obvious, others are not. They will be chosen to obtain the correct balance in the equations and conditions. Dimensionless numbers that will appear are compared to ε, the ratio between the micro and the macro scale i.e.

ε = L (44)
where is the pore scale length and L the macro reservoir/experimental scale. Typically for dolomites and layered clays = O(10 -7 m), yielding a permeability from 1 to 10 millidarcy (md) (see for instance [START_REF] Bear | Dynamics of Fluids in Porous Media[END_REF]). Under the scaling x → x/L, Ω, Ω p f , Ω p s and Γ p become Ω, Ω ε f , Ω ε s and Γ ε , respectively. Note that after the scaling, the characteristic dimensionless pore size becomes ε.

Let p c , u c , τ c , v c , η c and Λ c denote, respectively, a characteristic pressure, deformation, time, fluid velocity and elastic constant. With these values one defines dimensionless variables in the obvious way: p := p/p c . . . etc. Throughout this paper we take Λ c = Λ (Young's modulus) and η c = η(0) (viscosity at reference temperature T 0 ).

Let us start with condition (42). Introducing the characteristic values and using [START_REF] Hughes | Lagrangian-Eulerian finite element formulation for incompressible viscous flows[END_REF] gives

(-p c (p -p h x 3 )I + 2η c v c L η(ϑ)e(w + u c v c τ c ∂ ∂t u)) • n = (Λ u c L Ae(u) - 4 -2ν (1 + ν)(1 -2ν) Λβ s ∆T ϑI) • n (45)
at Γ ε . Balancing the first term on the left and right gives

p c = Λ u c L and thus p h := ρ 0 f gL p c = ρ 0 f gL 2 Λu c . ( 46
)
where p h is a measure for the dimensionless hydrostatic pressure. Leaving u c and τ c unspecified for the moment, we take for v c the characteristic Darcy velocity

v c = 2 η c p c L = ε 2 Lp c η c = ε 2 Λu c η c =: u c τ T , (47)
where

τ T = η c ε 2 Λ [s]. ( 48
)
is the time scale related to the deformation of the medium in response to a change in pore pressure. Dividing (45) by p c from (46) gives

(-(p -p h x 3 )I + 2ε 2 η(ε)e(w + τ T τ c ∂ ∂t u)) • n = (Ae(u) -K S ϑI) • n at Γ ε . ( 49
)
Here

K S = 4 -2ν (1 + ν)(1 -2ν) β s ∆T L u c . ( 50
)
Using the same scaling in the momentum balance equation for the solid structure [START_REF] Zimmerman | Coupling in poroelasticity and thermoelasticity[END_REF] gives

-div(Ae(u) -K S ϑI) = L 2 Λu c b =: -K g e 3 in Ω ε s , (51)
where

K g = L 2 ρ s g Λu c = ρ s ρ 0 f p h . ( 52
)
Next we consider Stokes' equations [START_REF] Wick | Fluid-Structure Interactions using Different Mesh Motion Techniques[END_REF]. Applying the scaling with ( 46), ( 47) and (48) gives

∇p -2ε 2 div(η(ϑ)e(w + τ T τ c ∂ ∂t u)) = K B ϑe 3 in Ω ε f , (53)
where

K B = K g ρ 0 f ρ s β f ∆T. ( 54
)
The incompressibility condition becomes

div(w + τ T τ c ∂ ∂t u) = 0 in Ω ε f . ( 55
)
It remains to scale the temperature equations ( 39) and (41). These equations contain the heat capacities ρ 0 f c f and ρ s c s , both considered constant here. It is convenient at this point to introduce the average heat capacity ρc = Φρ

0 f c f + (1 -Φ)ρ s c s , (56)
where Φ = |Ω ε f |/|Ω| denotes the porosity of the porous structure. The scaled equation for the fluid temperature becomes

C f ∂ϑ ∂t + v c τ c L C f w • ∇ϑ = D f ∆ϑ + hf in Ω ε f , (57)
where

           C f = ρ 0 f c f ρc , D f = κ f τ c ρcL 2 = D f τ c L 2 hf = τ c ρc∆T h f . . (58)
Here D f = κ f / ρc denotes the dimensional thermal diffusivity of the fluid. Note that Péclet's number in (57) is

Pe := v c τ c LD f C f = ρ 0 f c f v c L κ f (59)
The equation for the temperature in the solid structure reads

∂ ∂t C s ϑ + K T div u = D s ∆ϑ + hs in Ω ε s , (60) with
                 C s = ρ s c s ρc , K T = 4 -2ν (1 + ν)(1 -2ν) Λβ s T 0 ρc∆T u c L , D s = κ s τ c ρcL 2 = D s τ c L 2 hs = τ c ρc∆T h s . . ( 61
)
Here D s = κ s / ρc denotes the dimensional thermal diffusivity of the solid. Note that

ΦC f + (1 -Φ)C s = 1.
The temperature flux condition (43) becomes

D f ∇ϑ • n = D s ∇ϑ • n at Γ ε .
Remark 2. When writing (57) and (60) we silently replaced the heat conduction tensors by positive constants κ f and κ s .

Remark 3. In the quasi-static setting we are mainly interested in a filtration dominated regime.

Then intuitively one would think of a decoupling, where first the filtration through a rigid porous structure is calculated and then the deformation of the solid skeleton is added a posteriori. In such a setting the characteristic Darcy velocity is given by (47), with p c being the pressure drop. It originates from either the gravity forces, with p c = ρ 0 f gL and u c = L 2 ρ 0 f g/Λ or from the applied pressure at the outer boundary. The linear elasticity equations hold if the relative displacement is sufficiently small, say

u c /L = τ T /τ c = p c /Λ < O(10 -2 ) .

Discussion on the choice of scales

First we recall in Table 1 some dimensional solid and fluid characteristics (order of magnitude). These values are taken from the literature, (e.g. [START_REF] Lee | Thermal consolidation in porous media by homogenization theory-I. Derivation of macroscale equations[END_REF], [START_REF] Mctigue | Thermoelastic Response of Fluid-Saturated Porous Rock[END_REF], [START_REF] Zimmerman | Coupling in poroelasticity and thermoelasticity[END_REF]). Table 1: Order of the magnitude of data

Motivation

We take L = 1m as the macro (Darcy) scale and = 10 -7 m as the pore scale. Note that this leads to a permeability of O(10 -14 m 2 ) which is representative for various rock and soil types. L.

Further we choose τ c = L/v c in (57) and (58). This gives for Péclet's number (59)

Pe = ρ 0 f c f κ f L 2 τ c .
Using (47) we note that

τ T τ c = u c L .
Remark 4. If we compare our assumptions to [START_REF] Lee | Thermal consolidation in porous media by homogenization theory-I. Derivation of macroscale equations[END_REF], it turns out that these authors assumed p c /Λ = O(ε 2 ), see their formula [START_REF] Mctigue | Thermoelastic Response of Fluid-Saturated Porous Rock[END_REF] on page 129. Such a strong assumption is not necessary in our approach.

With the choices and values from Table 1 we obtain for the other characteristic quantities and dimensionless numbers the values (with ∆T = 50 0 K).

SYMBOL QUANTITY ROCK VALUE CLAY VALUE τ T = η C /(Λε 2 )
characteristic deformation time 10 s 10 4 s p h = ρ 0 f gL/pc ratio hydrostatic vs pore pressures 10 -3 10 -1 vc = uc/τ T velocity of the structure deformation 10 -4 m/s 10 -6 m/s

K S = βs∆T L/uc solid expansion coefficient 5 • 10 -2 5 • 10 -3 Kg = ρsL 2 g/(Λuc) solid gravitational coefficient 2.5 • 10 -3 2.5 • 10 -1 K B = Kgρ 0 f β f ∆T /ρs buoyancy coefficient 5 • 10 -6 5 • 10 -4 τc transport time 10 4 s 10 6 s D f = D f τc/L 2 scaled fluid diffusivity 10 -3 10 -1 Ds = Dsτc/L 2 scaled solid diffusivity 10 -2 1 K T = ΛβsT 0 uc/( ρc∆T L) thermoelastic coupling parameter 3 • 10 -4 6 • 10 -7
Table 2: Data description Remark 5. We observe that for both type of materials the factor K T is very small and can be disregarded for any practical or computational purposes. Note that one can write

K T = p c β s T 0 ρc∆T .
Hence for considerably larger pressures this factor could gain importance. Therefore we keep it in the two-scale expansions. The smallness of the thermoelastic coupling parameter K T was already noticed in [START_REF] Mctigue | Thermoelastic Response of Fluid-Saturated Porous Rock[END_REF].

Remark 6. Let u ε denote the dimensionless displacement. Then

u ε = u 0 + εu 1 + . . .
Hence the actual displacement can be written as

u c u ε = u c u 0 macro displacement + u c L u 1 displacement in the pore + . . .
Observe that the pore displacement is much smaller than the macroscopic displacement.

3 Upscaling through two-scale expansions

Geometry

We start by some comments on the geometry:

• Porous media structures met in nature are hard to describe. In order to model their shape it is necessary to pose some hypotheses. In general it is supposed that there are two connected phases: a fluid phase and a solid phase, the latter being either rigid or deformable. Furthermore, the porous medium is supposed to be heterogeneous at the microscopic (pore) level but statistically homogeneous at the macroscopic level.

• The characteristic length of the scaled pores is ε. Since the theory for the true velocities, pressures and other quantities is too complicated, one prefers working with the averaged quantities over characteristic volumes being of order ε 3 .

A representative example of the geometry is the periodic porous medium with connected fluid and solid phases. It is obtained by a periodic arrangement of the pores. After dividing by the characteristic length L, the porous medium becomes a box Ω = (0, R)

3 in R 3 , R > 0.
The description goes along the following lines: First define the geometrical structure inside the unit cell Y = (0, 1) 3 . Let Y s (the solid part) be an open subset of Y and Y f = Y \ Ȳs (the fluid part). The liquid/solid interface is S = ∂Y s \ ∂Y . It is assumed, see for instance the unit cells in Figure 1, that Y s is such that (i) if we repeat Y s all over R 3 we obtain a periodic solid skeleton E S which is connected and (ii) the complement of ĒS in R 3 is the pore space E F which is connected as well.

Now Ω = (0, R) 3 is covered with a regular mesh of size ε, each cell being a cube

Y ε i = ε(Y + i), with 1 ≤ i ≤ N (ε) = |Ω|ε -3 [1 + 0(1)]. We define Y ε si = ε(Y s + i) and Y ε fi = ε(Y f + i). For sufficiently small ε > 0 we consider T ε = {k ∈ Z 3 |Y ε s k ⊂ Ω} and define Ω ε s = k∈Tε Y ε s k , Γ ε = ∂Ω ε s \ ∂Ω, Ω ε f = Ω \ Ω ε s .
The domains Ω ε s and Ω ε f represent, respectively, the solid and fluid parts of a porous medium Ω. For simplicity, we will suppose that L/ε is an integer. Note that |Y f | = Φ is the porosity of the reference configuration.

Remark 7. The sketched geometry is the reference configuration in the Lagrangian description of the deformation of the elasticity part.

Pore level equations

Recalling the equations from Subsection 2.3, we consider the dimensionless system of partial differential equations for t ∈ (0, T ), T > 0 :

-divσ f,ε = (K B ϑ ε -p h )e 3 in Ω ε f , (62) div
(w ε + τ T τ c ∂u ε ∂t ) = 0 in Ω ε f , ( 63
) divσ s,ε -K g e 3 = 0 in Ω ε s , (64)
w ε = 0 on Γ ε , (65)
σ f,ε = -(p ε -p h x 3 )I + 2ε 2 η(ϑ ε )e(w ε + τ T τ c ∂u ε ∂t ) in Ω ε f , (66)
σ s,ε = Ae(u ε ) -K S ϑ ε I in Ω ε s , (67)
σ s,ε n = σ f,ε n on Γ ε , (68)
C f ∂ϑ ε ∂t + w ε • ∇ϑ ε = D f ∆ϑ ε + hf in Ω ε f , (69)
∂ ∂t C s ϑ ε + K T div u ε = D s ∆ϑ ε + hs in Ω ε s , (70)
D f n • ∇ϑ ε = D s n • ∇ϑ ε and [ϑ ε ] = 0 on Γ ε . (71)
Note that τ c is the heat convection time (transport time) and τ T /τ c = u c /L.

The two-scale asymptotic expansion

Due to the geometrical structure, our problem has two natural scales: the macro-scale x (the core scale) and the micro-scale y = x ε . We search for {u ε , ϑ ε , w ε , p ε } in the following two-scale expansion form

u ε (x, t) = u 0 (x, y, t) + εu 1 (x, y, t) + ε 2 u 2 (x, y, t) + . . . , (72)
ϑ ε (x, t) = ϑ 0 (x, y, t) + εϑ 1 (x, y, t) + ε 2 ϑ 2 (x, y, t) + . . . , (73)
w ε (x, t) = w 0 (x, y, t) + εw 1 (x, y, t) + ε 2 w 2 (x, y, t) + . . . ,

p ε (x, t) = p 0 (x, y, t) + εp 1 (x, y, t) + ε 2 p 2 (x, y, t) + . . . ,

First we show that the effective displacement does not oscillate. To this end we substitute the expansion into equations ( 64) and (67) to get:

K g e 3 = ε -2 div y (Ae y (u 0 )) + ε -1 div y (Ae y (u 1 ) -K S ϑ 0 I + Ae x (u 0 )) + div x (Ae y (u 0 )) + ε 0 div y (Ae y (u 2 ) + Ae x (u 1 ) -K S ϑ 1 I) + div x (Ae y (u 1 ) -K S ϑ 0 I+ Ae x (u 0 ) + O(ε) in Y s × Ω × (0, T ). ( 76
)
Expanding interface condition (68) gives

σ s,ε • n| Γ ε = ε -1 Ae y (u 0 ) • n + (Ae y (u 1 ) -K S ϑ 0 I + Ae x (u 0 )) • n+ ε(Ae y (u 2 ) -K S ϑ 1 I + Ae x (u 1 )) • n | Γ ε = σ f,ε • n| Γ ε = O(1). (77)
At the order O(ε -2 ), for any (x, t) ∈ Ω × (0, T ) we have

div y Ae y (u 0 ) = 0 in Y s , Ae y (u 0 ) • n = 0 on S, u 0 is 1-periodic. ( 78
)
Since the solid structure is connected, as we are in R 3 , the only possible rigid motion is a constant translation with respect to y. Hence we have

u 0 (x, y, t) = u 0 (x, t) in Y s × Ω × (0, T ). ( 79
)
Therefore the displacement of the solid structure doesn't depend on the fast variable at order zero. Its extension to the pore space is also non-oscillatory.

STEP 2

Next we show that the effective pressure does not oscillate in the pores and write the Stokes system for w 0 in the pores. Insert now the two-scale expansions into equations ( 62), ( 63) and ( 65)-(66). It yields

σ f,ε = -(p 0 -p h x 3 )I + ε -p 1 + 2η(ϑ 0 )e y (w 0 + τ T τ c ∂u 0 ∂t ) + O(ε 2 ) in Y f × Ω × (0, T ), (80) ε -1 ∇ y p 0 + ε 0 ∇ y p 1 + ∇ x (p 0 -p h x 3 ) -2 div y (η(ϑ 0 )e y (w 0 )) -(K B ϑ 0 -p h )e 3 = O(ε) (81) in Y f × Ω × (0, T ). At order O(ε -1 ) we get p 0 (x, t, y) = p 0 (x, t) in Y f × Ω × (0, T ). ( 82
)
Extending the pressure by p 0 = 0 in Ȳs × Ω × (0, T ) we have

p 0 (x, t, y) = χ Y f (y)p 0 (x, t),
where χ Y f denotes the characteristic function of the set Y f . Furthermore, the mass conservation equation becomes

ε -1 div y w 0 + div x w 0 + τ T τ c ∂u 0 ∂t + div y w 1 + τ T τ c ∂u 1 ∂t + O(ε) = 0 (83)
and at order O(ε 0 ) we have the following Stokes problem

∇ y p 1 + ∇ x (p 0 -p h x 3 ) -2 div y (η(ϑ 0 )e y (w 0 )) = (K B ϑ 0 -p h )e 3 in Y f × Ω × (0, T ), ( 84
)
div y w 0 = 0 in Y f × Ω × (0, T ), ( 85
)
w 0 = 0 on S × Ω × (0, T ), ( 86
)
{w 0 , p 1 } are 1 -periodic in y. (87)
Comparison between (77) and (80) yields

ε -1 Ae y (u 0 ) • n + (Ae y (u 1 ) -K S ϑ 0 I + Ae x (u 0 ) + (p 0 -p h x 3 ))I) • n + ε(Ae y (u 2 ) -K S ϑ 1 I+ Ae x (u 1 ) + p 1 I -2η(ϑ 0 )e y (w 0 )) • n + O(ε 2 ) = 0 on S × Ω × (0, T ). ( 88
)
Now we are able to go back to the expansion for the elasticity system and study the problem for u1 .

STEP 3

Using (76) and (88) yields at order O(ε -1 )

       div y A(e y (u 1 ) + e x (u 0 ) -K S ϑ 0 I) + div x (A e y (u 0 ) =0 ) = 0 in Y s × Ω × (0, T ),
A(e y (u 1 ) + e x (u 0 )) At this point we need information concerning ϑ 0 . Inserting the two-scale expansions into equations (69)-(71) gives

-K S ϑ 0 I • n = -(p 0 -p h x 3 )n on S × Ω × (0, T ), u 1 is 1-periodic.
-ε -2 D f ∆ y ϑ 0 + ε -1 C f w 0 • ∇ y ϑ 0 -D f ∆ y ϑ 1 -2D f div x ∇ y ϑ 0 + ε 0 C f ∂ϑ 0 ∂t + C f (w 0 • ∇ y ϑ 1 + w 0 • ∇ x ϑ 0 + w 1 • ∇ y ϑ 0 ) -D f ∆ y ϑ 2 -2D f div x ∇ y ϑ 1 -D f ∆ x ϑ 0 -hf = O(ε) in Y f × Ω × (0, T ), ( 90
) -ε -2 D s ∆ y ϑ 0 + ε -1 K T div y u 0 -D s ∆ y ϑ 1 -2D s div x ∇ y ϑ 0 + ε 0 C s ∂ϑ 0 ∂t + K T ( div y u 1 + div x u 0 ) -D s ∆ y ϑ 2 -2D s div x ∇ y ϑ 1 -D s ∆ x ϑ 0 -hs = O(ε) in Y s × Ω × (0, T ), (91)
ε -1 D f ∇ y ϑ 0 • n -D s ∇ y ϑ 0 • n + D f (∇ y ϑ 1 + ∇ x ϑ 0 ) • n -D s (∇ y ϑ 1 + ∇ x ϑ 0 ) • n+ ε D f (∇ y ϑ 2 + ∇ x ϑ 1 ) • n -D s (∇ y ϑ 2 + ∇ x ϑ 1 ) • n = O(ε 2 ) on S × Ω × (0, T ). ( 92
)
We recall that ϑ j is continuous across S × Ω × (0, T ) for every j. The leading order terms in the two-scale expansions (90)-(92) yield the following problem at order O(ε -2 ):

-D f ∆ y ϑ 0 = 0 in Y f , (93)
-D s ∆ y ϑ 0 = 0 in Y s , (94)
D f ∇ y ϑ 0 • n = D s ∇ y ϑ 0 • n and [ϑ 0 ] = 0 on S, ( 95
)
ϑ 0 is 1 -periodic, (96)
for every (x, t) ∈ Ω × (0, T ).

Only possible solution for problem (93)-( 96) is a constant in y. Hence ϑ 0 = ϑ 0 (x, t).

The next order terms in the expansions (90)-(92) yield the following problem at order O(ε -1 ):

-D f div y (∇ y ϑ 1 + ∇ x ϑ 0 ) + C f w 0 • ∇ y ϑ 0 =0 -D f div x ∇ y ϑ 0 =0 = 0 in Y f , (97)
K T ∂ ∂t div y u 0 =0 -D s div y (∇ y ϑ 1 + ∇ x ϑ 0 ) -D s div x ∇ y ϑ 0 =0 = 0 in Y s , (98)
D f (∇ y ϑ 1 + ∇ x ϑ 0 ) • n = D s (∇ y ϑ 1 + ∇ x ϑ 0 ) • n and [ϑ 1 ] = 0 on S, ( 99
)
ϑ 1 is 1 -periodic, ( 100
)
for every (x, t) ∈ Ω × (0, T ). Fredholm's alternative implies that problem (89) has a solution. This solution is unique if we impose

Y ϑ 1 dy = 0.
For ϑ 1 a separation of scales is possible, which is achieved through the temperature auxiliary problems, for j = 1, 2, 3,

-D f div y (∇ y Θ j + e j ) = 0 in Y f , (101)
-D s div y (∇ y Θ j + e j ) = 0 in Y s , (102)
D f (∇ y Θ j + e j ) • n = D s (∇ y Θ j + e j ) • n and [Θ j ] = 0 on S, (103)
Θ j is 1 -periodic, Y Θ j dy = 0. ( 104
)
Then ϑ 1 is given by the scale separation formula

ϑ 1 (x, y, t) = 3 j=1 Θ j (y) ∂ϑ 0 ∂x j (x, t). ( 105
)
In order to have an equation for ϑ 0 , we need the next order terms in the expansions (90)-(92). They yield the following problem at order O(1):

-D f div y (∇ y ϑ 2 + ∇ x ϑ 1 ) + C f ∂ϑ 0 ∂t + C f (w 1 • ∇ y ϑ 0 =0 + w 0 • (∇ y ϑ 1 + ∇ x ϑ 0 )) -D f div x ∇ y ϑ 1 -D f ∆ x ϑ 0 = hf in Y f , (106)
-D s div y (∇ y ϑ 2 + ∇ x ϑ 1 ) -D s div x ∇ y ϑ 1 -D s ∆ x ϑ 0 + C s ∂ϑ 0 ∂t + K T ∂ ∂t ( div y u 1 + div x u 0 ) = hs in Y s , (107)
D f (∇ y ϑ 2 + ∇ x ϑ 1 ) • n = D s (∇ y ϑ 2 + ∇ x ϑ 1 ) • n and [ϑ 2 ] = 0 on S, ( 108
)
ϑ 2 is 1 -periodic, ( 109
)
for every (x, t) ∈ Ω × (0, T ). By Fredholm's alternative, problem (106)-( 109) admits a solution if and only if 110) Equation ( 110) can be useful only if we are able to eliminate ϑ 1 . Substituting (105) into the fourth term (divergence of the effective heat flux) of ( 110) yields

K T ∂ ∂t ( Ys div y u 1 dy + (1 -Φ) div x u 0 ) + (C f Φ + C s (1 -Φ)) =1 ∂ϑ 0 ∂t + C f Y f w 0 • (∇ y ϑ 1 + ∇ x ϑ 0 ) dy- div x D f Y f ∇ y ϑ 1 dy + D s Ys ∇ y ϑ 1 dy + (D f Φ + D s (1 -Φ))∇ x ϑ 0 = hf Φ + hs (1 -Φ). (
D f Y f ∇ y ϑ 1 dy + D s Ys ∇ y ϑ 1 dy + (D f Φ + D s (1 -Φ))∇ x ϑ 0 = 3 j=1 (D f Y f ∇ y Θ j dy+ D s Ys ∇ y Θ j dy) ∂ϑ 0 ∂x j (x, t) + (D f Φ + D s (1 -Φ))∇ x ϑ 0 =: B∇ x ϑ 0 ,
where

B ij = D f Y f ∂ ∂y i Θ j dy + D s Ys ∂ ∂y i Θ j dy + (D f Φ + D s (1 -Φ))δ ij , i, j = 1, 2, 3. (111)
B is the effective heat dispersivity. It is easy to see that it is a symmetric positive definite matrix (see for instance [START_REF] Sanchez-Palencia | Non-Homogeneous Media and Vibration Theory[END_REF] or [START_REF] Mei | Homogenization methods for multiscale mechanics[END_REF]). We also eliminate ϑ 1 from the transport term in (110). This gives

Y f w 0 • (∇ y ϑ 1 + ∇ x ϑ 0 ) dy = ( Y f w 0 dy) • ∇ x ϑ 0 + Y f w 0 • ∇ y ϑ 1 dy = ( Y f w 0 dy) • ∇ x ϑ 0 - Y f div y w 0 =0 ϑ 1 dy + S w 0 • n =0 ϑ 1 dS y = ( Y f w 0 dy) • ∇ x ϑ 0 .
Hence equation (110) takes the form

K T ∂ ∂t ( Ys div y u 1 dy + (1 -Φ) div x u 0 ) + ∂ϑ 0 ∂t + C f ( Y f w 0 dy) • ∇ x ϑ 0 - div x B∇ x ϑ 0 = hf Φ + hs (1 -Φ). ( 112
)
After evaluating the average over Y s of div y u 1 , it will represent the homogenized equation for the temperature.

STEP 5

Now we return to equation (89) for the solid displacement. Since this problem is linear, we decompose u 1 through the scale separation formula

u 1 (x, y, t) = 3 i,j=1 (e x (u 0 (x, t))) ij a ij (y) + (p 0 (x, t) -p h x 3 -K S ϑ 0 (x, t))b 0 (y), ( 113
)
Here the displacement fields a ij and b 0 depend only on the geometry and are defined through the auxiliary problems for i, j = 1, 2, 3 :

   div y A(e y (a ij ) + 1
2 (e i ⊗ e j + e j ⊗ e i )) = 0 in Y s , A(e y (a ij ) + 1 2 (e i ⊗ e j + e j ⊗ e i ))

• n = 0 on S, a ij is 1-periodic (114) and div y A(e y (b 0 ) = 0 in Y s , Ae y (b 0 ) • n = -n on S, b 0 is 1-periodic. ( 115
)
With the auxiliary problems defining u 1 we associate the linear poroelasticity effective coefficients. They are given by the Gassman tensor

G klij = Ys A( 1 2 (e i ⊗ e j + e j ⊗ e i ) + e y (a ij )) dy kl , i, j, k, l = 1, 2, 3, (116)
and the Biot matrix

B H = Ys Ae y (b 0 ) dy, (117)
with a ij and b 0 defined by the problems ( 114) and ( 115), respectively. The tensor G and the matrix B H are symmetric and G is positive definite, see ( [START_REF] Sanchez-Palencia | Non-Homogeneous Media and Vibration Theory[END_REF], p. 129-190) .

Next we go back to the equations for the displacement (76), (88). More precisely we consider the terms of order O(1) in (76) and of order O(ε) in ( 77) and (80). This results in the following boundary value problem for u 2 -div y (Ae y (u 2 ) + Ae x (u 1 ) -K S ϑ 1 I) -div x (Ae y (u 1 ) + Ae x (u 0 ) -K S ϑ 0 I) = -K g e 3 in Y s × Ω × (0, T ), ( 118)

-(A(e y (u 2 ) + e x (u 1 )) -K S ϑ 1 I) • n = p 1 n -2η(ϑ 0 )e y (w 0 ) • n on S × Ω × (0, T ), ( 119
)
u 2 is 1-periodic in y. ( 120
)
Fredholm's alternative gives the compatibility condition

div x {(1 -Φ)(Ae x (u 0 ) -K S ϑ 0 I) + Ys Ae y (u 1 ) dy} -K g (1 -Φ)e 3 = S (p 1 I -2η(ϑ 0 )e y (w 0 )) • n dS y . Since S (p 1 I -2η(ϑ 0 )e y (w 0 )) • n dS y = - Y f (∇ y p 1 -2 div y (η(ϑ 0 )e y (w 0 ))) dy,
(note that the normal vector n points from the solid part to the fluid) and using equation (84) yields

div x {(1 -Φ)(Ae x (u 0 ) -K S ϑ 0 I) + Ys Ae y (u 1 ) dy} -K g (1 -Φ)e 3 = Y f (-K B ϑ 0 e 3 + ∇ x p 0 ) dy. (121)
To use this condition we have to evaluate the integral over Y s . Here we use the decomposition formula (113):

Ys {A(e y (u 1 ) + e x (u 0 )) -K S ϑ 0 I} dy

= 3 i,j=1 e x (u 0 (x, t))) ij A(a ij (y)+ 1 2 (e i ⊗ e j + e j ⊗ e i ) dy + (p 0 (x, t) -p h x 3 ) Ys Ab 0 (y) -K S ϑ 0 (x, t) Ys (Ab 0 (y) + I) dy = Ge x (u 0 ) + (p 0 -p h x 3 )B H -K S ϑ 0 B H -K S |Y s |ϑ 0 I. (122)
Writing equation (121) in the form

div x { Ys (Ae x (u 0 ) + Ae y (u 1 ) -K S ϑ 0 I) dy -Φp 0 I} = -K B Φϑ 0 e 3 + K g (1 -Φ)e 3 . (123)
and substituting (122) into (123) gives

div x {Ge x (u 0 ) + (p 0 -p h x 3 )B H -K S ϑ 0 B H -K S |Y s |ϑ 0 I -Φp 0 I} = -K B Φϑ 0 e 3 + K g (1 -Φ)e 3 .
Hence we have obtained the system of partial differential equations for the effective solid displacement u 0 :

-div x {Ge x (u 0 )} + div x {(ΦI -B H )(p 0 -p h x 3 ) + K S ϑ 0 (B H + (1 -Φ)I)} = K B Φϑ 0 e 3 -p h Φe 3 -K g (1 -Φ)e 3 . (124)
STEP 6

Now we return to the problem (84)-(87). Since p 0 and ϑ 0 do not depend on y, we can use the scale separation from the classical derivation of the Darcy law via homogenization (see e.g. [START_REF] Allaire | One-Phase Newtonian Flow[END_REF] or [START_REF] Sanchez-Palencia | Non-Homogeneous Media and Vibration Theory[END_REF]). It uses the well known auxiliary problems in the homogenization of the filtration through rigid porous media and gives

w 0 (x, y, t) = 3 j=1 q j (y) η(ϑ 0 ) (K B ϑ 0 δ j3 - ∂p 0 ∂x j (x, t)), (125)
p 1 (x, y, t) = 3 j=1 π j (y)(K B ϑ 0 δ j3 - ∂p 0 ∂x j (x, t)), (126)
where {q j , π j } satisfy -∆ y q j + ∇ y π j = e j and div y q j = 0 in Y f , (127)

q j = 0 on S, {q j , π j } is 1-periodic in y. (128)
Then

K ij := Y f q i j (y) dy (129)
is the classical permeability tensor in Darcy's law, which is symmetric and positive definite (see for instance [START_REF] Sanchez-Palencia | Non-Homogeneous Media and Vibration Theory[END_REF]).

Consequently, the Darcy law reads

v D = Y f w 0 dy = - K η(ϑ 0 ) ∇ x p + (p h -K B ϑ)e 3 .
It remains to find an effective equation for the pressure p 0 . We start with the equation at order O(1) in the expansion of mass conservation condition (83). It reads

div x w 0 + τ T τ c ∂u 0 ∂t + div y w 1 + τ T τ c ∂u 1 ∂t = 0. Integration of this equation in Y f yields div x Y f w 0 dy + Φ τ T τ c ∂u 0 ∂t = - τ T τ c Y f div y ∂u 1 ∂t dy = τ T τ c ∂ t Ys div y u 1 (x, y, t) dy. (130)
Recalling, see [START_REF] Mikelić | On the interface law between a deformable porous medium containing a viscous fluid and an elastic body[END_REF],

B H ij = Ys (Ae y (b 0 )(y)) ij dy = 1 2 Ys Ae y (b 0 )(y) : (e i ⊗ e j + e j ⊗ e i ) dy = - Ys Ae y (b 0 )(y) : e y (a ij ) dy = Ys div y a ij (y) dy, (131)
introducing the coefficient, see again [START_REF] Mikelić | On the interface law between a deformable porous medium containing a viscous fluid and an elastic body[END_REF] 1 M := - (e x (u 0 (x, t))) ij Ys div y a ij (y) dy+

(p 0 (x, t) -p h x 3 -K S ϑ 0 (x, t)) Ys div y b 0 (y) dy, (133)
we replace Ys div y a ij (y) dy by B H ij and Ys div y b 0 (y) dy by -1/M . After inserting (133) into (130), we find the the pressure equation of the quasi-static thermal Biot model:

∂ t 1 M p 0 - K S M ϑ 0 + (ΦI -B H ) : e x (u 0 ) + τ c τ T div x K η(ϑ 0 ) (K B ϑ 0 e 3 -∇ x p 0 ) = 0. ( 134
) STEP 7
After the computation from the previous step, we are able to write more explicitly the temperature equation. We substitute (133) into (112) to finally obtain:

K T ∂ ∂t - 1 M (p 0 -p h x 3 ) + (B H + (1 -Φ)I) : e x (u 0 ) + (1 + K S K T M ) ∂ϑ 0 ∂t + C f ( Y f w 0 dy) • ∇ x ϑ 0 -div x B∇ x ϑ 0 = hf Φ + hs (1 -Φ). ( 135
)
Remark 8. (Non-isotropic media): In case of non-isotropic tensors A and M one uses equations ( 22) and ( 23) as starting point for the two-scale expansions. Then the upscaled temperature equation ( 112) is replaced by the expression

K T ∂ ∂t Ys M : e y (u 1 ) dy + (1 -Φ)M : e x (u 0 ) + ∂ϑ 0 ∂t + C f ( Y f w 0 dy) • ∇ x ϑ 0 - div x B∇ x ϑ 0 = hf Φ + hs (1 -Φ). ( 136
)
The scale separation formula for u 1 becomes more involved. An additional auxiliary problem of the form div y A(e y (b

1 ) = 0 in Y s , Ae y (b 1 ) • n = Mn on S, b 1 is 1-periodic, (137)
is needed. This is used in the separation formula which now reads

u 1 (x, y, t) = 3 i,j=1 (e x (u 0 (x, t))) ij a ij (y) + (p 0 (x, t) -p h x 3 )b 0 (y) + K S ϑ 0 (x, t)b 1 (y). ( 138
)
Then arguing as is Step 5, we obtain the modified (non-isotropic) Biot equation

-div x {Ge x (u 0 )} + div x {(ΦI -B H )(p 0 -p h x 3 ) -K S ϑ 0 (M H -(1 -Φ)M)} = K B Φϑ 0 e 3 -p h Φe 3 -K g (1 -Φ)e 3 , (139)
where M H = Ys Ae y (b 1 ) dy. Similarly, the mass conservation equation (134) changes into

∂ t p 0 M + K S ϑ 0 ( Ys div y b 1 (y) dy) + (ΦI -B H ) : e x (u 0 ) + τ c τ T div x K η(ϑ 0 ) (K B ϑ 0 e 3 -∇ x p 0 ) = 0. ( 140
)
Finally, the modified temperature equation is obtained by substituting expression (138) into equation (136). Details are omitted.

The upscaled model and a corresponding Lyapunov functional

Let the Darcy velocity (i.e. the upscaled fluid velocity relative to the solid skeleton) be given by

v D = Y f w 0 dy. ( 141
)
With this definition the homogenization procedure from Section 3 gives the following set of equations (for (x, t) ∈ Ω × (0, T )) : Darcy's law

v D = - K η(ϑ 0 ) ∇ x p + (p h -K B ϑ)e 3 , (142)
where p denotes the dimensionless actual fluid pressure (i.e. including the hydrostatic component) and p h -K B ϑ is proportional to the dimensionless density of the fluid phase. Momentum balance

-div x σ = K B Φϑe 3 -p h Φe 3 -K g (1 -Φ)e 3 , (143)
where the total stress σ is given by

σ = Ge x (u) -(ΦI -B H )p -K S ϑ(B H + (1 -Φ)I). ( 144
)
Mass balance

∂ t ζ + τ c τ T div x v D = 0, ( 145
)
where the change in volume fluid content ζ is given by

ζ = p -K S ϑ M + (ΦI -B H ) : e x (u); (146)
Energy balance

K T ∂ ∂t -ζ + div x u + ∂ϑ ∂t + C f v D • ∇ x ϑ -div x B∇ x ϑ = hf Φ + hs (1 -Φ). (147)
Note that the upscaled system (142)-(147) preserves the basic conservation law properties of the micro-scale system. Furthermore, it has the following important stability property:

Proposition. Let 2W (e, p, ϑ) = Ge : e + 1 M (p -K S ϑ) 2 + K S K T ϑ 2 , ( 148
)
where e is taking values in the space of symmetric square matrices on R 3 and p and ϑ are scalars. Then, in the absence of the gravity and other external forces and for small ∆T , W is a Lyapunov functional for the solutions of the system (142)-(147), i.e. Finally, equation ( 147) is multiplied by K S /K T and then tested by ϑ. This gives

-K S Ω ∂ t p -K S ϑ M ϑ dx + K S Ω (B H + (1 -Φ)I) : e x (∂ t u)ϑ dx + K S K T Ω ∂ t ϑϑ dx+ K S K T Ω C f v D • ∇ x ϑϑ dx + K S K T Ω B∇ x ϑ • ∇ x ϑ dx = 0. ( 151
)
Summing up the above three equations, we observe that most of the cross terms cancel and there remains

1 2 ∂ t Ω {Ge x (u) : e x (u) + p 2 M - 2K S pϑ M + ( K 2 S M + K S K T )ϑ 2 } dx+ K S K T Ω B∇ x ϑ • ∇ x ϑ -C f K η(ϑ) ∇ x p • ∇ x ϑϑ + τ c K T τ T K S K η(ϑ) ∇ x p • ∇ x p dx = 0. ( 152
) Since p 2 M - 2K S pϑ M + ( K 2 S M + K S K T )ϑ 2 = (p -K S ϑ) 2 M + K S K T ϑ 2 ,
one can write (152) as 1 2

d dt Ω W (e x (u), p, ϑ) dx+ K S K T Ω B∇ x ϑ • ∇ x ϑ -C f ϑ K η(ϑ) ∇ x p • ∇ x ϑ + τ c K T τ T K S K η(ϑ) ∇ x p • ∇ x p dx = 0. ( 153
)
It remains to study the influence of the terms, containing the gradients. Observe that these terms are a quadratic form in ∇ x ϑ and ∇ x p. Manipulating this quadratic form, one observes that a sufficient condition for its nonnegativity is the absence of real roots for the polynomial

ξ → ξ 2 -C f ||BK -1 ϑ|| 1/2 ∞ ξ η min + τ c K T τ T K S η max
. Here η min and η max are, respectively, the minimal and maximal values of η(ϑ). Since τ c K T τ T K S = ΛT 0 u c ρcL(∆T ) 2 , the nonnegativity holds for ∆T sufficiently small.Text 2

Remark 9. Note that the ratio K S /K T is fundamental in both the definition of the free energy and in the control of the dissipation.

In [START_REF] Coussy | From mixture theory to Biot's approach for porous media[END_REF], a free energy is introduced which is a function of the variables e, ζ and ϑ. Rewriting (148) in terms of these variables gives

2W (e, ζ, ϑ) = Ge : e + M (ζ -(ΦI -B H ) : e) 2 + K S K T ϑ 2 .
This expression reduces exactly to Biot's free energy [START_REF] Biot | Mechanics of Deformation and Acoustic Propagation in Porous Media[END_REF] in the isothermal case, where ϑ = 0, and it has qualitatively the same form as the reduced potential in [START_REF] Coussy | From mixture theory to Biot's approach for porous media[END_REF], page 4634.

σ = Ge x (u) + M (ΦI -B H )(ΦI -B H ) : e x (u) -M (ΦI -B H )ζ -K S ϑI (154)
and

p = M ζ + K S ϑ -(ΦI -B H ) : e x (u), (155)
we are able to make a comparison with the consitutive laws ( 106)-( 107) from [START_REF] Coussy | From mixture theory to Biot's approach for porous media[END_REF] (page 4634). We find that both σ and p contain similar terms, of which now the matricial coefficients are computable in terms of the microstructure.

Comparing with [START_REF] Lee | Thermal consolidation in porous media by homogenization theory-I. Derivation of macroscale equations[END_REF], the constitutive law for the total stress σ is in agreement with equation (144) in its elastic and pressure parts, but not in the temperature part. In [START_REF] Lee | Thermal consolidation in porous media by homogenization theory-I. Derivation of macroscale equations[END_REF] an additional auxiliary cell problem for the thermal coefficient was introduced. It turns out that the result is identical to our matricial coefficient for the temperature in (144). Concerning the constitutive equation ( 146) for the volume content ζ, we observe that in [START_REF] Lee | Thermal consolidation in porous media by homogenization theory-I. Derivation of macroscale equations[END_REF] (equation ( 126), page 135) the expression contains in addition a term involving the gradient of porosity. Furthermore, the temperature equation in [START_REF] Lee | Thermal consolidation in porous media by homogenization theory-I. Derivation of macroscale equations[END_REF] (equation ( 143), page 137) has the transport velocity that contains the effective thermal diffusivity.

For future reference we present the upscaled model in dimensional form. Then it reads

∂ ∂t P ΛM - 4 -2ν (1 + ν)(1 -2ν) β s T M + div x ΦI -B H U + div x V D = 0, (156)
V D = K 2 η(T ) ρ 0 f β f g(T -T 0 - 1 β f )e 3 -∇ x P , (157) ρc
+ β 2 s ΛT 0 M ∂T ∂t + Λβ s T 0 ∂ ∂t - P ΛM -(Φ -1)I -B H : e x (U) - div x B term ∇ x T ) + ρ 0 f c f V D • ∇ x T = h f Φ + h s (1 -Φ), (158)
-div x GΛe x (U) -ΦI -B H P -(1 -Φ)I + B H (3λ + 2ν)β s (T -T 0 ) = ρ 0 f β f Φg(T -T 0 )e 3 -(ρ 0 f Φ + ρ s (1 -Φ))ge 3 , (159)
where

B therm = ρc ρ 0 f c f κ f PeB (the homogenized thermal dispersion coefficient).
All coefficients and variables appearing in (156)-( 159) are defined in Table 3. We conclude this section by two generalizations of the dimensionless upscaled model ( 141)-(147). We make them the content of the following two remarks.

Remark 11. (Weak compressibility of the fluid)

Allowing for weak fluid compressibility, the mass balance reads, see also [[24], eq. (2.

3)],

α ∂p ∂t + div x v f = 0,
where α denotes the fluid compressibility which depends on the fluid bulk modulus and density.

K B = C 0 K S = K S Kg ρ 0 f β f ∆T K S ρs buoyancy O(K S ) ΛcM > 0 compressibility P a
Table 3: Description of parameters and unknowns.

In dimensionless form we find

κ c ∂p ε ∂t + div(w ε + τ T τ c ∂u ε ∂t ) = 0, (160)
where κ c = αp c . This equation replaces (63) in Ω ε f and in the expansions. With minor modifications we obtain the same upscaled equations, except equation (145) which now takes the form

∂ t (Φκ c p + ζ) + τ c τ T div x v D = 0. (161)
Remark 12. (Nonlinear heat capacity)

The upscaled model ( 141)-(147) was derived under the assumption that the temperature T does not deviate much from the reference temperature T 0 . Furthermore we introduced constant heat capacities for fluid and solid. We recall that the heat capacity is defined as T ∂S ∂T . This quantity proves to be virtually independent of the temperature, see for instance [START_REF] Vaxasoftware | [END_REF].

However in our derivation we use the approximation T 0 ∂S ∂T . Hence for large deviation from T 0 , we need to correct the heat capacity by a factor T 0 T . To cover the whole range of T -values we introduce the temperature dependence

R(T ) = T 2 0 + δ 2 T 2 + δ 2 , ( 162
)
where δ is small parameter. Note that R(T ) ≈ 1 for T close to T 0 and R(T ) ≈ 1 T for large T . In the dimensionless setting this leads to the heat capacities

C i (ϑ) = C i R(ϑ) for i = f, s, (163)
where

R(ϑ) = 1 + δ 2 T 2 0 1 + δ 2 T 2 0 + 2∆T ϑ T 0 + ( ∆T ϑ T 0 ) 2 , 1 δ > 0, T 0 > 0, ∆T > 0. ( 164
)
Note that also in the presence of nonlinear heat capacities the Boussinesq approximation remains valid, see e.g. [START_REF] Farina | On the equations governing the flow of mechanically incompressible, but thermally expansible, viscous fluids[END_REF]. Again, with minor modifications we obtain the same upscaled equations, except equation (147) which now reads

K T ∂ ∂t -ζ + div x u + R(ϑ) ∂ϑ ∂t + C f R(ϑ)v D • ∇ x ϑ -div x B∇ x ϑ = hf Φ + hs (1 -Φ). ( 165
)
5 A pore network model: calculation of the homogenized coefficients

In this section, we compute the upscaled matricial coefficients that appear in the equations ( 142)-(147). Specifically, these are

• the thermal diffusion matrix B (111),

• the Gassman tensor G (116),

• the Biot matrix B H (117),

• the coefficient 1/M (132),

• and the permeability tensor K (129).

They are obtained by solving the corresponding auxiliary cell problems

• Cell problem 1: temperature (101)-(104),

• Cell problem 2: solid in terms of elasticity for a ij (114),

• Cell problem 3: solid in terms elasticity for b 0 (115),

• Cell problem 4: fluid in terms of Stokes (127)-(128).

and for the non-isotropic case, we obtain:

B =  
0.950041 6.32571e -06 6.32571e -06 6.32571e -06 0.982624 6.11721e -06 6.32571e -06 6.11721e -06 0.982624   Note that the matrices are (almost) diagonal, as expected from the geometries. The geometrical non-isotropy has only a minor effect, as can be seen from the diagonal terms in the matrix. This is due to the relatively small differences in porosity and the limited contrast D s /D f = 10.

Clearly, the effect of non-isotropy could be enhanced by taking other geometries.

Cell problem 2: Gassmann tensor The Gassmann's tensor is written in Voigt's representation as a 6 × 6 symmetric matrix, involving twenty-one independent coefficients. For the isotropic case, we obtain: As for the thermal diffusivity, the matrices are almost diagonal and again the effect of nonisotropy is minor. Note that the matrix ΦI -B H ≤ 0 in both the isotropic and non-isotropic case.

G =         1
Cell problem 4: Permeability: Solving the Stokes problem, we obtain for the isotropic case the following permeabilities: Note that for the permeability the effect of anisotropy is well noticeable (as expected).

K =   1
We finally remark that we have developed our codes for sequential programming and uniform grid refinement. Therefore, in these 3D simulations, the discretization error may be relatively large. For this paper, the obtained simulation values provide a very good understanding and most importantly substantiate our theory. To enhance the precision (i.e., reducing the discretization error) future studies may address local mesh adaptivity and parallel computing in which the last two authors have recent work (e.g., [START_REF] Lee | Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model[END_REF]) on very related problems in fractured porous media.

Conclusion

In this paper we derive the thermal Biot equations by the method of two-scale expansions. Starting at the pore level with the momentum and the energy balance equations in the fluid and in the solid phase, together with the appropriate conditions at the fluid-solid interface. We derive an upscaled formulation. It is given in terms of Darcy's law, and the balances for momentum, mass and energy. We also construct a Lyapunov functional for the upscaled system, which implies stable behavior of the corresponding solutions. We also provide details when comparing our results to those obtained by Coussy et al [START_REF] Coussy | From mixture theory to Biot's approach for porous media[END_REF] and Lee and Mei [START_REF] Lee | Thermal consolidation in porous media by homogenization theory-I. Derivation of macroscale equations[END_REF].

Finally we consider, as an example, a periodic porous structure in the form of a pore network model. We evaluate numerically the matricial coefficients in the upscaled equations by solving the corresponding auxiliary problems. Two specific cases are considered: an isotropic structure and a non-isotropic structure, in which there is a preference in the single direction. The numerical results clearly demonstrate the effect of non-isotropy of the porous structure on the coefficients of the upscaled model. This is part 1 of a series of papers in preparation on the subject. Future work focuses on rigorous mathematical analysis of the equations, on the analysis of benchmark problems and on numerical simulations of concrete nuclear waste and CCS examples.

  connected domain, filled with an incompressible fluid to which Stokes' equations apply.

Figure 1 :

 1 Figure 1: (left) Unit cell of isotropic pore network model, (right) unit cell of anisotropic pore network model. The fluid part is in blue. In Sections 3 -5, this geometry corresponds to the notation Y = (0, 1) 3 (unit cell), where Y f is the fluid part and Y s is the solid part.

  QUANTITY VALUE FOR ROCKS VALUE FOR CLAYS Young's modulus Λ = 10 10 Pa Λ = 10 7 Pa Solid thermal expansion coefficient βs = 0.5 • 10 -5 / 0 K βs = 10 -6 / 0 K Fluid thermal expansion coefficient β f = 10 -4 / 0 K β f = 10 -4 / 0 K Solid thermal diffusivity Ds = 10 -6 m 2 /s Ds = 10 -6 m 2 /s Fluid thermal diffusivity D f = 10 -7 m 2 / s D f = 10 -7 m 2 / s Rock mass density ρs = 2.5 • 10 3 kg/m 3 ρs = 2.5 • 10 3 kg/m 3 Fluid mass density ρ 0 f = 10 3 kg/m 3 ρ 0 f = 10 3 kg/m 3 Water viscosity ηc = 10 -3 Pa s ηc = 10 -3 Pa s Macroscopic length L = 1 m L = 1 m Pore scale size = 10 -7 m = 10 -7 m Small parameter ε = /L = 10 -7 ε = /L = 10 -7 Average heat capacity ρc = 10 6 J/ (m 3 0 K) ρc = 10 6 J/ (m 3 0 K)

  For rock type material we choose a characteristic pressure p c = 10 7 Pa, giving (from (46)) for u c = O(10 -3 m). For clay type material we choose p c = 10 5 Pa, yielding u c = O(10 -2 m). With these values we have a clear separation of scales in the sense that u c

  ), p 0 does not depend on y on S and S p 0 n dS y = 0. Therefore, by the Fredholm alternative 1 , problem (89) has a solution. The solution is unique if we impose Ys u 1 dy = 0. STEP 4

Ys div y b 0

 0 dy = Ys Ae y (b 0 ) : e y (b 0 ) dy > 0 (132) and evaluating Ys div y u 1 (x, y, t) dy = 3 i,j=1

ΩW

  (e x (u), p, ϑ)(t) dx is monotone decreasing in time.Sketch of the proof. For simplicity we suppose homogeneous Dirichlet's conditions for u and homogeneous Neumann conditions for p and ϑ. First we multiply equation (143) by ∂ t u and integrate over Ω. This yieldsΩ Ge x (u) : e x (∂ t u) dx -Ω ΦI -B H : e x (∂ t u)p dx-Ω (1 -Φ)I + B H : e x (∂ t u)K S ϑ dx = 0. (149)Next we test equation (145) by p and obtain Ω p -K S ϑ M ∂ t p dx + Ω ΦI -B H : e x (∂ t u)p dx -Ω v D • ∇ x p dx = 0. (150)

  Note that ∂ ∂e (W (e, ζ, ϑ) + K S ϑ(ζ -I : e)) = σ and ∂ ∂ζ (W (e, ζ, ϑ) + K S ϑ(ζ -I : e)) = p. Remark 10. (Comparison to literature) Recasting expressions (144) and (146) into the form

  .24019e -05 -5.26e -22 -2.05e -23 -1.14e -21 1.24019e -05 -3.21e -21 3.18e -23 -4.14e -21 1.24019e -05   and for the non-isotropic case:

  17 -2.97e -17 -1.36e -17 -4.69e -18 -1.94e -16 1.33e -17 0.408828 -1.08e -16 6.85e -19 -3.78e -18 -7.45e -18 -2.97e -17 -1.08e -16and for the non-isotropic case (an orthotropic effective material), we calculate: Cell problem 3: Biot matrix and coefficient 1/M : Here, we compute Biot's matrix B H and the coefficient 1/M again for the two specific geometries. The calculation for the isotropic case, yields:

			.21187	0.428817	0.428817	-4.69e -18 -3.78e -18	7.0e -17	
			0.428817 0.428817	1.21187 0.428817	0.428817 1.21187	-1.94e -16 -7.45e -18 -9.92e -18 1.33e -0.408828 8.57e -17	      
			7.0e -17	-9.92e -18 -1.36e -17 6.85e -19	8.572e -17	0.408828
			1.15945		0.434037	0.434037	1.03e -17	3.5e -18	-3.83e -17	
	G =	       3.50433e -18 -1.11e -17 -1.31e -16 -1.3e -17 0.434037 1.2435 0.45122 -5.49e -16 -1.11e -17 -7.62e -18 0.434037 0.45122 1.2435 1.68e -17 -1.31e -16 1.13e -17 1.03e -17 -5.49e -16 1.68e -17 0.439163 -1.3e -17 -3.44e -18 0.426031 3.49e -17	      
			-3.83e -17 -7.62e -18 1.13e -17 -3.44e -18 3.49e -17	0.426031
							0.116782		2.93853e -17	1.12859e -16	
			B H =	 2.93853e -17	0.116782	-7.11201e -19	
					1.12859e -16 -7.11201e -19	0.116782
	and							1 M	= 0.182911.
	The computations for the non-isotropic case yield:
							0.125088	5.52671e -18 4.0376e -17	
				B H =	 5.52671e -18	0.175705	7.36425e -17 
							4.0376e -17 7.36425e -17	0.175705
	and							1 M	= 0.255538.

Comment for the non-mathematical reader: Fredholm's alternative gives a necessary and sufficient criteria for solvability of an equation, in the critical situation when we are in a spectrum. For linear algebraic system Ax = b, it says that if 0 is an eigenvalue of the matrix A, the system has a solution if and only if b is orthogonal to the eigenvectors of A that correspond to the eigenvalue 0. In all examples considered here 0 is an eigenvalue, with the eigenfunction being either a constant or a constant vector. Therefore the corresponding boundary value problem in the y-variable admits a solution if and only if the volume average of the right hand side with respect with respect to y is equal to the interface average of the flux. We refer to the textbooks as[START_REF] Wloka | Partial differential equations[END_REF], for the Fredholm theory for the partial differential equations.

* C.J.v.D. acknowledges support of the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR), as well as the support of the Center for Subsurface Modeling, Institute for Computational Engineering and Science (ICES), UT Austin. A.M. was partially supported by Darcy Center, The Netherlands and by The J. Tinsley Oden Faculty Fellowship Research Program at the Institute for Computational Engineering and Science (ICES), UT Austin. T.W. was partially supported by the Center for Subsurface Modeling and the Institute for Computational Engineering and Science (ICES) through the J. Tinsley Oden Faculty Fellowship Research Program.

Discretization

The discretization of the four cell problems is based on a Galerkin finite element scheme. As usually done, the strong forms are integrated over the domain Y = (0, 1) 3 to obtain variational formulations. The grid consists of hexahedral cells. We introduce H 1 conforming discrete spaces (see e.g., [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]), which consist of continuous trilinear functions Q c 1 for the thermal problem, and correspondingly vector-valued spaces for the two elasticity cell problems. The Stokes equations are discretized with the Taylor-Hood element in which triquadratic Q c 2 functions are adopted for the velocity and Q c 1 functions for the pressure, yielding an inf-sup stable finite element pair. The (linear) discrete systems can be solved with a conjugate gradient method for the thermal and elasticity cell problems and a Schur complement technique for the Stokes system. The basic programming code is based on the deal.II finite element library ([4], [START_REF] Bangerth | II -a General Purpose Object Oriented Finite Element Library[END_REF]) and uses parts of deal.II-step-22, deal.II-step-45, and the fluid-structure program [START_REF] Wick | Solving monolithic fluid-structure interaction problems in arbitrary Lagrangian Eulerian coordinates with the deal. II library[END_REF].

Numerical simulations

We now provide more details on the geometry, boundary conditions and material parameters and finally discuss our numerical results.

Geometries. We consider two specific porous media that have the form of the pore network models from Figure 1. Figure 1 (left) shows the base case, in which the geometry is isotropic and Figure 1 (right) shows a non-isotropic case, where the cross section of the channels in (y, z)-directions is different from the cross section in the x-direction. With respect to the base case we take the cross section of the channels in (y, z)-directions the same and in x-direction a factor 2 × 2 larger. In the computations we impose 1-periodicity of solutions to the cell problems.

The dimension of the fluid part in the unit cell is such that Φ = |Y f | = 0.0820312 in the isotropic case; 0.111328 in the non-isotropic case.

Material parameters. We suppose that the elastic material of the solid part is isotropic, with the elasticity tensor A given by [START_REF] Gallavotti | Foundations of Fluid Dynamics[END_REF]. Since the coefficients are divided by the Young's modulus Λ, we just need to set the Poisson ratio. Here we take ν = 0.25. This implies that the dimensionless Lamé first parameter λ and the shear modulus µ are both equal to 0.4 and that the corresponding dimensionless bulk modulus is K = 2/3. Since the ratio of D s /D f = 10 for both clay and rock, we run only one set of auxiliary problems where D f = 0.1 and D s = 1. The actual results for clay and rock are obtained by a trivial scaling.

Discussion of our findings. We discuss our findings in the previous order: temperature, elasticity, and Stokes.

Cell problem 1: Thermal diffusivity: The results are given in terms of B = B/D s (by scaling).

For the two specific geometries considered here, we calculate for the isotropic case the following values of the thermal diffusion matrix:


Sipho Ngobeni

Hein Venter

Ivan Burke

A FORENSIC READINESS MODEL FOR WIRELESS NETWORKS

Keywords:

Over the past decade, wireless mobile communications technology based on IEEE 802.11 wireless local area networks (WLANs) has been adopted worldwide on a massive scale. However, as the number of wireless users has soared, so has the possibility of cyber crime, where criminals deliberately and actively break into WLANs with the intent to cause harm or access sensitive information. WLAN digital forensics is seen not only as a response to cyber crime in wireless environments, but also as a means to stem the increase of cyber crime in WLANs. The challenge in WLAN digital forensics is to intercept and preserve all the communications generated by the mobile devices and conduct a proper digital forensic investigation. This paper attempts to address this issue by proposing a wireless forensic readiness model designed to help monitor, log and preserve wireless network traffic for digital forensic investigations. A prototype implementation of the wireless forensic readiness model is presented as a proof of concept.

Introduction

Wireless technologies have become very popular around the world. Wireless local area networks (WLANs) or "hotspots" blanket public places such as convention centers, airports, schools, hospitals, railway stations, coffee shops and other locations to provide seamless public access to the Internet [START_REF] Velasco | Wireless forensics: A new radio frequency based location system[END_REF]. These hotspots provide several advantages over hard-wired networks, including user mobility and flexible Internet access. However, due to their open nature, WLANs have become a major target for cyber criminals.

WLAN digital forensics involves the application of methodologies and tools to intercept and analyze wireless network events for presentation as digital evidence in a court of law [START_REF] Siles | Wireless forensics: Tapping the air -Part one[END_REF]. As such, WLAN digital forensics is complementary to intrusion prevention -when intrusion prevention fails, WLAN digital forensics is useful for obtaining information about the intrusion. However, the primary challenge in WLAN digital forensics is to acquire all the digital evidence related to a crime [START_REF] Newman | Computer Forensics: Evidence Collection and Management[END_REF]. This challenge arises from the fact that the devices participating in a WLAN environment are mobile. Furthermore, since the devices are not always connected to the network, it is difficult to attribute criminal activity to a particular device.

This paper proposes a wireless forensic readiness model for monitoring, logging and preserving wireless network traffic for digital forensic investigations. The wireless forensic readiness model builds on the work of Rowlingson [START_REF] Rowlingson | A ten step process for forensic readiness[END_REF] related to traditional forensic investigations. A prototype implementation of the readiness model is presented as a proof of concept.

Wireless Local Area Networks

WLANs represent a ubiquitous technology that provides seamless high-speed Internet connectivity at public locations. Unlike traditional LANs, WLANs connect computers to the network without physical (wired) connections. WLANs offer tremendous user mobility, enabling users to access files, network resources and the Internet [START_REF] Scarfone | Guide to Securing Legacy IEEE 802.11 Wireless Networks[END_REF].

Criminal Misuse of WLANs

The lack of a physical connection between a WLAN and its participating mobile devices causes crimes to remain discreet, especially since the mobile devices are potentially far removed. This fact needs to be considered when digital evidence is identified and collected in an investigation involving wireless traffic. Potential criminal misuse of WLANs include [START_REF] Turnbull | The 802.11 technology gap -Case studies in crime[END_REF][START_REF] Turnbull | Wi-Fi network signals as a source of digital evidence: Wireless network forensics[END_REF]: WLAN Detection and Connection: This type of misuse involves an intruder using the wireless medium as a tool to commit other criminal activities (e.g., unauthorized use of the WLAN or use of the WLAN as a launch pad for other criminal activities).

Concealment of Digital Evidence:

This type of misuse involves hidden wireless devices or hidden wireless networks (e.g., fake access points).

WLAN as an Attack Vector: This type of misuse involves attacks against the networked (mobile) devices originating from the wireless network, and attacks against the WLAN medium itself.

Sources of Digital Evidence

WLANs typically incorporate 802.11-based wireless devices. The locations where digital evidence is stored on devices and the extraction of evidence are dependent on the specific wireless device. However, the fundamental problem with 802.11-based wireless devices is their lack of a physical footprint, which is the most crucial issue in the identification of these devices [START_REF] Turnbull | Wi-Fi network signals as a source of digital evidence: Wireless network forensics[END_REF].

It is imperative to locate all the relevant wireless devices in a digital forensic investigation. Various open source and commercial tools (e.g., Wireshark, Kismet and AirCapture) may be used by a digital forensic investigator to identify wireless networks within range, the devices connected to the wireless networks and, possibly, the locations of the wireless devices [START_REF] Turnbull | Wireless forensic analysis tools for use in the electronic evidence collection process[END_REF]. The principal drawback of these tools is the high packet drop rate [START_REF] Broadway | Improving the analysis of lawfully intercepted network packet data captured for forensic analysis[END_REF]. The large volume of network traffic makes it difficult for the tools to accept and store all the packets; some packets may be dropped, resulting in the loss of evidence.

WLAN Digital Forensics

Digital forensics deals with the investigation of computers and other digital devices believed to be involved in criminal activities [START_REF] Francia | Computer forensics laboratory and tools[END_REF]. WLAN digital forensics involves the application of methodologies and tools to capture and analyze wireless network traffic that can be presented as evidence in a court of law [START_REF] Siles | Wireless forensics: Tapping the air -Part one[END_REF]. A WLAN digital forensic methodology is a digital forensic process; the tools are software systems that intercept and analyze network traffic. A digital forensic process is a procedure that is followed to investigate a particular criminal activity involving digital evidence [START_REF]Handbook of Computer Crime Investigation: Forensic Tools and Technology[END_REF]. Every digital forensic investigation must go through the following phases of the digital forensic process: Define the scope and goals of the investigation.

Determine the work and materials.

Acquire the images of the devices to be examined.

Perform the digital forensic analysis.

Prepare the report.

Currently, EnCase and FTK are the most popular tools used in digital forensic investigations. The phases of the digital forensic process for EnCase are preview, imaging or acquisition, verification, recovery and analysis, restoration and archiving; the phases for FTK are detection, 1 lists the phases for EnCase and FTK along with those for the wireless forensic readiness model (WFRM), which is described in the next section.

According to Table 1, only the analysis phase is common to EnCase, FTK and WFRM. The preservation and analysis phases are common to FTK and WFRM. However, it is worth noting that the digital forensic processes for FTK and EnCase are essentially the same as far as the general digital forensic process is concerned. This suggests that the phases correlate although they are named differently. The inconsistent naming of phases is due to the fact that the digital forensic processes for forensic tools are not standardized. In this paper, we adopt the general digital forensic process described by Casey [START_REF]Handbook of Computer Crime Investigation: Forensic Tools and Technology[END_REF].

Researchers have studied various issues related to wireless network forensics. Yim, et al. [START_REF] Yim | The evidence collection of DoS attack in WLAN by using WLAN forensic profiling system[END_REF] proposed a WLAN forensic profiling system for collecting digital evidence after denial-of-service attacks on WLANs. Turnbull and Slay [START_REF] Turnbull | Wi-Fi network signals as a source of digital evidence: Wireless network forensics[END_REF] consider the potential sources of digital evidence in 802.11-based wireless networking environments. Then [START_REF] Then | Examining wireless access points and associated devices[END_REF] discusses methods for examining wireless access points to determine if the devices of interest are connected or were connected to a wireless network. While these efforts and others are useful, a digital forensic readiness approach for WLANs has yet to be articulated.

Digital Forensic Readiness

The purpose of digital forensic readiness is to reduce the effort involved in performing an investigation while maintaining the level of credibility of the digital evidence being collected [START_REF] Endicott-Popovsky | A theoretical framework for organizational network forensic readiness[END_REF]. The decrease in effort includes reductions in the time and the cost of incident response. An organization that is "forensically ready" can respond to an attack rapidly and efficiently. In general, reducing the time involved in incident response reduces the cost of the investigation.

Tan [START_REF] Tan | Forensic readiness: Strategic thinking on incident response[END_REF] discusses an incident in which an intruder took approximately two hours to launch an attack, but digital forensic experts required 40 billable hours to respond to the incident. The response took such a long time because the organization was not forensically prepared for the incident. Organizations deploying WLANs that are at a high risk of cyber attack should be ready to collect digital evidence before an incident occurs. The model presented in the next section addresses the concept of digital forensic readiness in WLANs.

Wireless Forensic Readiness Model

The most salient characteristic of the wireless forensic readiness model (WFRM) is that it monitors wireless network traffic at access points. The monitored traffic is stored in a log file and the integrity of the stored information is preserved. Thus, the information needed by digital forensic investigators is readily available should the need arise. The availability of the information reduces the cost of conducting the digital forensic investigation because a major portion of the digital forensic process (monitoring, logging and preservation) has already been conducted based on the WFRM.

Figure 1 shows the five phases of the digital forensic process corresponding to the WFRM. As listed in Table 1, the phases are monitoring, logging, preservation, analysis and reporting.

Phase 1 (monitoring) shows several mobile devices (MDs) connected to a WLAN through different access points (APs). The mobile devices use the access points to connect to the Internet. In addition to providing Internet connectivity, the access points are modified (for the purposes of this model) to monitor all the traffic generated by the mobile devices. For security reasons, the monitoring component uses a firewall to filter inbound and outbound wireless traffic. Filtering is the process of controlling access to the WLAN by screening packets based on the content of their headers [START_REF] Velasco | Wireless forensics: A new radio frequency based location system[END_REF].

Phase 2 (logging) records all the traffic monitored by the access points. Each access point has its own capture unit (CU) that logs the traffic passing through it. The log file is divided into separate storage areas, each consisting of (for example) 1 MB of data. When the buffer of a capture unit is full, a fixed-size block of data is moved to permanent storage. For example, B1 in Phase 2 represents a block of data consisting of 4 MB.

In Phase 3 (preservation), the capture unit sends the accumulated blocks of data to the evidence store (ES). The capture unit computes a hash value for each block of data, which is saved in the hash store (HS) for integrity checking purposes. Phase 4 involves the analysis of the stored data and Phase 5 involves the creation of a report.

WFRM Simulation

The WFRM prototype was simulated using AnyLogic Professional (version 6.0) [START_REF] Coensys | AnyLogic 6: Multi-Paradigm Simulation Software[END_REF], a Java-based, multi-paradigm, hybrid simulation tool capable of modeling systems as a combination of discrete events, system dynamics and agents. The simulation is designed to validate the use of the WFRM for implementing digital forensic readiness in WLAN environments.

Figure 2 shows a graphical representation of the WFRM before the simulation starts. The MobileDevice component generates random simu- CaptureUnit contains the variables currentCache, cacheSize and cur-rentCacheUsage. The currentCache variable represents the log file in our model, which works like a buffer; currentCache can store up to eight packets (based on the cacheSize in Figure 3). The eight captured packets are put together to form a single message; this represents a created block of data in our model. CaptureUnit computes the hash value of the formed message and stores the value in the HashStore; also, it passes the formed message to the EvidenceStore for storage. The variable current-CacheUsage keeps track of the number of times that currentCache was filled with the eight packets that are combined to form a single message.

Figure 4 presents sample data in EvidenceStore. The message in Row 2 shows that an anonymous user with IP address 132.143.133.122 is attempting to log into a remote host via FTP. The fact that this machine is anonymous could be of interest to a digital forensic investigator. The message in Row 5 contains data such as "Please send your message," "hacker1@badSMTP.com" and "zombieNET@hostPC.org." This data seems suspicious and constitutes potential digital evidence.

The message in Row 9 shows that a machine is using HTTP to access Google; this does not appear to indicate any malicious activity. The message in Row 13 shows that the machine with IP address X is asking the machine with IP address Y if it knows the machine with IP address Z; this does not appear to be malicious. However, if Machine Y responds to Machine X that it is not aware of Machine Z, then it is possible that Machine Z is not part of the network and could be an intruder who intends to sniff network traffic between Machines X and Y. The message in Row 15 shows that a machine with IP address 132.143.133.169 is accessing a suspicious website named "www.illegalSite.com;" an error message "Do not go here" pops up when this website is accessed. The message in Row 19 shows that a machine is providing its login details to a website and downloading a suspicious file named secrets.bat.

Figure 5 presents the HashStore corresponding to the captured simulated messages. Every time a message is captured, a copy of the original captured message is hashed and transferred to the HashStore. The main reason for hashing the captured information and keeping a separate copy of the original information is to verify the integrity of the captured information and to determine whether or not it was tampered with. The integrity of any message can be verified by extracting and computing the hash value (y) of the message stored in the EvidenceStore. The hash value (y) for the particular message block is then retrieved from the HashStore. If the hash values x and y match, the content of the original captured message was not tampered with and the integrity of the captured message in the EvidenceStore is verified. The integrity checking mechanism was built into the prototype because the integrity of evidence is a crucial requirement in any digital forensic investigation [START_REF]Handbook of Computer Crime Investigation: Forensic Tools and Technology[END_REF].

Discussion

In the simulation described in the preceding section, the simulated packets were logged (by CaptureUnit) and preserved (by EvidenceStore and HashStore). However, note that traffic monitoring was not implemented in the prototype because it is performed by the access point mainly for security reasons. After the traffic generated by the mobile devices that have connected to the WLAN has been captured and preserved, the data is ready for analysis in a digital forensic investigation. Because this data is forensically ready and forensically sound, the time and cost involved in conducting the digital forensic investigation are reduced considerably. In fact, the data needed for the investigation is readily available and the bulk of the digital forensic process (i.e., monitoring, logging and preservation) has been completed.

One disadvantage of the WFRM simulation is that the traffic is preserved in the EvidenceStore and HashStore, which potentially requires a large amount of storage space. This is not a serious problem be-cause storage is becoming ever cheaper. Nevertheless, we are working on compression techniques that will facilitate the preservation of the entire stream of wireless network traffic. Since this might not be an optimal long-term solution to the problem, further research is needed to address the storage issue.

Finally, we note that the digital forensic processes for EnCase, FTK and WFRM (Table 1) are essentially equivalent. However, since our emphasis in this paper is the design of a readiness model, the practical implementation of the digital forensic process employed for WFRM is different from the conventional digital forensic process models for En-Case and FTK.

Conclusions

The wireless forensic readiness model helps address the twin challenges of intercepting and preserving all the communications generated by mobile devices in WLANs. In general, WLANs are not forensically prepared to gather digital evidence for use in ensuing investigations. The forensic readiness model focuses on the monitoring, logging and preservation of wireless network traffic. This covers the bulk of the general digital forensic investigation process, reducing both the time and the cost of forensic investigations.

Our future research will focus on several issues. One issue, as mentioned above, is the efficient storage of data in the hash store and evidence store. Another key issue is the analysis of potentially large amounts of data gathered as a result of the application of the wireless forensic readiness model. Other issues involve evidence management and the consideration of infrastructure requirements, admissibility requirements and retention requirements.

Figure 1 .

 1 Figure 1. Wireless forensic readiness model (WFRM).

Figure 2 .

 2 Figure 2. WFRM before the simulation.

Figure 3 .

 3 Figure 3. WFRM during the simulation.

Figure 3

 3 Figure3shows the simulated messages flowing through the network during the simulation. The simulated messages correspond to the packets that pass through the network from various devices in the WLAN.CaptureUnit contains the variables currentCache, cacheSize and cur-rentCacheUsage. The currentCache variable represents the log file in our model, which works like a buffer; currentCache can store up to eight packets (based on the cacheSize in Figure3). The eight captured packets are put together to form a single message; this represents a created block of data in our model. CaptureUnit computes the hash value of the formed message and stores the value in the HashStore; also, it passes the

Figure 4 .

 4 Figure 4. Evidence store.

Figure 5 .

 5 Figure 5. Hash store.

Table 1 .

 1 Digital forensic phases for EnCase, FTK and WFRM.

	EnCase	FTK	WFRM
	1. Preview	1. Detection	1. Monitoring
	2. Imaging	2. Identification 2. Logging
	3. Verification	3. Analysis	3. Preservation
	4. Recovery and Analysis 4. Preservation	4. Analysis
	5. Restoration	5. Reporting	5. Reporting
	6. Archiving
	identification, analysis, preservation and reporting. Table