Identification of viscoelastic properties from numerical model reduction of pressure diffusion in fluid-saturated porous rock with fractures

This paper deals with the computational homogenization and numerical model reduction of deformation driven pressure diffusion in fractured porous rock. Exposed to seismic waves, the heterogeneity of the material leads to local fluid pressure gradients which are equilibrated via pressure diffusion. However, a macroscopic observer is not able to measure the diffusion process directly but senses the intrinsic attenuation of an apparently monophasic viscoelastic solid. The aim of this paper is to establish a reliable, yet numerically efficient, computational homogenization method to identify the viscoelastic properties of the macroscopic substitute model. Inspired by the Nonuniform Transformation Field Analysis, we incorporate a Numerical Model Reduction procedure. The proposed method is validated for several scenarios ranging from pressure diffusion in an unfractured poroelastic matrix, via localized pressure diffusion in interconnected fractures embedded in an impermeable matrix, to the fully coupled pressure diffusion both in fractures and the embedding poroelastic matrix.


Introduction
Fractures that are present in subsurface formations dominate not only their mechanical properties but also their hydraulic properties [2,31]. When interconnected, they can increase the hydraulic conductivity of a rock by several orders of magnitude [20,28]. This makes the detection and characterization of fracture networks of enormous interest for applications such as the geological sequestration of CO 2 , production of deep geothermal energy, nuclear waste storage, and the explo- Institute of Earth Sciences, University of Lausanne, 1015 Lausanne, Switzerland ration and production of oil and gas. On the one hand, this interest is associated to the necessity of extremely low permeabilities over time in nuclear waste storage sites as well as in cap rocks on top of CO 2 injection sites. On the other hand, high permeability zones are the targets in geothermal and in non-conventional hydrocarbon reservoirs.

Pressure diffusion in fluid-saturated fractured rock
Squirt flow in interconnected micro-cracks causes dissipation of acoustic waves due to viscous friction [19]. This dissipation mechanism is governed by fluid pressure diffusion. An identical phenomenon occurs at a larger spatial scale in hydraulically interconnected mesoscopic fractures [26], here referred to as squirt-type flow. Two interconnected fractures are differently deformed by the propagating wave mainly as a consequence of their different orientation. Hydromechanical coupling induces different fluid pressures in the fractures and, thus, fluid pressure diffusion occurs from one fracture into the other one, causing dissipation. The resulting dissipation is frequency dependent and associated with dispersion of the elastic moduli. These effects of fracture interconnectivity on seismic waves represent a great potential for inferring hydraulic characteristics of subsurface formations through the interpretation of seismic data. A few numerical upscaling approaches have been used in the recent years to compute, at the level of a Representative Volume Element (RVE), frequency-dependent attenuation and the corresponding dispersion of the elastic moduli caused by squirt-type flow in interconnected fractures [22,23,26,30]. They describe the squirt-type flow within the fractures either by using Biot's quasi-static equations of poroelasticity or quasi-static, linearized Navier-Stokes equations, whilst the embedding rock matrix is described either as a poroelastic medium or as a monophasic linear-elastic medium based on Hooke's law.
Numerical simulations of wave propagation can be used to predict the effect of fracture interconnectivity on surface seismic data, when the subsurface geometry and properties are known. However, in this case, model domains much larger than the RVE considered for the previously mentioned numerical upscaling approaches are necessary which, if feasible at all, requires tremendous computational efforts. A solution for this problem is to replace the heterogeneous medium at the macro-scale by an equivalent homogeneous viscoelastic medium [4] exhibiting identical frequency-dependent attenuation and stiffness modulus dispersion. This presumes that the fluid pressure diffusion observed in the heterogeneous medium at the mesoscopic scale is not observed at the macroscopic scale [13]. In [15], the attenuation-dispersion behavior of a heterogeneous poroelastic medium, consisting of a two-dimensional periodic array of circular inclusions, was roughly approximated with a Maxwell-Zener model containing one single Maxwell chain. This procedure matched the maxima of the corresponding frequency-dependent P-and S-wave attenuation curves but failed to predict the high-frequency limits of the stiffness moduli. Hence, it is an open challenge to find a more fundamental approach for the identification of a suitable viscoelastic model that is able to provide a precise match of attenuation-dispersion behavior. Additionally, such an approach should be able to handle more realistic description of heterogeneities and stochastic fracture networks on the meso-level.

Computational homogenization and numerical model reduction
In this paper we develop a Numerical Model Reduction (NMR) technique based on the concepts of computational homogenization. We, therefore, refer to the length scales involved in the problem as shown in Fig. 1. We distinguish between, firstly, the macro-scale with the characteristic length L where propagation and attenuation of seismic waves is observed, secondly, the meso-scale with the characteristic length l L where the squirt-type flow in fracture networks as well as pressure diffusion in the embedding porous rock takes place and, thirdly, the micro-scale with the characteristic length λ l as the scale of discrete grains and pores. In order to connect the micro-and the meso-scale we employ the Theory of Porous Media (TPM) [6,7] to model the fluid-saturated porous rock as a biphasic poroelastic medium with possibly heterogeneous material properties. In fact, the resulting mesoscopic description of the rock matrix corresponds to Biot's equations of quasi-static linear consolidation [3]. In addition, the mesoscopic model is enriched by a fracture network, where we treat the fractures as fluid-filled mechanically and hydraulically open sharp interfaces embedded in the poroelastic matrix. Hence, fluid pressure diffusion occurs in the poroelastic matrix as well as along the fractures. We denote the resulting model that couples fluid pressure diffusion on those two topology levels as "hybrid-dimensional interface model" [30].
It is important to remark that quasi-static considerations are admissible and reasonable on the SVE-level due to separation of scales. Hence, the macroscopic wavelength at seismic frequencies is much larger than the SVE size. Moreover, the inertia forces at those low frequencies remain, at the SVE-level, much smaller than the viscous drag forces in the SVE, see [24].
The primal interest in this contribution, however, is the scale transition from the heterogeneous meso-level towards the homogenized macro-scale. Here, we promote the concept of computational homogenization based on volume averaging techniques. This concept, often referred to as FE 2 approach as coined in [8], is well established in the literature for various applications. A comprehensive overview of the concept can be found for example in [17,32]. Applying the FE 2 approach to pressure diffusion in fractured rock requires that we define a mesoscopic volume element considered to be representative for the entire meso-structure. This Representative Volume Element (RVE) represents one single material point on the macro-scale. It is important to remark that, for practical applications, the volume element that can be used for simulations with reasonable numerical efforts is usually much smaller than a true RVE. We, therefore, introduce the concept of Statistical Volume Elements (SVE) as proposed for example in [21]. In order to correlate macroscopic and mesoscopic quantities we execute volume averaging of selected mesoscopic quantities to compute their macroscopic counterparts. However, we assume that squirttype flow and pressure diffusion in the embedding matrix induced by seismic waves are local phenomena. In other words, fluid pressure diffusion takes place inside the SVE, but any macroscopic transport is neglected. In this case the macroscopic substitute material behaves like a monophasic viscoelastic solid. This phenomenon was discussed in detail in [13] for heterogeneous poroelastic media. The influence of the chosen SVE size for pressure diffusion in fractured porous media was studied in [5] similarly to the classical homogenization problem investigated in [14].
However, the FE 2 approach requires the nested solution of one macroscopic and numerous mesoscopic initial boundary value problems, which leads to tremendous computational costs thus restricting practical applicability to rather coarse meso-scale problems. In this paper, we therefore propose a NMR strategy that significantly reduces the numerical efforts in solving the SVE problems. Our approach is inspired by the Nonuniform Transformation Field Analysis (NTFA), which was initially established for elasto-viscoplastic materials [18,25] and recently extended towards solids with cohesive interfaces [10], generalized standard media [9,27], poroelas-tic composites [12] and transient heat flow [1], to name only a few application fields. The advancement presented in this paper is that we significantly extend the NMR-methodology that was initially proposed for undamaged poroelastic media in [12] towards pressure diffusion in a hybrid-dimensional formulation.
The key ingredient of our scheme is that we approximate the fluid pressure in the hybrid-dimensional domain as a linear combination of so-called pressure modes which we identify via a Proper Orthogonal Decomposition (POD) from training computations on the SVE level. The pressure modes span a reduced basis for the spatial approximation of the mesoscopic pressure diffusion problem and, as a consequence, allow for the identification of the viscoelastic evolution equations of the macroscopic substitute model. Hence, the reduced FE 2 scheme is split into two stages: The numerically expensive training computations on SVE level define the first stage. They need to be executed only once in advance to derive the homogenized substitute model. Therefore, the first stage is often referred to as "offline" phase. The macroscopic substitute model is used in the second stage, called "online" phase, to solve the macroscopic initial boundary value problems. This allows us to access the full sub-scale information during the "online" phase without explicitly solving the underlying SVE problems and, thus, we succeed to reduce the computational costs tremendously.
The paper is organized as follows: In Sect. 2 we describe the hybrid-dimensional modeling approach on the SVElevel. We, therefore, employ a sharp interface model for pressure diffusion along fractures embedded in a poroelastic rock matrix. Section 3 is devoted to the computational homogenization framework and the definition of the mesoscopic initial boundary value problem for the hybrid-dimensional description in a FE 2 sense. In Sect. 4 we develop the novel NMR technique which allows us to derive viscoelastic substitute models for arbitrary diffusion problems. Finally, we validate our methodology in several numerical experiments presented in Sect. 5.
Throughout the manuscript, vector and tensor quantities are written as bold types. Taking into account Einstein's sum convention we write out simple and double contractions as a · b = a i b i and A : B = A i j B i j . In Voigt notation, vector and matrix quantities are written as underlined italic types, for example a and A.

Hybrid-dimensional interface model for fractured poroelastic media
We introduce a cubic SVE in the mesoscopic domain Ω with the volume V := |Ω | = l 3 , where l is the edge length l of the SVE. The SVE contains n fluid-filled thin frac- The internal surface separating the matrix material from the fracture space is denoted ∂Ω i M . We define the surface normal vector n such that it points from ∂Ω i M into Ω F . For seismic attenuation and pressure diffusion in fracture networks the case of very thin fractures is of major interest. In fact this means that we consider the case that V F V M . Since it would be numerically extremely costly to discretize the volume of such thin fractures we aim to simplify the geometrical description of the fracture space. We, therefore, condense the fracture volume towards a set of n possibly intersecting planar interfaces ∂ F k , k = 1, 2, . . . , n, see Fig. 2b, whilst ∂ F := ∪ n k=1 ∂ F k . In the fracture space, we assume all quantities to be homogeneous in thickness direction. In other words, we neglect any gradients perpendicular to ∂ F k . This hybrid-dimensional interface model is discussed in detail in [29].
Hence, we can state that a volume element dv in the fracture Ω F,k computes as dv = τ da. Here, we use the interface element da of ∂ F k and introduce the (current) fracture open- 1 Note that this assumption is invoked here for the sake of simplicity. Existence of a symmetry plane is not ultimately necessary to set up the interface model.
In order to distinguish between positive and negative fracture surfaces ∂Ω i,± M we assign the fracture normal vector n F perpendicular to the interface ∂ F k . By convention, n F shows in the direction of the positive fracture surface ∂Ω i,+ M,k and, therefore, n + := −n F and n − := n F .
The fracture opening is defined as the apparent normal distance of the fracture surfaces and can be computed, for all k = 1, 2, . . . , n, as Here, we evaluate the solid displacement u of the poroelastic matrix material at the fracture surfaces and connect the displacements via the fracture jump operator We denote the quantity τ 0 (x) = τ (x, t = 0) the initial fracture opening and treat it as a geometry parameter that is constant in time. However, it is important to remark that we do not resolve τ 0 geometrically. In the following, we address the material models used for the sub-scale description.

Poroelastic rock matrix
The fluid-saturated porous rock matrix is treated as a biphasic mixture consisting of the solid skeleton and the saturating fluid phase. We follow the notation proposed in [7] and denote quantities referring to solid and fluid phase respectively with We adopt the concept of effective stress and use Biot's quasi-static theory of linear consolidation [3] in a displacement-pressure formulation. This leads us to the strong form of the coupled equation system with the boundary conditions The subscripts ∂ D and ∂ N refer to Dirichlet and Neumann boundary conditions, respectively, pertaining to the solid displacement u and the pore pressure p M . Quantities with index * are given boundary values on the boundary ∂Ω M . In (3), σ is the total stress in the mixture, and it is defined via the constitutive relation with the shear and the bulk moduli G and K of the try solid rock skeleton, the Biot-Willis parameter α and the pore pressure p M . We assume small deformations and define the linear solid strain tensor ε = (u ⊗ ∇) sym . We execute the volumetric-deviatoric decomposition and define the volumetric strain tensor ε vol = e I/3 with e = tr ε and the deviatoric strain tensor ε dev = ε − ε vol . In (4), we identify q M as the seepage velocity, whereas Φ M is the amount of fluid stored in the saturated pore space. The corresponding constitutive relations are chosen as Defining the initial value of Φ M completes the initial boundary value problem. The material parameters used for the poroelastic model are summarized in Table 1.

Pressure diffusion along fractures
The fluid pressure diffusion along the fractures ∂ F k , k = 1, 2, . . . , n, is described in terms of the continuity equation The operator ∇ F evaluates the gradient in the plane that is tangential to ∂ F k . The strong form (9) is completed with the initial value of the fluid storage Φ F and Dirichlet and Neumann boundary conditions pertaining to the fluid pressure Here, we introduce the fluid velocity q F along ∂ F k and the fluid storage Φ F associated with ∂ F k with the volumetric fracture strain e F . K f is the fluid bulk modulus. Moreover, we introduce the leak-off q L or, in other words, the outflux of fluid from ∂ F k into the surrounding porous matrix, and we define the mass supplyq k from one fracture to another, in the case they intersect, aŝ whereq kl = −q lk is the influx from fracture l to fracture k along their intersection ∂ F k ∩ ∂ F l . Here, we use the Diracfunction 2 δ kl at the intersection of fractures k and l. The volumetric fracture strain is defined as whereby the volume change of the fracture due to elongation in tangential direction is ignored. For the fluid velocity q F along the fracture we suppose that the velocity profile has a quadratic shape across the thickness of the fracture. The Poiseuille flow assumption allows for computation of the effective seepage velocity q F in ∂ F k in form of Darcy's law as with the apparent fracture permeability τ 2 0 12 . 2 We define the Dirac-function as the distribution satisfying

Hydro-mechanical coupling between matrix and fracture space
Finally, we specify the interaction between the poroelastic matrix and the fluid-filled fractures. We, therefore, compute the volumetric fracture strain rateė F with (1) and (13) and define the leak-off q L as the outflux of fluid from the fracture space into the poroelastic matrix. Hence, we writė Vice versa, the fluid pressure in the fracture acts on the internal surface ∂Ω i,k M . The assumption of a continuous fluid pressure leads to the interface conditions k = 1, 2, . . . , n. Henceforth, we simplify notation and write p instead of p M and p F .

Scale transition
In the subsequent section, we aim to embed pressure diffusion in fluid-saturated fractured rock in the multiscale framework of computational homogenization and derive a macroscopic substitute model for the hydro-mechanically coupled problem on the meso-level. To this end, we define a SVE for the mesoscopic pressure diffusion problem and impose separation of scales in space and time. This means that macroscopic gradients have a wavelength that is much larger than the SVE size (L l). Similarly, due to the scaling of spatial gradients, processes within the SVE will occur at much higher rates than those occurring on the macroscopic level. We conclude that, under these conditions, transient pressure diffusion, brought on by macroscopic waves, occurs on the meso-scale only. In standard fashion, we assume the SVE to be periodic and apply periodic boundary conditions. Thereby, we a priori ensure that the amount of fluid stored within the SVE is conserved. It remains to discuss pressure diffusion through the SVE: Firstly, due to separation of time scales and negligible macroscopic time rates, pressure diffusion through the SVE will be stationary. Secondly, due to the linearity of the problem, it does not interfere with the transient processes and is, therefore, ignored. See [13] for a more thorough discussion.
Altogether, the macroscopic substitute medium can be described as a monophasic material. In fact, the attenuation properties of the macro-model represent viscoelastic damping of the Maxwell-Zener type.
We tackle the scale transition between meso-and macroscale with the aid of the volume averaging operator Here¯ represents the macroscopic counterpart of an appropriate subscale quantity . Taking into account the dimensional reduction Ω F → ∂ F, we obtain the hybrid-dimensional form of the volume averaging operator as Periodic boundary conditions on the SVE surface are imposed as  [16]. The surface tractions and the fluid outflux across ∂Ω M and ∂∂ F ,k are defined as t = σ ·n, q M = q M ·n and q F = q F ·n with the outer surface normal vector n. With this choice, it is obvious that ∂Ω M q M da + n k=1 ∂∂ F ,k q F τ 0 ds = 0. Thus, the fluid mass in the SVE is conserved. Moreover, we observe that the macroscopic strainε is the only loading applied on the SVE. This corresponds to the modeling assumption of a monophasic solid overall material model.

Weak form and macro-homogeneity condition
In the next step, we establish the weak representation of the balance equations (3), (4) and (9) using the notation for the generalized variational format presented in [16]. Hence, we seek solutions in the trial spaces U M and P = P M × P F of admissible displacements and fluid pressure fields 3 that are sufficiently regular in Ω M and ∂ F. Furthermore, we introduce the corresponding trial spaces of self-equilibrated fluxes T M , Q M and Q F that are sufficiently regular on ∂Ω + M and ∂∂ F + . We write the equations for finding u, which hold for any admissible test functions δu, δ p, δt, δq M , and and Using (12) and the constraintq kl = −q lk , we see that n k=1 e p F,k (q k , δ p) = 0. Moreover, using (15) 2 , we conclude that d p M ( p, δ p)+d p F ( p, δ p) = 0. In order to abbreviate notation we introduce q = (q M , q F ) and the corresponding set Q = Q M × Q F . For a unified description of the fully coupled pressure diffusion process on SVE level we add (24) and (25), and we write the SVE continuity equation solving for p, q ∈ P × Q as which holds for any admissible test functions δ p, δq ∈ P × Q . In (45) we introduce the SVE forms For the sake of completeness, we derive Hill's principle of macro-homogeneity for the given problem. We, therefore, choose δu →u and δ p → p. We add (23) and (45) and write in expanded form Hereby, we employ the periodic boundary conditions in the form of (20), (21) and (22). On the right-hand side of the macro-homogeneity condition (51) we can see that, under the chosen periodic boundary conditions, the macroscopic substitute medium is represented by a monophasic Cauchy continuum model. We use the divergence theorem and rewrite the macroscopic stress tensorσ as

Numerical model reduction for the RVE problem
In this section we aim to enhance the computational homogenization concept in a way that enables us to explicitly derive the macroscopic viscoelastic material properties of the pressure diffusion problem on the sub-scale. We, therefore, extend the NMR technique introduced in [12] towards the deformation-driven hybrid-dimensional pressure diffusion problem.

Decomposition of the sub-scale field quantities
The key ingredient to our NMR approach is the series expansion of the pressure field p(x, t) in terms of The pressure modes p a , a = 1, 2, . . . , N , form a linearly independent reduced basis of the space P . We show in the sequel that even rather small numbers N result in very good approximations of the pressure field and the depending quantities. The pressure modes satisfy the condition N a=1ξ a p a = 0 only by the trivial choiceξ a = 0, a = 1, 2, . . . , N . The quantitiesξ a are called mode activity coefficients. They represent internal variables of the macroscopic substitute medium. Altogether, the current state of the hybriddimensional medium on the meso-level at time t is defined by the overall strainε(t) and by the internal variablesξ a (t), a = 1, 2, . . . , N .
The pressure modes p a are generated via a Proper Orthogonal Decomposition (POD) of a set of S snapshotsp s , s = 1, 2, . . . , S, of the pressure distribution obtained from training computations. As described in [11,25], training computations take place on the SVE level and are driven by a macroscopic loadinḡ Hereby, the tensors B i , i = 1, 2, . . . , 6, represent the irreducible orthonormal basis of the symmetric strain tensor ε. The training function γ i (t) may, for example, prescribe stress-relaxation tests or a frequency sweep. The pressure snapshots are used to compute the system correlation matrix We solve the eigenvalue problem (g st − λ δ st ) v t = 0 and arrange the resulting eigenvalues in the order of decreasing absolute values |λ s |. We reduce the set of eigenvalues to the N S persisting eigenvalues λ a , a = 1, 2, . . . , N , for which holds |λ a | > 1e-6 |λ 1 |. Finally, we compute the N pressure modes as Due to the orthonormality of the eigenvectors v a s , the pressure modes are orthogonal such that In the hybrid-dimensional case we have to take into account that, usually, |Ω F | |Ω M |. In other words, the matrix contribution to the volume integral in (55) exceeds the fracture contribution by far. By consequence, the diffusion processes in the fracture space might be underrepresented in the POD procedure. One way to avoid this problem is to introduce scalar and dimensionless weights w M , w F > 0 to compute the correlation matrix as with w M w F . A self-evident choice would, for example, be based on the volume fractions such that w M = |Ω F | |Ω | and w F = |Ω M | |Ω | . In this manuscript, however, we choose another yet comparably simple concept and execute the POD for the different topologies Ω M and ∂ F separately. We, therefore, use S M snapshots covering the pressure diffusion processes in the matrix material and S F snapshots covering the pressure diffusion in the fracture space. We define the correlation matrices accordingly as Thus, the N = N M + N F members of the reduced set of basis modes consist of N M modes that are derived from (59) and N F modes that are derived from (60). It is important to remark that this procedure might, in general, lead to a non-orthogonal reduced basis and therewith to a basis that is larger than necessary. However, the detailed investigation of optimal identification procedure of the reduced basis is outside the scope of the present investigation. After having computed the reduced modal basis according to (56), we now decompose further mesoscopic fields according to (53) and write

ε(x,ε(t),ξ(t))
which hold for x ∈ Ω M . We have introduced the 4th rank strain and the 3 rd rank displacement localization tensors whereas C(x) is the elastic stiffness tensor for the dry solid skeleton. Moreover, we use the mode-dependent strain, displacement and stress tensors ε a , u a and σ a . All resulting fields on the meso-level depend linearly on the prescribed overall strainε and the mode activity coefficientsξ a , a = 1, 2, . . . , N . The contributions associated with the strain localization tensor E 0 or, in other words, depending on the macroscopic strain tensorε represent the instantaneous response of the dry linear-elastic solid skeleton under kinematic loading and zero mode activity, that isξ a = 0, a = 1, 2, . . . , N . We, therefore, solve for u i and t i , i = 1, 2, . . . , 6, from for all admissible test functions δu, δt ∈ U M × T M . The unknown fields ε a , u a and σ a are driven by the activity coefficientsξ a . They are computed by solving b = 1, 2, . . . , N linear-elastic eigenstress problems corresponding to the unit loading cases whilstε = 0. Since, p a is known, we solve for u a and t a from the system Finally, we compute the total macroscopic stressσ as Equation (70) allows us to identify the macroscopic effective stressσ eff and the macroscopic pressure-dependent stress σ p . It is particularly remarkable thatσ p is not necessarily a spherical tensor. This is in sharp contrast to the pressuredependent stress σ p in the meso-scale formulation (7). This is a consequence of the heterogeneity of the SVE: Any fluid pressure distribution on the sub-scale causes eigenstresses in the SVE which are, in general, anisotropic. Finally, it remains to identify a relation for the temporal evolution of the activity parametersξ a which we interpret in the following as internal variables of the macroscopic substitute medium.

Numerical experiments
We now aim to demonstrate the capability of our approach to predict the viscoelastic properties of the macroscopic substitute model in a reliable way. We have, therefore, implemented the hybrid-dimensional SVE problem in the Finite Element Software COMSOL Multiphysics and investigate three different cases. In example 1, we study pressure diffusion in the case of a so-called patchy saturation. In this case, pressure diffusion occurs only in the poroelastic matrix material in the absence of fractures. In example 2, we choose a SVE consisting of a simple fracture network embedded in an impermeable matrix. In this case, on the other hand, pressure diffusion only occurs in the fractures. Finally, we investigate the full hybrid-dimensional diffusion problem including pressure diffusion in both the matrix and the fracture network.

Example 1: patchy saturation
We start our considerations with the 2D SVE shown in Fig. 3, where we use the plain-strain assumption. The SVE consists of a homogeneous poroelastic material which is saturated with two different fluids: Water and gas. The material parameters are given in Table 2.
Because the SVE contains no fractures, the volume averaging operator is simplified as¯ = := M . Consequently, the weak form of the diffusion equation (25) vanishes, and it follows in (82) -(84) thatĀ F,ab = 0, B F,ai = 0,M F,ab = 0. The macroscopic stress response computes asσ = σ = σ M . We execute training com- putations in terms of stress relaxation tests on the SVE level.
In view of (54), we choose a maximum strainε max = −0.01, a ramp time t ramp = 1e-5 s and prescribe the three loading cases in two dimensions as In Fig. 4 we show examples of snapshots of the pore pressure distribution during the stress relaxation computation for the caseε 11 (t) = γ 1 (t) whilst all other macroscopic strain components equal zero. We observe an instantaneous pressure gradient between the water-saturated patch (high pressure) and the gas-saturated matrix (low pressure). This is due to the given bulk moduli K f of the involved fluid constituents. As time passes, the heterogeneous pressure causes local redistribution of pore fluid, and the pressure gradient is equilibrated via pressure diffusion. We execute the POD as described in Sect. 4 and identify a reduced basis with 9 members. After having computed the system matrices constituting the viscoelastic evolution Eq. (86) we can identify the characteristic frequencies of the involved Maxwell chains as C aa = [0, 1.15e + 0, 5.91e + 0, 1.61e + 1, 4.58e + 1, 1.30e + 2, 4.0e + 2, 1.75e + 3, 4.66e + 3] T 1/s.
Hereby, we use the plain-strain vector representation of the macroscopic strain defined asε = [ε 11 ,ε 22 , 2ε 12 ] T . Obviously, the macroscopic shear deformation does not cause any pressure diffusion. This can be explained by the fact that the solid phase material properties are those of a homogeneous porous rock with isotropic properties. Hence, the macroscopic shear strain causes a purely deviatoric deformation state in the entire SVE which does not interact with the fluid pressure according to (8).
In the next step, we use the reduced viscoelastic substitute model to recompute the macroscopic stress relaxation experimentε 11 (t) = γ 1 (t). In Fig. 5 we show the temporal evolution of (a) the internal variablesξ =Rχ and (b) χ. We observe that, as expected, the evolution of the viscoelastic variablesξ is strongly coupled whilst the variables χ show the typical relaxation behavior of uncoupled viscoelastic variables in a generalized Maxwell-Zener model. In order to validate the reduced viscoelastic substitute model, we consider the stress response computed for the SVE problem with full resolution to represent a reference solution. In Fig. 6 we plot the pressure dependent stress contribution σ p i j as defined in (70) for the reference solution and for the reduced viscoelastic model. We find an excellent agreement with the reference solution in the case of the stress relaxation experiment. Finally, we modify function γ 1 (t) and use a wavelet-type loading with two different central frequencies. The total macroscopic stress response of the SVE is successfully validated in Fig. 7 against the reference computation.

Example 2: fractures in impermeable matrix
In example 2, we investigate a simple 3D network of two intersecting perpendicular fractures saturated with water and embedded in a linear-elastic and impermeable matrix material, see Fig. 8. The used material parameters are given in Table 3. Similar to the previous example, this second example represents a special case of the modeling concept. Here, pressure diffusion occurs in the fractures only. Hence, the weak form of the diffusion equation (24) in Ω M can be neglected, several coupling terms in the remaining weak forms vanish.
In particular, it holds b u M ( p, δu) = 0 and d p F ( p, δ p) = 0. Moreover, we need to condense the system matrices that con- In analogy to the previous example and (87) we compute 6 stress relaxation experiments as training cases. Hereby, we choose t ramp = 1e-1 s andε max = −1e-6. We show snapshots of the pressure field at selected time instances during the relaxation experiment underε 11 (t) = γ 1 (t) in Fig. 9. We observe an instantaneous pressure increase in fracture 2 which is oriented perpendicular to the macroscopic perpendicular loading whilst fracture 1 remains nearly unaffected. As time passes by, fracture 2 and 1 exchange fluid mass which allows for a diffusive equilibration of the pressure gradients.

Example 3: fractures in permeable matrix
Finally, we investigate the full hybrid-dimensional diffusion problem, whereby we utilize the same SVE geometry as in example 2, see Fig. 8. However, the matrix is defined as a poroelastic material. Both, poroelastic matrix and fractures, are saturated with water. The material parameters are specified in Table 4. As in the previous examples we compute six stress relaxation experiments where we choose t ramp = 1e-4 s andε max = −1e-6. In Fig. 12 we show snapshots in time of the evolving pressure field during the uniaxial stress relaxation experimentε 11 (t) = γ 1 (t). We see that fracture 2 is compressed instantaneously whilst fracture 1 as well as the matrix material remain nearly unaffected. In other words, we observe a strong pressure gradient between fracture 2 and fracture 1 as well as between fracture 2 and the surrounding porous matrix. With increasing time the pressure gradient in the fracture space is equilibrated first. Afterwards, pressure diffusion occurs between fractures and matrix.
We use snapshots from all training computations for a POD and identify a reduced basis with 18 members. The characteristic frequencies are spread over 3 decades and read 1.14e + 1, 1.16e + 1, 3.05e + 1, 3.71e + 1] T 1/s (94) Again, mode a = 1 represents the equilibrium state of a homogeneous pressure distribution. It is important to remark that the basis in example 3 (18 modes) is significantly larger than the basis in example 2 (6 modes). The reason for this increase is that the pressure diffusion process in example 3, compared to example 2, is far more complex as it involves pressure diffusion in the fractures, in the embedding poroe-lastic matrix as well as the exchange of fluid mass between fractures and matrix. This complexity increase is mirrored in the number of basis modes.
In Fig. 13 we show example pressure modes. Whilst the modes a = 4, 6 are related to processes in the fracture space, modes a = 13, 15 describe processes in the matrix.
Finally, we validate the reduced viscoelastic basis against a reference computation of a uniaxial stress relaxation test on the basis of a fully resolved SVE. The resulting macroscopic total stresses are shown in Figs. 14 and 15, where we observe a very high accuracy of the reduced viscoelastic model.

Discussion
We establish a NMR procedure that enables us to identify macroscopic viscoelastic substitute models for strain-driven pressure diffusion problems in fractured porous rock in a reliable and numerically highly efficient way. The overall viscoelastic response is caused by squirt-type fluid flow which causes fluid pressure diffusion in interconnected fractures and in the embedding porous rock within a mesoscopic SVE. The fluid mass stored in the SVE is conserved throughout the entire process. This allows us to treat fluid pressure diffusion as a local process which is, on the macro-scale, only observable indirectly via overall frequency-dependent properties such as dispersion of the stiffness moduli and attenuation. Hence, the macroscopic substitute medium represents a monophasic viscoelastic solid. For the mesoscopic model we employ a hybrid-dimensional interface model treating fractures as fluid-filled open conduits embedded in a poroelastic matrix material. Pressure diffusion occurs along the fractures as well as in the poroelastic matrix. Moreover, exchange of fluid mass between fractures and matrix is incorporated. We employ volume averaging techniques to derive a computational homogenization scheme that links the het-erogeneous meso-scale to the viscoelastic macro-model. We propose a NMR technique that is based on the approxi-mation of the space-time dependent fluid pressure in the SVE by a linear combination of pressure modes. Hereby, the pressure modes span a finite-dimensional reduced basis. They are identified via a POD from a set of training computations. Evaluation of the mesoscopic continuity equation together with the pressure decomposition and accordingly decomposed further quantities allows us to derive the evolution equation of the desired viscoelastic substitute model which turns out to be of the Maxwell-Zener type. Finally, we investigate the performance of the proposed method in three different numerical examples. Throughout the examples, we observe a convincing agreement between NMR and reference solutions.
The main achievement of the present contribution is that we establish a general, yet reliable, method that enables us to compute seismic wave propagation on the macro-level with full access to the mesoscopic pressure diffusion effects but without the numerical drawbacks of the FE 2 concept. The numerically expensive solution of mesoscopic boundary value problems is restricted to 6 transient training computations, to 6 computations of the linear-elastic response of the dry solid skeleton and to the solution of N linear-elastic eigenstress problems related to the identified pressure modes. All mesoscopic computations are "off-line" precomputations and, therefore, do not burden the "online" computation of the macroscopic wave propagation.
In our ongoing research we will use the proposed method to execute forward simulations of propagation of seismic waves in a cross-hole tomography setting. The goal is to gain a better understanding for the correlation between seismic attenuation and the interconnectivity of stochastic fracture networks.