A poro-viscoelastic substitute model of ﬁne-scale poroelasticity obtained from homogenization and numerical model reduction

Numerical model reduction is exploited for computational homogenization of the model problem of a poroelastic medium under transient conditions. It is assumed that the displacement and pore pressure ﬁelds possess macro-scale and sub-scale (ﬂuctuation) parts. A linearly independent reduced basis is constructed for the sub-scale pressure ﬁeld using POD. The corresponding reduced basis for the displacement ﬁeld is constructed in the spirit of the NTFA strategy. Evolution equations that deﬁne an apparent poro-viscoelastic macro-scale model are obtained from the continuity equation pertinent to the RVE. The present model represents an extension of models available in literature in the sense that the pressure gradient is allowed to have a non-zero macro-scale component in the nested FE 2 setting. The numerical results show excellent agreement between the results from numerical model reduction and direct numerical simulation. It was also shown that even 3D RVEs give tractable solution times for full-ﬂedged FE 2 computations.


Introduction
A large class of natural and engineering materials are porous with a pore-fluid interacting with the solid skeleton. These materials are often micro-heterogeneous in the sense that there is a pronounced difference in the mechanical properties between the different micro-constituents, e.g. solid grains and pore fluid. Moreover, heterogeneity may occur on a much larger length scale. A typical example is porous rock which might be saturated by different pore fluids (e.g. patchy saturation models [25,28]) or which might be composed of zones with varying material properties of the solid skeleton (e.g. double porosity models [12,22]). In both cases, the mechanical properties, such as compressibility of the constituents and permeability of the pore space, are highly heterogeneous. The The straightforward approach is to model the problem on one single scale using Direct Numerical Simulation (DNS). However, resolving the substructural features is computationally very expensive and, therefore, some sort of homogenization strategy is preferable. When the intrinsic material properties are nonlinear and/or the fine-scale problem is inherently transient, it is necessary to resort to nested macro-and sub-scale computations (FE 2 ), whereby the sub-scale computations are carried out on a so-called Representative Volume Element (RVE) 1 in each quadrature point in the macro-scale domain, possibly within a given time step. Clearly, the purpose is to obtain macro-scale properties of engineering interest; hence, whether it is possible to avoid resolution via DNS of the fine-scale problem and accept the homogenized solution can only be assessed via some sort of goal-oriented error quantification.
However, it is widely realized that straight-forward application of the FE 2 -strategy can be computationally intractable for a fine macro-scale mesh, particularly in 3D. Therefore, there is significant interest in reducing the cost of solving the individual RVE-problems by introducing some kind of reduced basis, here denoted Numerical Model Reduction (NMR). In particular, we note strategies that are based on the superposition of "modes" that are characteristic for the RVE-solution fields. In the case of sub-scale small strain (visco)plasticity, various attempts have been made to approximate the inelastic strain field by a reduced basis consisting of so-called "inelastic modes". One of the early proposals is the so-called "eigendeformation-based reduced order homogenization" technique by Fish and coworkers [7,23], which in its turn relies on the Transformation Field Analysis (TFA) that was originally proposed by Dvorak and Benveniste [4]. A similar approach, coined Nonuniform Transformation Field Analysis (NTFA), was proposed by Michel and Suquet [20,21,26]. Recent developments are e.g. by Fritzen and coworkers [8][9][10][11], who extended the idea further within the framework of Proper Orthogonal Decomposition (POD) and applied it to visco-elasticity and, more generally, to a subclass of Standard Dissipative Materials.
Quite importantly, however, is the obvious fact that the richness of the reduced basis will determine the accuracy of the RVE-solution, which calls for error control. Here, we consider the full FE-solution as the exact one. An example of error estimation due to model reduction, although not in a homogenization context and for a PGD-basis, is Ladeveze and Chamoin [18]. More recently, Ekre et al. [6] investigated the transport of error from the micro-to the macro-scale model in a reduced FE 2 approach.
In recent years, several attempts have been made to establish multi-scale procedures for poroelastic media, e.g. Su et al. [27]. A recent contribution is by Khoei and Hajiabadi [17]. However, the straight-forward FE 2 -strategy turns out, as mentioned above, to be intractable for realistic large-scale problems in 3D. Jänicke et al. [13,15] studied the computational homogenization of heterogeneous poroelastic media under locally undrained conditions, i.e. fluid redistribution is only allowed locally whilst global fluid transport is suppressed. This "selective" homogenization approach allows to reduce the macro-scale problem to that of a single-phase viscoelastic material, whereby the "'mode coefficients" play the role of classical internal variables of the macro-scale model. More recently, this selective homogenization towards a viscoelastic macro-scale material was extended to account for Darcy-type seepage in fluid-filled fracture networks [14].
In this paper, we introduce the FE 2 approach with Numerical Model Reduction based on the complete, i.e. non-selective, upscaling of the poroelastic fine-scale model towards a poro-viscoelastic macro-scale model. The scale transition is formulated in terms of a Variationally Consistent Homogenization scheme with weakly periodic boundary conditions on the RVE problems. Inspired by the NTFA approach, we use training computations and POD to identify a reduced basis of pore pressure modes that allows us to shift all fine-scale computations from the "online" FE 2 solution phase towards an "offline" training phase. The resulting macro-scale model contains both the viscoelastic properties associated with local fluid redistribution on the RVE level as well as the hydro-mechanical properties of a macroscopic poroelastic medium. Therefore, it is called poro-viscoelastic. We benchmark the predictions emerging from the poroviscoelastic NMR procedure against results from DNS in several numerical examples and demonstrate the numerical efficiency of the model in a large-scale 3D example.
Throughout this paper, meager type is used to denote scalars, whereas bold type is used to denote vectors and (higher order) tensors. Single contraction is denoted by ·. For example, if a, b are vectors and A is a second order tensor, we have a ·b = a i b i , (A·b) = A i j b j , where the Einstein summation convention is used.
As to homogenization in the spatial domain, volume and surface averaging of an intensive field are denoted where the domain occupied by the RVE is Ω with boundary Γ . The macro-scale representation of is denoted¯ , and frequently it holds that¯ := . The paper is organized as follows: Sect. 2 gives a short review of the poroelastic fine-scale model. Section 3 is devoted to the two-scale modeling based on computational homogenization. The NMR strategy is presented in Sect. 4. Section 5 presents the numerical results for 1D, 2D and 3D settings including benchmark experiments. Finally, Sect. 6 concludes the paper with a summary and an outlook to future work.

Preliminaries
We give a brief summary of poroelasticity under geometrically linear and quasi-static conditions. Initially introduced by Biot [1,2], poroelasticity as a homogenized continuum model of solid-fluid mixtures with hydro-mechanical coupling can be seen as a special case of the more general class of materials described by the Theory of Porous Media (TPM) as proposed e.g. by De Boer [3] and Ehlers [5]. The primary global fields are the displacement of the solid skeleton, denoted u, and the intrinsic pore fluid pressure, denoted p.

Constitutive relations for the sub-scale constituents
The constitutive relations for the equilibrium stress σ and the seepage velocity w become: where ε is the strain 2 and ζ is the pressure gradient. The material properties are as follows: E is the fourth order elasticity modulus tensor which is defined in the simplest case of isotropy by the constant parameters G and K , which are the shear modulus and compression modulus, respectively. K is the second order constant permeability tensor, defined (in the simplest case of isotropy) by the scalar permeability coefficient k. Finally, α is the Biot coefficient. We also introduce the storage function Φ that represents the current volume fraction of fluid where φ is the initial porosity of the solid phase, and β is the intrinsic compression compliance of the pore fluid. Note that φ is taken as constant and independent of the deformation. As to the spatial variation of these parameters, they may be strongly heterogeneous on the fine scale.

Strong and weak format of the fine-scale problem
Restricting to quasi-statics, we thus seek the displacement u(x, t) : Ω × R + → R 3 and the pore pressure p(x, t) : together with the constitutive relations in (2) and (3) and the initial condition Clearly, (4a) represents static equilibrium, whereas (4b) expresses mass conservation of the two-phase solid-fluid mixture.
The standard space-variational format corresponding to (4) reads: Find (u(•, t), p(•, t)) in the appropriately defined spaces U × P that solve where U 0 and P 0 are the appropriately defined test spaces.

First order homogenization in the spatial domain
We replace the single-scale problem in Sect. 2 by a two-scale problem upon introducing (i) running averages in the weak format and (ii) scale separation via first-order homogenization. The first step is to replace the space-variational problem in (6) by that of finding (u(•, t), p(•, t)) ∈ U FE 2 × P FE 2 that solves where the pertinent space-variational forms in (7) are given as Here, (v; q) represent arbitrary but admissible displacement and pore pressure fields. Obviously, these space-variational forms represent running averages on square/cubic domains Ω with the side length L and located at each macro-scale spatial pointx ∈ Ω. We note that the introduced averaging results in a change of the original problem. The weak format thus serves as the basis for variationally consistent homogenization as discussed below. Next, in order to define the two-scale solution and test spaces U FE 2 , V FE 2 and P FE 2 , Q FE 2 associated 3 with FE 2 , we introduce scale separation and first-order homogenization. The standard approach is then to decompose the sub-scale fields u and p into a macro-scale part, u M and p M , and a micro-scale or fluctuating part, u μ and p μ , within each RVE such that We are now in the position to define the two-scale trial and test spaces whereas The trial and test spaces for the macro-scale problem are denotedŪ,V and for the micro-scale problem U μ ,i , V μ ,i , and P μ ,i , Q μ ,i . We introduce the notation Hence, there exist components (u μ , p μ ) and corresponding spaces (U μ ,i , P μ ,i ) for each single RVE Ω ,i .
We now restate the two-scale problem (7) as follows: Find where we used the interpretation δu μ = δu μ i and δ p μ = δ p μ i inside Ω ,i .

Macro-scale problem
We obtain the macro-scale (homogenized) problem from (13) upon choosing purely macroscopic test functions; hence or, more explicitly, where we introduced the variationally consistent macro-scale (homogenized) fields In addition,t p andw p are defined as the suitably homogenized quantities on the Neumann boundary parts. Note that, according to (15b) and (4b), the total macroscopic seepage velocityw can be computed as Also note that the homogenized stress can be expressed in standard fashion as the surface integral

Micro-scale problem on a representative volume element
In the previous section, the macro-scale problem was defined, and it was established that the macro-scale fields in (16) depend on (u, p) inside the pertinent Ω . We shall now consider the local problem on a single RVE Ω ,i to determine the implicit relation (u To this end, we consider the problem in (13) for a single test function δu For brevity, we henceforth drop index i since we now restrict to one single RVE.
The decomposition in (9) presumes that the solution (u, p) can be uniquely decomposed into macroscopic parts (u M , p M ) and fluctuations (u μ , p μ ). This imposes constraints on U μ and P μ in (13). In order to be explicit, we shall consider an expanded format of the problem. To this end, we introduce the richer spaces as follows: where H 1 (•) is the Sobolev space of functions with square integrable derivatives.
In order to ensure a uniquely solvable RVE-problem, we next adopt the following additional model assumptions: -We assume that u μ and p μ are periodic, i.e.
where we introduced the "difference operator" is the corresponding "mirror point", see Fig. 1. Upon introducing, for the sake of brevity, the RVE-forms we may express the variational (weak) statements 4 of the micro-periodicity constraints in (20) as whereby it is obvious that the Lagrange multiplier δλ in (22a) represents a field of virtual boundary tractions on Γ + , whereas the Lagrange multiplier δμ in (22b) represents a field of virtual boundary fluxes, see Remark 2 below. Here, we introduce the test spaces where L 2 (•) is the space of square integrable functions. Upon introducing the expressions for u M and p M from (9) on the RHS in (22a) and (22b), respectively, we obtain the explicit expressions where it was used that Γ [x −x] dΓ = 0. Hence, we conclude that it is only the datah andζ that will have an effect on the periodicity conditions.
-We assume thatū andp physically play the role of surface-averaged quantities as follows: which conditions are expressed weakly as The space-variational RVE-problem 5 can now be stated as follows: For a given historyū(t), It is possible to identify the tractions t and the flux w on the entire RVE-boundary Γ as follows: (28b) In particular, it appears that the spaces T and Q represent "self-equilibrating" tractions and fluxes, respectively.

Remark 1
Setting δu ∈ R 3 in (27a), we obtain the solution λ = 0. Nevertheless, it is still necessary to satisfy the constraint (27e) in order to prevent Rigid Body Modes (RBM), in this case translation. However, the macro-scale stressσ is invariant to any constant translationū within the RVE, which means thatū can be prescribed to any given value, saȳ u = 0. Moreover,σ is invariant to the skew-symmetric part ofh since ε[u] =ε + ε[u] μ ; hence it is onlyε :=h sym that has to be fed to the RVE problem as data.

Remark 3
Setting δu =u, δ p = p and summing up (27a) and (27b), we are able to derive the relation Moreover, taking (16)-(18) into account together with (28) and choosingλ = 0,μ = −Φ, we can derive the identity This expression is precisely the macro-homogeneity condition; therefore, it is a "built-in" property of the problem formulation.

Decomposition of sub-scale fields
The standard decomposition of the sub-scale fields were given in (9). However, in order to utilize NMR for the purpose of solving the RVE problems in an efficient manner, it is convenient to introduce the alternative decompositions 6 where the fields (u stat , Remark We note that u stat = u M and p stat = p M in general; hence,ũ μ = u μ andp μ = p μ . Moreover,λ stat = 0 and μ stat = 0, and it is concluded thatλ =λ μ = 0, whereas μ =μ μ = 0. Although it is possible to excludeλ stat and μ stat a priori from the formulation of the sensitivity problem, these variables are retained henceforth in order to preserve symmetry of the problem formulation. Since the problem in (33) is linear for each given time, t serves only as a load parameter, we may compute the solu- while noting that u M and p M can be expanded in a similar fashion as The following set of problems are then solved "offline" as a preliminary to the model reduction. It is possible to sequentially solve for ( p stat , μ stat ,μ stat ) and (u stat , λ stat ,λ stat ), in terms of the sensitivities, from the equations (33b), (33d), (33f) and from (33a), (33c), (33e) respectively: which has the trivial solutionp (ε) -Solve for sensitivity fields (û (p) ,λ (p) ,λ (p) ) from -ε = 0,p = 0,ζ = 0: -Solve for sensitivity fields (p

Numerical model reduction based on the NTFA strategy
We first note that the fluctuation fieldsũ μ ( Once again,λ μ is kept as an unknown variable for the sake of symmetry although it is known that the solution must bẽ λ μ = 0. The NMR-strategy for the solution of (42) is then outlined as follows: in the spatial domain, where {p a (x)} M R a=1 is a set of linearly independent global basis functions called "pressure modes". The pressure modes are identified in a training phase by solving (42) with different loadingsu stat andṗ stat with a POD of extracted snapshots ofp μ . The "exact" solutionp μ ∈ P 0 is then approximated bỹ p μ R ∈ P 0 ,R as follows 7 : where {ξ a (t)} M R a=1 are "mode activity" parameters. For computational efficiency, M R should be a finite, i.e. reasonably small, number. It turns out that the parameters ξ a (t) represent internal variables which define the viscoelastic contribution to the macro-scale model. Most importantly, however, the response of the poroelastic RVE depends on the loading in terms of the macroscale quantities {ε,p,ζ } as well as on the parameter set {ξ a (t)} M R a=1 . -Establish NMR-approximations for the other fluctuation fields in the spirit of NTFA, i.e. introduce the expansions whereby it is noted that the same mode activity parameters are employed to represent all the fields. Now, the key property is that it is possible to determine the mode basis sets introduced in (44) in terms of the pressure basis {p a (x)} M R a=1 which is identified in the training phase via POD. In order to do so, let us first consider the system defined by (42a), (42c), (42e) and insert the NMR-expansions for the fields {p μ R ,ũ μ R ,λ μ R ,λ μ R }. Since these equations must hold for any t, the coefficients of each mode activity parameter ξ a (t) must vanish; hence, the problem reduces to that of finding the modesû a ∈ U ,λ a ∈ T ,λ a ∈ R 3 for a = 1, 2, . . . , M R from the system Now, since U ,R ⊆ U and T ,R ⊆ T , we may choose δu =û b in (45a), δλ =λ b in (45b), and δλ =λ b in (45c) and combine the results, to give the useful (as will be shown later) relation Noting that Remark Since it is only the NMR-modesû a andp a that occur in (48), there is no need to actually compute the basis modeŝ u a ,μ a as part of the "online" algorithm.
It is convenient to rewrite (48) in the abbreviated form or, in matrix format, as where the matrix entries are defined aŝ Moreover, the RHS of (49), representing the macroscopic loading, is given aŝ When (50) has been solved for ξ a (t), we may express the pertinent macro-scale fields as Finally, it is convenient to derive the spectral form of the evolution equation (50). Thus, we execute a basis shift {ξ } →f forŜ andM as presented in [14]. The matrixR contains the eigenvectors of the generalized eigenvalue problem. We may writė

Special cases
Viscoelastic macro-scale model: An important special case of the model above is defined by the constraintsμ = 0 and ζ = 0, which correspond to the situation of an "undrained" RVE, i.e the value of the storage functionΦ is stationary. Fluid redistribution is only possible "locally", i.e. within the RVE. The conditionμ = 0 infers that the appropriate value ofp = p is part of the solution in the sense that it is computed as post-processing when p is known. In addition, the constraint equation (27f) becomes obsolete. Altogether, these model choices give rise to the single-phase viscoelastic macro-scale model proposed in [13,14]. The procedure may, therefore, be called "selective" homogenization. The macroscale problem is then reduced, as compared to the coupled problem for (ū,p) in (15), to the equilibrium problem that is pertinent to a single-phase continuum whereσ , as given in (54a), simplifies tō The fieldsσ (ε) i j andσ a are still given in (55a) and (55d), respectively.
The corresponding NMR-reduced RVE-problem is still that in (49) with the RHS simplified aŝ The problem of determining the stationary solution is simplified accordingly.
Hence, the poro-viscoelastic character of the macro-scale model is preserved.

Computational cost of the NMR model
To conclude this section, we address the computational cost of the NMR model. As mentioned in the previous sections, the key advantage of the method is that all computational efforts on the micro-scale can be executed as pre-computations. The computational cost for the training can be estimated as where the first contribution to (64) refers to the n training transient training computations for the snapshot generation whilst the second contribution refers to the linear-elastic computations needed to determine the system matrices in (59). In the 3D case, on may choose n training = 10 (6 strain, 1 pressure, and 3 pressure gradient loading cases). The number n it refers to the number of iteration steps needed to solve the transient problems, n S is the number of degrees of freedom of the RVE problem, and M R is the number of modes that span the reduced basis. Once the NMR model is trained, it can be re-used to solve an arbitrary number of macroscopic initial boundary value problems on Ω. The "online" cost of a macroscopic computation can be estimated as where n M is the number of macroscopic degrees of freedom.
In each Gauss point, the M R evolution equations need to be updated. For simplicity reasons, we assume that the number of Gauss points corresponds to n M . The gain in computational efficiency of the NMR model can be be quantified if the computational cost of the DNS solution is taken into account as where n DNS is the number of degrees of freedom of the DNS model. Moreover, we estimate the computational cost of a conventional FE 2 model with as nested solution of microand macroscopic boundary value problems as For the sake of simplicity, we assume n it to be of the same order of magnitude in all three cases. With that, the speedup from DNS to conventional FE 2 can be estimated as where n M n DNS 1, n M n S n DNS ≤ 1, and n S n DNS 1. However, the latter enters (68) only with power 1. By contrast, the speedup from DNS to NMR can be estimated as In particular in the 3D case, "much smaller" means in fact "several orders of magnitude smaller".
Finally, one may quantify the gain in computational efficiency from fully-nested conventional FE 2 to NMR by evaluating Altogether, we conclude that the computational savings of using NMR compared with DNS and FE 2 scale with the complexity of the problem and, thus, can be huge. In Sect. 5, we investigate computational costs in terms of examples at increasing complexity more explicitly.

Numerical results
In this section, we validate the proposed NMR method and demonstrate its performance in a large-scale reduced FE 2 computation. The material parameters used to define the poroelastic fine-scale model are listed in Table 1. All simulations were carried out running the Finite Element package COMSOL Multiphysics via Matlab livelink on a standard laptop computer.
In the subsequent simulations, mechanical loads are applied by means of a linear ramp function and kept constant after t = t ramp . Thus, we define the loading function (71)

Consolidation of a layered poroelastic medium
The first numerical example aims to provide a proof of concept of the proposed poro-viscoelastic NMR procedure. For simplicity reasons, we choose a periodically layered poroelastic medium. The 1D RVE shown in Fig. 3a represents the unit cell of the layered medium with L = 1 m and d = 0.25 m. Strains in the lateral directions equal zero. The material parameters are chosen according to set 1 given in Table 1  Remark -We observed that it is difficult to find a reduced basis with a reasonably low approximation error if the POD includes snapshots of all loading cases (ε,p,ζ ). The reason is that the pressure snapshots for the the loading cases (ε,p) are strongly different from the snapshots for theζ loading case. As a simple work-around, we execute two separate PODs, the first based on the snapshots from the (ε,p) loading cases, whilst the second is based on the snapshots from theζ loading case. The output of both PODs is combined to one unified reduced basis. Obviously, this procedure does not guarantee orthogonality of the reduced modes. Nevertheless, the approximation error turns out to be negligibly small as demonstrated in the subsequent validation experiments. -The characteristic frequencies given in (72) correspond to relaxation times of a Maxwell-Zener model. As shown subsequently, the chosen reduced basis allows for a high-fidelity approximation of the RVE properties incorporating (a) relaxation processes over many decades in time and (b) three strongly different macroscopic loading scenarios (ε,p,ζ ). The complexity of the RVE properties explains the maybe surprisingly high number of 22 modes.
Examples of pressure modes are shown in Fig. 3b. Finally, all sensitivities and system matrices that define the poro-viscoelastic substitute medium are computed from the definitions presented in Sect. 4. We are now in the position to utilize the poro-viscoelastic substitute model for the simulation of the macroscopic consolidation experiment shown in Fig. 3c as a reduced FE 2 computation. For validation purposes, we also compute a DNS of the macroscopic consolidation experiment which serves as reference. The reference computation is performed with full geometrical resolution, i.e. all layers contained in the macroscopic sample are discretized. In accordance with  Firstly, we compare the pore pressure distributionp(x) in the macroscopic computation domain at different times during the consolidation experiment, see Fig. 4a. At t = 1e+0 s, the reference computation shows pressure gradients between the layers of the fully resolved sample that are equilibrated in the course of time. Hence, at t = 1e+2 s, the global pressure decay due to the drained boundary conditionp = 0 atx = 0 becomes dominant. By construction, the reduced FE 2 computation, using the poro-viscoelastic NMR model, is unable to resolve the fluctuations (wavelength L ) of the pore pressure field at t = 1e+0 s. Instead, the NMR model correctly predicts the pore pressure that is measured at the RVE boundaries. This is in line with (25) 2 where we constrain the macroscopic pressure to equal the surface average of the microscopic pressure. Also, at t = 1e+2 s and t = 1e+3 s, the poro-viscoelastic computation matches the reference computation with very high accuracy.
Secondly, we investigate the displacementū measured at the surfacex = 0, see Fig. 4b. We observe three distinct phases of the consolidation experiment: (i) The loading phase for t < 1e-4 s, (ii) the phase where the inter-layer pressure gradients are equilibrated by local redistribution of the pore fluid for t ∈ [1e-2, 1e+1] s ("viscoelastic creeping"), and (iii) the phase of macroscopic drain-off due to the boundary conditionp = 0 atx = 0 at t > 1e+1 s ("poroelastic creeping"). In all phases, the poro-viscoelastic NMR model predicts the material behavior with very high accuracy.
The macroscopic drain-offw(x = 0) during the consolidation test as well as the change in fluid contentΦ − φ stored in the pore space at several evaluation points are plotted in Fig. 5. In both cases, the poro-viscoelastic NMR model is in excellent agreement with the reference computation.

Role of the RVE size
An important issue in computational homogenization is to decide whether the chosen volume element is large enough to be considered as sufficiently representative. In standard firstorder homogenization of single-phase materials, the size of an RVE is determined mainly in terms of stochastic representativity and boundary layer effects, e.g. due to over-stiff boundary conditions, whereas computational homogenization of poroelastic media involves an additional imprinted internal length scale, namely the diffusion length.
In the subsequent study, we aim to investigate the effect of RVE size on the macroscopic poro-viscoelastic properties of the homogenized model. We consider a 2D poroelastic medium (plain strain) consisting of stiff, circular inclusions which are embedded in a softer matrix. The material parameters are chosen according to set 1 in Table 1. The RVEs are constructed such that they contain one single inclusion with radius r = 0.175 m per 2D volume l 2 . Hence, the volume fraction of the inclusions is for n = 0.385. The different RVE sizes, |Ω | = {1, 4, 25} m 2 , are shown in Fig. 6a. For each RVE size, we generate an ensemble of five periodic RVEs with randomly distributed inclusions. The minimal distance between the inclusions is constrained to be larger than 0.1 r . For each individual RVE, we perform training computations with a POD of the extracted snapshots and compute the properties of the corresponding poro-viscoelastic NMR model. We show the pressure modes associated with the 16 lowest characteristic frequencies in Fig. 7  Note that the NMR model used for the subsequent simulations contains additional pressure modes with higher characteristic frequencies. Typically, we employ 27 pressure modes for each RVE. For the sake of brevity, the higher pressure modes are omitted in Fig. 7.
First of all, we remark that choosing a periodic RVE containing one single inclusion results in a strong scattering of the transient behavior, see Figs. 8, 9, and 10a. This might be somewhat surprising since the only difference between these RVEs with |Ω | = l 2 is the position of the inclusion. In the case of linear elasticity, it is clear that the homogenized properties of a periodic RVE under periodic boundary conditions that contain one single inclusion do not depend on the position of the inclusion. Obviously, this is not true in the present case. The explanation for this behavior is the presence of the pressure diffusion mechanism. Since only the pore pressure fluctuation in(20) 2 is periodic but not the prescribed macroscopic partp M , the position of the inclusion controls the diffusion behavior close to the RVE boundary which results in the scattering behavior as seen in Figs. 8, 9, and 10a.
The scatter is significantly reduced if we increase the RVE size. This holds, in particular, for the macroscopic stressσ 11 under the loadingε 11 (t), see Fig. 8, and the porosity increasē Φ − φ under the loadingp(t), see Fig. 10. To a lesser extent, this also holds for the effective seepage velocityw 1 underζ 1 (t), see Fig. 10. Moreover, for the larger RVE size we observe that the stationary state (t → ∞) is reached much later than for the smaller RVE size. This is related to the larger diffusion length associated with a larger RVE size. We show in the subsequent experiments that, nevertheless, all RVE sizes are suitable for predicting the transient behavior. Most importantly, however, we want to highlight that all NMR predictions (lines in Figs. 8,9,10) are in excellent agreement with the reference computations for the chosen RVEs (dots in Figs. 8,9,10).
Finally, we want to employ the identified NMR rules to execute reduced FE 2 simulations of the macroscopic consolidation experiment shown in Fig. 6b for the two different sets of material parameters given in Table 1 use as reference solution. The reference sample for DNS is generated according to Fig. 6b. It is periodic inx 2 -direction and contains 50 stochastically distributed inclusions. The reference computation is carried out with full geometrical resolution of the inclusions. The boundary conditions on the macroscopic problem aret = [t p , 0] T atx 1 = 0 with the prescribed tractiont p = 1e+6 γ (t) and t ramp = 1e-4 s. Moreover,p = 0 atx 1 = 0 (drained boundary). All other boundaries undergo symmetry conditions, i.e.ū · n = 0 and w = 0 (undrained boundary). The results obtained for parameter set 1 are shown in Fig. 11. As for the 1D case, we observe a pronounced fluctuation of the pore pressure fieldp for t = 1e-1 s that decays rapidly over time. The associated process is local redistribution of pore fluid between the inclusions and the matrix material in their vicinity. With the passage of time, the drained boundary conditionp(x 1 = 0) = 0 induces a macroscopic transport of pore fluid over the full length of the sample. Note that the reduced FE 2 solution shown in Fig. 11a is computed using the largest RVE size |Ω | = 25 l 2 . The agreement with the reference computation is very good. Nevertheless, it is important to remark that a certain deviation from the reference computations is to be expected due to the fact that the reference geometry itself is stochastic.
In Fig. 11b we plot the displacementū 1 (x = 0). Here, we compare ensemble averages of the RVE sizes |Ω | = {1, 4, 25} l 2 . Whilst the displacement related to local redistribution at t = 1 s is not very pronounced in this example, the NMR predictions of all RVE sizes show a very good agreement with the reference computation of the transient behavior. The only systematic deviation from the reference computation is the slight underestimation of the settlement in the fully equilibrated state by the smallest RVE size.
The same consolidation test has been performed with parameter set 2. The pressure distribution and the settlement are plotted in Fig. 12. Because the material properties G, K , and β are chosen to represent a significantly softer material than in parameter set 1, the amplitude of the pressure fluctuation at t = 1e-1 s is, compared to the overall pressure level, much larger. This is also reflected in the displacement u 1 (x = 0) where the local redistribution process at t = 1 s ("viscoelastic creeping") is much more pronounced than for parameter set 1. In this case, all NMR computations have been carried out on the basis of the ensemble consisting of five RVE realizations of the largest size. Taking into account that also the reference geometry itself is stochastically disturbed, we find again a very good agreement of the reduced FE 2 computation with the reference computation.

Consolidation experiment with inhomogeneous boundary conditions
We conclude our investigations with a 3D consolidation experiment with inhomogeneous boundary conditions. We assume a poroelastic RVE that contains 8 spherical inclusions, see Fig. 13a, and use the material parameters of set 2, see Table 1. The macroscopic problem is displayed in Fig. 13b. The loading at the top surface is applied via the prescribed tractiont = [0, 0,t p ] T on a rigid plate. We choosē t p = −γ (t) 1e+6 Pa and t ramp = 1e-4 s. The remaining part of the top surface is traction-fee and drained (p = 0). Symmetry conditions apply to all other boundaries. Hence, u · n = 0 andw = 0. We follow the NMR procedure and perform training computations to identify all relevant sensitivities of the poroviscoelastic substitute model. Examples of pressure modes are displayed in Fig. 14. Finally, we are able to execute the reduced FE 2 computations. In Fig. 15a we display the effective seepagew in the sample at different times. The arrows represent the seepage velocity vector. In Fig. 15b, we show, as a post-processing result, the local seepage velocity w as observed in an RVE adjacent to the evaluation point P for the same instants of time. We see that, at t = 60 s, the seepage is active essentially in the horizontal direction in the part of the sample which is exposed to the strong inhomogeneity of the loading. By contrast, the seepage at t = 3600 s is active mainly in the vertical direction towards the drained boundary. This behavior is also reflected by the streamlines in the RVE, see Fig. 15b.
It is important to remark that, for the given example, the fully resolved macroscopic geometry would con- tain 2,000,000 inclusions. A DNS of the sample would be far beyond typically available computing and memory resources. By contrast, the reduced FE 2 simulation based on the poro-viscoelastic NMR approach requires only a short computation time with a low memory utilization. The computational effectiveness is discussed in more detail in the subsequent section.

Computation time
To investigate computational effectiveness, we compare the explicit computation times between DNS and NMR solutions for the previous examples. All simulations were carried out on a standard laptop computer. The computation times are reported in Table 2 for the 1D case (cf. Sect. 5.1), the 2D case (cf. Sect. 5.2) and the 3D case (cf. Sect. 5.3). The computation times are evaluated for the solution of (i) one load case on one single RVE and of (ii) the macroscopic consolidation experiments. We observe that, with increasing complexity of the problem, the NMR model becomes more and more favorable. In practice, the computational effort to solve the RVE problem by NMR does not depend on the dimensionality of the problem. The speed-up factor for solving one single 3D RVE problem is 500. With NMR, also the 3D macroscopic consolidation experiment can be solved within a few minutes, although the relaxation behavior of the RVE problem gives rise to a considerable number of evolution equations. By contrast, the DNS of the 3D consolidation test, containing 2,000,000 inclusions, exceeds the capacities of the used computer resources by far. It is, therefore, impossible to quantify the speed-up in that case. These results point to the high computational effectiveness of the proposed method.

Concluding remarks: outlook
To provide increased modeling capability at the computational homogenization of fine-scale poroelasticity, we have considered the situation where the pore pressure field possesses macro-scale as well as sub-scale parts. The corresponding RVE problems will become more expensive as part of the FE 2 algorithm. Accordingly, it was expected that the NMR approach is much more efficient. Indeed, the presented numerical results show that the NMR method provides excellent agreement with the reference solution obtained from DNS for a quite small number of pressure modes. Moreover, it was shown that even 3D RVEs employed as part of the FE 2 algorithm give tractable solution times. As to future work, it is of significant interest to evaluate the importance of including the macro-scale viscoelastic variables, i.e. to discern when it is sufficient to assume that the RVE problem is reduced to what is coined the stationary problem. Such a reduction would mean that the macro-scale model is purely poroelastic. Another aspect that requires further investigation is the effect of the compressibility of the pore fluid on the NMR algorithm. Whilst in rock mechanics applications the bulk moduli of fluid and solid phase are usually of the same order of magnitude, the assumption of an, in practice, incompressible pore fluid is common for the description of "softer" porous media such as biological tissues or soft clay with full water saturation.