Anisotropic dual-continuum representations for multiscale poroelastic materials: Development and numerical modelling

Dual-continuum (DC) models can be tractable alternatives to explicit approaches for the numerical modelling of multiscale materials with multiphysics behaviours. This work concerns the conceptual and numerical modelling of poroelastically coupled dual-scale materials such as naturally fractured rock. Apart from a few exceptions, previous poroelastic DC models have assumed isotropy of the constituents and the dual-material. Additionally, it is common to assume that only one continuum has intrinsic stiffness properties. Finally, little has been done into validating whether the DC paradigm can capture the global poroelastic behaviours of explicit numerical representations at the DC modelling scale. We address the aforementioned knowledge gaps in two steps. First, we utilise a homogenisation approach based on Levin's theorem to develop a previously derived anisotropic poroelastic constitutive model. Our development incorporates anisotropic intrinsic stiffness properties of both continua. This addition is in analogy to anisotropic fractured rock masses with stiff fractures. Second, we perform numerical modelling to test the dual-continuum model against fine-scale explicit equivalents. In doing, we present our hybrid numerical framework, as well as the conditions required for interpretation of the numerical results. The tests themselves progress from materials with isotropic to anisotropic mechanical and flow properties. The fine-scale simulations show anisotropy can have noticeable effects on deformation and flow behaviour. However, our numerical experiments show the DC approach can capture the global poroelastic behaviours of both isotropic and anisotropic fine-scale representations.


Introduction
Numerical modelling of multiscale, poroelastically coupled materials can be challenging due to inherent length scale heterogeneities and multiphysics behaviours. Explicit modelling approaches allow one to account for each length scale directly within a model. This representation can therefore provide accurate and detailed descriptions. However, the number of degrees of freedom needed for direct models of multiscale, poroelastic materials can make simulation computationally prohibitive. Further, particularly within the subsurface, the data needed to populate such explicit approaches may be sparse.
Implicit models alleviate the problems associated with explicit models, at the expense of abstraction of local scale physics. One such modelling concept is the dual-continuum (DC) model, originally attributed to Barenblatt et al. (1960). This implicit approach has been used successfully within the context of flow modelling in a variety of subsurface engineering settings (Gerke and Van Genuchten 1993;Wu et al. 2002;Reimus et al. 2003;March et al. 2016). In the DC paradigm, one continuum represents a high storage, low permeability material (e.g. matrix), whilst the other represents a low storage, high permeability material (e.g. fractures).
This work concerns the dual-continuum modelling of multiscale, poroelastic geomaterials. We address two knowledge gaps associated with the DC modelling paradigm within the context of poroelasticity: First, we develop the poroelastic DC modelling approach. We introduce the underlying modelling assumptions, whilst considering the material symmetry and mechanical properties of the constituents in the process. With respect to the latter, previous poroelastic DC models have, for the most part, assumed isotropy of the continua and bulk material (Berryman and Wang 1995;Khalili and Valliappan 1996;Loret and Rizzi 1999;Choo and Borja 2015). However, rock formations are well known to exhibit anisotropic properties (Snow 1969;Price and Cosgrove 1990;Babuska and Cara 1991). Recent work by Zhang et al. (2019) showed that anistotropic permeabilities can have measurable impacts on the flow-patterns in poroelastic dual-continuum materials. Further to anisotropy, and in the case of fractured materials, the fractures themselves can have intrinsic mechanical properties owing to local asperities and/or bridging material between fracture faces (Olsson and Barton 2001;Lemarchand et al. 2009;Jaeger et al. 2009). Intrinsic mechanical properties of both continua have been considered for isotropic materials in the works of Elsworth and Bai (1992), Berryman (2002), Berryman and Pride (2002) and Nguyen and Abousleiman (2010). In this work we further explore the impact of the anisotropic elasticity, in addition to permeabilities, on dual-continuum responses.
Incorporating anisotropic and intrinsic properties can be done at the constitutive modelling stage. In the following, we add to a micromechanically derived anisotropic constitutive model by Dormieux et al. (2006). Contrary to the model by Dormieux et al. (2006), we incorporate linear (poro-) elastic properties for the low storage, high permeability continuum at the microscale. In this case both continua have intrinsic stiffness properties. Following homogenisation, the resulting model, complete with expressions for the effective parameters, is an anisotropic, dualstiffness constitutive model. Previous isotropic constitutive models reviewed in Ashworth and Doster (2019b) can then be recovered under isotropy and void-space assumptions on the general anisotropic, dual-stiffness constitutive model derived herein.
Second, with the derived poroelastic constitutive model, we proceed to numerical modelling. We investigate whether the DC representation is able to capture the global poroelastic behaviours of a fine-scale explicit model at the DC modelling scale. Whilst work has gone into testing and validating the DC concept for the flow problem (e.g. Lewandowska et al. 2004;Egya et al. 2019), little has been done to asses validity of the poroelastically coupled DC approach. Further, we discuss several considerations to ensure meaningful interpretations between the two modelling approaches.
To summarise, our aims are twofold. First, in Section 2, we use a homogensation approach and develop a previously introduced anisotropic dual-continuum constitutive model. For this development we allow both continua to have (anisotropic) intrinsic mechanical properties. Second, we perform numerical modelling of the poroelastic dual-continuum concept, investigating its validity against fine-scale explicit representations. To do so we introduce the hybrid numerical framework used to perform the numerical study in Section 3. We present the numerical tests, modelling considerations and test results in Section 4. For our study we consider numerical test cases as conceptualisations of naturally fractured rock samples that satisfy certain representative elementary volume (REV) requirements. Our results show that the DC model is capable of capturing the global poroelastic behaviours of isotropic and anisotropic fine-scale equivalents. Finally, we offer conclusions and recommendations for future work in Section 5.

Homogenisation of the dual problem
In the following we develop the anisotropic, dual-stiffness constitutive model for a poroelastic DC material. To do so we expand the homogenisation approach originally proposed by Dormieux et al. (2006) by including intrinsic mechanical properties for the low storage, high permeability continuum.
Whilst this work strictly assumes linear poroelasticity, the inclusion of stiffness properties are necessary for extensions to non-linear modelling of materials . For example, it is well known that mechanically weak materials, such as fractures, show non-linearly elastic, or inelastic, hardening behaviours even at small deformations (Bemer et al. 2001;Deude et al. 2002;Lemarchand et al. 2009;Bidgoli et al. 2013). We acknowledge the simplifying assumptions used in the current work, with a view to incorporating more realistic deformation behaviours on the basis of the modelling concepts developed herein.
To keep notation brief, we refer to the low permeability storage continuum as the matrix continuum, and the low storage, high permeability transport continuum as the fracture continuum. However, this work is sufficiently general such that other multiscale materials can be considered e.g. soil aggregates (under the assumption of infinitesimal deformations) (Choo and Borja 2015;Choo et al. 2016).

Volume averaging
We define the averaging operation, and assumptions therein, required in the homogenisation approach. A dual-continuum representation can be justified if an REV can be taken from a large macroscopic structure. Identification of an REV requires the satisfaction of the scale separation principle summarised as (Bear and Bachmat 2012), where s, S and L denote the characteristic lengths at the local heterogeneity, REV and macroscopic body scales respectively. Eq. (1) should honour length scale requirements for the physical system, both geometrically ( Fig. 1), and with respect to the wavelengths of the physical process (Auriault 2002;Geers et al. 2010). Accordingly, the REV represents the scale at which relationships between averaged quantities are defined (Fig. 1).
Defining an REV over fractured media is a subject of much debate due to the challenge of establishing criteria for scale separation (Long et al. 1982;Neuman 1988;Min and Jing 2003;Berre et al. 2019). In the following, however, we suppose a material for which an REV can be defined, such as densely fractured rock masses (Berkowitz 2002). To proceed we assume statistical homogeneity of the underlying material and thus make use of the volume average over the REV (Nemat-Nasser and Hori 1993), where z is an arbitrary tensor field, x is a position vector locating an REV within the macroscopic composite body, and |Ω| is the volume of an arbitrary REV within the body.

Preliminaries
To develop the homogenisation problem introduced by Dormieux et al. (2006) we consider the domain, Ω, over an REV in which there exists a porous matrix continuum, Ω m ⊂ Ω, and porous fracture continuum Ω f ⊂ Ω. We assume linear poroelasticity for each continuum (Biot 1941, Coussy 2004, and that the continua are saturated by the same slightly compressible fluid. Further, we assume isothermal evolutions and zero initial stress and pressure conditions. An important assumption is that we consider microscopic fluctuations in pressures are negligible with respect to the macroscopic (average) continuum pressures (Dormieux et al. 2006). As a result, fluids are assumed to be in steady state, but at different equilibrium pressures, within the respective continua in the REV. Accordingly, we model solid-fluid interactions at the microscale using uniform continuum pressures (Van den Eijnden et al. 2016).
With the given assumptions, the local constitutive model for a continuum, α, is then where C α [α = m, f ] is the intrinsic fourth-order stiffness tensor for continuum α, and the second-order tensors, σ α , α , b α , are the microscopic Cauchy stress and linearised strain tensors, and intrinsic Biot coefficient for continuum α respectively. Parameter n −1 α is the inverse of the Biot modulus, P α is the macroscopic fluid pressure, and dϕ α = ϕ α − ϕ 0 α is the evolution of the local Lagrangian porosity from the reference state (denoted by superscript 0), all written in terms of continuum α. The local Lagrangian porosity is the ratio of the continuum pore volume |Ω p α |, to the bulk volume of the undeformed continuum configuration, |Ω 0 α |. As is customary, we take stress and strain as positive in the tensile direction.
It is useful to re-write Eqs.
(3) to (4) in a unified way as follows (Dormieux et al. 2006), where C(x), and the prestress tensor distribution related to the fluid pressure (Chateau and Dormieux 2002), σ p (x), are given by respectively. The essence of the homogenisation approach is to define a boundary value problem on the REV, the solution to which allows for the determination of macroscopic constitutive properties. Accordingly the conservation of momentum boundary value problem is defined as where ∂Ω is the boundary of Ω, u is the microscopic displacement vector and E is the macroscopic (or surface prescribed) strain tensor. Quantities denoted byˆare boundary assigned values i.e. u =û on ∂Ω.
Using the averaging operation, Eq.
To give the link between microscopic fields, in this case strain, and macroscopic counterparts we consider a mapping between (x) and E. Owing to the linearity of Eq. (8) we can define a linear mapping so that where A(x) is the fourth-order mapping tensor (Hill 1963). Finally, combining Eqs. (11) to (13) it can be shown from which we can see where I is the fourth-order identity tensor.

Recovery of the constitutive system
From the superposition property in linear systems, Eqs. (8) to (10) can be decomposed into two subproblems. Subproblem I can be interpreted as a drained poroelastic problem: with where the macroscopic stress tensor, Σ = σ (Hashin 1972), and C * = C : A : E is the upscaled stiffness tensor for the dual-material. Subproblem II defines a constrained material, E = 0, subject to loading via the prestress field, σ p : with From Dormieux et al. (2006) one can show where Σ p = σ p : A. Eq. (24) is a part of a classical result in micromechanics referred to as Levin's theorem (Levin 1967). That is, the macroscopic constitutive equation follows the form of the linear local constitutive relation, Eq. (3), where we make use of the linearity of the problem to superpose subproblems I and II. Owing to the definition of C(x), and from Eqs. (14) to (15), the homogenised stiffness tensor of the composite dual-material is defined as Similarly, the homogenised prestress tensor is given as Intuitively, Eq. (27) can be interpreted as a weighted sum of the continuum pressures. In the work of Borja and Koliji (2009), the authors derive a pore fraction weighting formulation that is thermodynamically consistent. Such an approach was also proposed in Coussy (2004). Given the thermodynamic consistency, it would be interesting to see how one could recover a pore fraction weighted formulation within the general framework of microporomechanics. To proceed,using Eq. (27), and with the result from Eq. (26), we can identify the first of the macroscopic constitutive parameters, that is the effective Biot coefficients, From the energy approach to poromechanics (Coussy 2004), the dual-continuum model requires state equations for the evolutions of macroscopic Lagragian porosity (Ashworth and Doster 2019b). Accordingly, for subproblem I where we have made use of Eq. (4) in defining Eq. (30). Given subproblem II we have where we have used Eq. (11) together with the fact II = 0 to eliminate v f f . To advance we must substitute for v m II m . Following Dormieux et al. 2006, v m II m can be expressed as With Eq. (33) in Eqs. (31) to (32) we recover where the effective constitutive parameters N −1 α and Q −1 α are defined as Provided the storage continuum is isotropic, where 1 is the second-order identity tensor (Dormieux et al. 2006).
Finally, through superposition of subproblems I and II for the macroscopic variables Σ and dφ α , we recover the anisotropic, dual-stiffness constitutive model for the dual-scale, poroelastic material as where expressions for the effective constitutive parameters

Model equivalencies
Under certain conditions, the parameter models just referenced reduce to other mechanical property based parameter models proposed in literature. For example, in the case of soils, the high permeability (transport) continuum is all void space (Koliji et al. 2008). As a result C f = 0, and we recover the original anisotropic parameter models proposed by Dormieux et al. (2006). For an isotropic material, the constitutive system can be written in terms of scalar invariants of the tensorial quantities. Accordingly, with a change from the mixed-compliance constitutive formulation to a pure-stiffness formulation, we recover the dual-stiffness models introduced by Berryman (2002) (Ashworth and Doster 2019b). Ashworth and Doster (2019b), using parameter models by Berryman (2002), show some possible situations in which the inclusion of intrinsic fracture stiffness properties may be important. Further, in Kim et al. (2012), the authors show use of Berryman (2002) type coefficient models lead to well-posed mathematical problems, an important consideration for numerical modelling. Finally, combining the void space transport and isotropic dual-material assumptions, allows us to recover parameter models originally proposed by Berryman and Wang (1995) and Khalili and Valliappan (1996).
Under long-term drainage, P m = P f , dual-continuum models should reduce to single-continuum equivalents (Berryman and Wang 1995). As a result, we recover the following compatibility relations: where C s is the solid-grain stiffness tensor and φ = φ m + φ f . Eqs. (43) to (44) hold provided C s is the same for both the matrix and fracture continua. Accordingly, applying the long-term drainage condition, and contracting Eqs. (40) to (42) we recover the single-porosity constitutive model originally proposed in Biot (1941), albeit for anisotropic materials. Alternatively, we could recover the single-porosity model by setting v f = 0, and thus C f = 0 with C * = C m .

Numerical framework
Here we introduce the computational framework used for modelling the coupled DC problem. We start by introducing the strong form of the DC poroelastic problem and then progressing to its fully discrete counterpart.

Strong form
In addition to Eqs. (40) to (42), we require constitutive relations for intra-and inter-continuum mass flux terms, w α and γ α respectively. Intra-continuum mass flux is given according to Darcy's law, where q α is the volumetric flux vector associated with continuum α, ρ l and µ l are the intrinsic fluid density and fluid viscosity respectively, g is the gravity vector and k * α is the macroscopic continuum permeability tensor. The inter-continuum mass flux, γ α , is given according to a first-order transfer term originally proposed by Warren and Root (1963), where k denotes the interface permeability, taken here as the intrinsic matrix permeability (Barenblatt et al. 1960;Choo and Borja 2015), and κ is a parameter referred to as the shape factor (Warren and Root 1963). In this work we use an analytically derived κ for an isotropic matrix given according to Lim and Aziz (1995), where N is a dimension parameter related to the number of fracture sets, and s is the characteristic spacing length of the fracture continuum ( Fig. 1). Finally, we give the compatibility between the macroscopic linearised strain tensor and macroscopic displacement, U , as where we introduce notation ∇ to denote the symmetric gradient operator on U .
The conservation equations considered for the DC poroelastic problem are the momentum equation and the continuity equations for each continuum Notations ρ andγ in Eq. (49) are the bulk density of the dual-material, and a momentum source arising from the inter-continuum mass transfer respectively. For the remainder we assumeγ ≈ 0, with respect to the other force density terms in Eq. (49). Notation m l,α in Eqs. (50) to (51) is the fluid mass content associated with continuum α. The fluid mass content is given by m l,α = ρ l φ α . We consider the conservation equations over a domain, Ω D ⊂ R 2 , bounded by ∂Ω D . The domain boundary is separated into disjoint boundary segments corresponding to Dirichlet and Neumann boundary conditions for the mechanical and flow problems. For the mechanical problem this implies displacement (Γ U ) and traction (Γ T ) boundary conditions. To ensure well-posedness Γ U ∪ Γ T = ∂Ω D and Γ U ∩ Γ T = ∅. For the flow problem the boundary conditions for a given continuum are pressure (Γ P α ) and flux (Γ Q α ). Again, for a well-posed problem we have The strong form is finally defined as: Find U , P m and P f that satisfy Eqs. (49) to (51) subject to boundary conditions: with initial conditions for all (X, t) ∈ (Ω D × t = 0). Notation X is the macroscopic position vector. The single porosity linear poroelastic model can be recovered from Eqs. (49) to (58) under the assumption P = P m = P f and combining Eq. (50) and Eq. (51). With the contraction to a single continuum system, DC constitutive parameters reduce to single porosity equivalents, Eqs. (43) to (44).

Weak form
The weak formulation of the strong form introduced previously requires the definition of the appropriate function spaces. Accordingly, solution spaces for continuum pressure and the displacements are S Pα = L 2 (Ω D ) and S U = {U ∈ H 1 (Ω D ) d : U =Û on Γ U } respectively, where L 2 and H 1 are the typical square integrable and first-order Sobolev function spaces. Weighting function spaces are then defined as W Pα = L 2 (Ω D ) and To progress we substitute the constitutive equations, Eqs. (40) to (42) and Eqs. (45) to (46), and macroscopic strain compatibility relation, Eq. (48), into Eqs. (49) to (51). We adopt the material and fluid assumptions introduced in Section 2, whilst also neglecting gravitational effects. Assuming isotropic matrix material results in Q −1 m = Q −1 f = Q −1 and B m = B m 1. Further, we restrict our anisotropic experiments to orthotropic materials. Finally, comparing trial functions against weight functions, the weak form is defined as: The bilinear form g(·, ·) in Eq. (59) is given by where Σ (U ) = C * : ∇U is the effective stress tensor. The term M −1 α in Eqs. (60) to (61) is given as where K l is the fluid bulk modulus.

Discrete block matrix form
The discrete counterpart to Eqs. (59) to (61) Ashworth and Doster (2019a) to dual-continuum materials. We use this hybrid framework in the current work due to its availability. However, recent works have shown the current modelling framework to be suitable for subsurface applications where complex geometrical structures can lead to irregular grids not easily handled by standard finite-element methods (Andersen et al. 2017b;Coulet et al. 2019). We partition our domain into disjoint elements (or cells). Accordingly, for the DC problem Ω D = ∪ n elem j=1 Ω D j , where n elem is the number of elements. Notation Ω D j denotes the dual-continuum element for which there are two pressure degrees of freedom, corresponding to each continuum.
We define the following discrete solution spaces for the DC problem as S h Pα ⊂ S Pα and S h U ⊂ S U . Discrete weighting spaces are given as W h Pα ⊂ W Pα and W h U ⊂ W U . Discrete continuum pressure fields, P h α ∈ S h Pα , and discrete displacement fields, U h ∈ S h U , are given according to the following interpolation relations respectively, where n node denotes the total number of vertices, andP j α andŨ b are pressure and displacement degrees of freedom respectively, with the corresponding basis functions denoted by I j and N b .
In FVM we considerP j α to be cell-centred quantities. Notation I j is then an indicator function for continuum α given as Further, we replace discrete pressure weight functions, ω h α ∈ W h Pα , by the indicator function whilst also using Eq. (64) such that Eqs. (60) to (61) can be interpreted as element-wise conservation statements. Using Gauss's theorem, element-wise divergence of flux volume integrals in Eqs. (60) to (61), are turned into face-wise surface integrals. In this work we use the two-point flux approximation to calculate these face-wise flux integrals (see Lie (2019) for further details).
The nodal basis function matrix, N b , takes the identity matrix 1 when located at node b and 0 at all other nodes. VEM is a Galerkin based method, thus the discrete displacement weight, η h ∈ W h U , is an interpolation of the type shown in Eq. (65). However in VEM, contrary to standard finite-element methods, the bilinear form with discrete fields can never be directly calculated as basis functions are never explicitly defined. Due to basis function independence, VEM can be interpreted as a generalisation of the finite-element method to arbitrary polygonal and polyhedral meshes. Such a property is desirable for subsurface modelling, where degenerate cells and hanging nodes are encountered (Andersen et al. 2017b). Instead, the idea in VEM is to approximate the bilinear form, such that where g h (η h , U h ) = n elem j=1 g h j (η h , U h ), and where details of the element-wise firstorder bilinear VEM approximation, g h j (η h , U h ), can be found in Gain et al. (2014) and Andersen et al. (2017b). Finally, as part of the VEM assembly, the constitutive relation, Eq. (62), need only be computed once, similar to a one-point quadrature finite-element scheme (Da Veiga et al. 2015).
Replacing solutions and weighting functions with their discrete counterparts, and using the time discretisation, the discrete residual equations from Eqs. (59) to (60) are where we make use of Voigt notation for tensor representation. Notation ∆z n+1 = z n+1 − z n , where n denotes the current time level. Details of the VEM calculations for the boundary and gradient terms involving N a in Eq. (68) can be found in Andersen et al. (2017b). Even though we assume linearity in the current work, poroelastic problems are generally non-linear due to material and geometric non-linearities. To provide a general numerical framework we therefore present the discrete equations describing the DC problem following application of Newton's method. In MRST this is handled naturally using an automatic differentation framework to generate the Jacobian. We give the discrete system of equations in block matrix form as where ] . The individual matrices comprising the Jacobian in Eq. (71) are given as where X i denotes the centroid of element Ω D i , and G ij,α is the transmissibility matrix for continuum α arising from the two-point flux approximation (Lie 2019).
Finally, Eq. (71) is solved using a fully coupled approach (Lewis and Schrefler 1998), although extensions to sequential solution strategies for DC materials have been shown in Kim et al. (2012) and Ashworth and Doster (2019a).

Numerical tests
With the framework in-hand, we present and conduct the numerical tests used to investigate whether the macroscopic dual-continuum poroelastic model can capture global flow and deformation behaviours of a fine-scale (FS) explicit model. In doing, we review several considerations for the interpretation of the results at the scale of the DC model.

Test cases
We introduce four numerical experiments to test the validity of the DC poroelastic concept. In each case we consider an idealised representation of a naturally frac-tured rock sample. Our idealisation comes in that we assume the fracture fabric to be periodic. To start we consider an undeformable isotropic material to understand the physics of the flow problem. We progress by introducing mechanics to the isotropic system, and then adding complexity by considering anisotropic material cases.
In every case we consider the dimension of the domain to be 1 m × 1 m. Each experiment then represents a thin 2D slice taken from a 3D sample such that, in the case of the mechanical problem, the plane-stress assumption applies.

Undeformable isotropic
For this test we study an (isotropic) undeformable matrix permeated by an isotropic undeformable fracture network. The test is setup as a uniaxial drainage problem, such that the top boundary is open to flow,P m =P f = 0, whilst the left, right and bottom boundaries are zero flux boundaries (Fig. 2a). Initial pressures for the continua are set at P 0 m = P 0 f = 2 MPa. Volume fractions for matrix and fracture material are v m = 0.998 and v f = 0.002 respectively, given a fracture spacing, s, of 0.1 m. Local porosities for the two continua are then prescribed as ϕ m = 0.1 and ϕ f = 0.9, where the volume fractions link the global and local Lagrangian porosities so that φ α = v α ϕ α . Intrinsic matrix permeability, k m , is taken as 0.01 md, whilst individual fracture permeability, k f , is calculated using the parallel plate model with a fracture aperture of a f ≈ 1.05 × 10 −4 m (Witherspoon et al. 1980). The resulting permeability is 950 d for each fracture 1 . Fluid properties are ρ l = 1000 kgm −3 , µ l = 1 cp, and K l = 2.5 GPa. Upscaling individual fracture permeability to a continuum permeability for use in the DC model is done using the cubic law (Witherspoon et al. 1980). The resulting isotropic fracture continuum permeabitity is k * f ≈ 1000 md. Finally, due to the dissociation by the fracture network, the macroscopic matrix permeability is zero.

Deformable isotropic
We now consider a deformable counterpart to the experiment described in Section 4.1.1 (Fig. 2b). Accordingly, the Young's moduli for the matrix and fracture materials are E m = 36 GPa and E f = 36 MPa respectively. The latter is chosen for illustrative purposes as E f = E m /1000. Both continua are assigned a Poisson's ratio of ν = 0.2. For an isotropic medium under the plane-stress assumption, the stiffness tensor written with Voigt notation is given as where parameter G = E/(2(1 + ν)) is the shear modulus. Entries for C m and C f can be calculated with Eq. (76) and the defined intrinsic parameter values. For C * , parameters must be calculated by homogenisation. In Ashworth and Doster (2019b) the author's suggest using the Hashin-Shtrikman lower bounds (Hashin and Shtrikman 1963), as an initial homogenisation approach for the estimation of the mechanical properties of densely fractured rock. For the bulk and shear moduli these lower bounds are quoted as (Hashin and Shtrikman 1963), where K α and G α are the 3D bulk and shear moduli for continuum α respectively.
We map between the 3D bulk modulus calculated in Eq. (77) and the 2D homogenised bulk modulus under plane-stress, K * , using the following relation (e.g. Torquato 2002), The Poissons ratio for the composite dual-continuum under plane-stress is given by Finally, the homogenised Young's modulus, E * , can be recovered as With Eqs. (80) to (81) the homogenised parameters are ν * = 0.2 and E * = 18.0 GPa. We assume the matrix and fracture skeletons to be made up of the same solid material. We then assign a solid modulus, K s , of 70 GPa for both continua.
For the coupled mechanics and flow problem we consider a different method of initialisation to Section 4.1.1. Instead of assigning initial continuum pressures, we define the starting point for the experiment to be the undrained, loaded configuration, P 0+ m , P 0+ f , U 0+ . This undrained state is induced following the application of an instantaneous load on an unpressurised and undeformed domain, P 0 m = P 0 f = 0 and U 0 = 0. Loading is prescribed as a vertical traction of −Σ · n y = −2 MPa on the top boundary. The domain is horizontally constrained at the boundaries, but remains free to move along the vertical axis apart from at the bottom boundary where the sample is fixed. The parameters for flow are as defined in Section 4.1.1

Geometry-induced anisotropy: explicit computation of C *
The third experiment is concerned with an anisotropic deformable material. Anisotropy has recently been studied in poroelastic DC materials in the context of flow properties (Zhang et al. 2019). However, here we consider the directional dependence of both mechanical and flow properties. Anisotropy is introduced geometrically by considering just a single vertical fracture set which is aligned with the second principal axis (Fig. 3a). The 2D domain is then orthotropic. Whilst anisotropy exists at the macroscale, the intrinsic mechanical parameters remain isotropic for each continuum and are as described in Section 4.1.2. The plane-stress stiffness tensor for an orthotropic material is given by Parameters of the homogenised stiffness tensor may be approximated explicitly for this geometry, using mixture theory. Accordingly, for the Young's moduli where having removed a fracture set, the volume fraction of the fracture continuum is now v f = 0.001 (resp. v f = 0.999). For the homogenised Poisson's ratio, ν * 21 , and shear modulus, G * 12 , mixture theory gives The other Poisson's ratio, ν * 12 , is readily determined by the symmetry in Eq. (82) which requires ν 12 E 2 = ν 21 E 1 . From Eqs. (83) to (84) and the aformentioned symmetry relation, the mechanical parameters are given as E * 1 = 18.0 GPa, E * 2 = 36.0 GPa, ν * 21 = 0.200, ν * 12 = 0.100 and G * 12 = 7.50 GPa. The anisotropic fracture continuum leads to an anistropic permeability tensor so that permeability in the x and y directions are k * f,x = 0 and k * f,y ≈ 1000 md respectively. For the matrix, the macroscopic permeability is also anisotropic with k * m,x = 0 and k * m,y ≈ 0.01 md. The remaining flow parameters, boundary conditions and initialisation are as described in Sections 4.1.1 and 4.1.2.

Material-induced anisotropy: numerical computation of C *
The final experiment is on an anisotropic material with two fracture sets aligned with each of the principal axes (Fig. 3b). Anisotropy is now introduced through the fracture material, with each fracture set having different intrinsic mechanical and flow properties. These property differences are in analogy to fractures containing different amounts of infill material. To represent this conceptually within the model we assign different intrinsic porosities to the individual fracture sets. Further, we separate the intrinsic Young's moduli and permeabilities of each fracture set by two orders of magnitude. For the horizontal fracture set we assign ϕ h f = 0.9,  9.5 d. Upscaling the fracture permeability remains trivial, with k * f,x ≈ 1000 md and k * f,y ≈ 10 md. However, homogenistion for the parameters in the homogenised stiffness tensor now cannot be done by explicit approximation. Instead we use a deformation-driven computational homogenisation approach: We generate unit strains for a sequence of linear displacement boundary conditions, and in doing, determine the entries of C * (Daniel et al. 1994). Linear displacements are chosen as they produce better estimates for effective stiffness tensors for materials with a stiff matrix and weaker inclusion material (Pecullan et al. 1999), as is the case here.
With the computational homogenisation approach, the mechanical parameters in C * are calculated as E * 1 = 32.7 GPa, E * 2 = 3.40 GPa, ν * 21 = 0.019, ν * 12 = 0.173 and G * 12 = 1.28 GPa. The overall volume fraction for the fracture continuum is the same as in experiment two. However, the intrinsic Lagrangian fracture porosity is now the arithmetic average of the two intrinsic fracture set porosities (ϕ f = 0.65). Fluid and matrix properties remain the same as those for the other experiments. Boundary conditions and initialisation are the same as in experiments two and three.

Modelling considerations
Here we review several considerations to enable the interpretation of the test results to follow.

On the REV
Our periodic assumption of the underlying microstructure eases the requirements on our definition for an REV. In this periodic case, all the necessary geometrical and physical process information is captured within an elementary cell that is the size of the heterogeneity (s) (Royer et al. 2002;Dormieux et al. 2006). The separation of scales is now defined as s L. The elementary cell definition of our REV will be useful for interpreting the discretisation choice of the DC problem.

Meshing
For the four tests we discretise the fine-scale problem with a 200 × 200 Cartesian mesh, that is locally refined around the fractures (Fig. 4). For the dual-continuum problem we discretise the domain using a 10 × 10 Cartesian mesh. In the latter case, each element then coincides with an elementary cell (in the geometrical sense) (Fig. 4).  Fine-scale and dual-continuum fields are compared at an observation point at the base of our samples. At this point we assume that our pressure solutions are sufficiently smooth, thus satisfying the physical process scale separation requirement.
Further, the observation point coincides with the macroscopic material point, in this case the element centroid (Fig. 4).

Quantities of interest
At each test observation point we consider the element-wise total and continuum volumetric strains, and continuum pressures. Element-wise total and continuum volumetric strains are defined as E v j = tr(E j ) = ∆|Ω D j |/|Ω D,0 j |, and E v j,α = ∆|Ω D j,α |/|Ω D,0 j,α |, respectively. We compare averaged results over the fine-scale to element level results from the dual-continuum. To enable the former, we define the following discrete continuum counterpart to Eq. (2) where z is a scalar field of interest. Continuum averaged pressures and volumetric strains can then be recovered using Eq. (85) with discrete microscopic fields p h or v,h in place of z h . Total volumetric strain is likewise obtained using Eq. (85) by replacing |Ω D j,α | with |Ω D j |. For the DC problem, pressures and element-wise total volumetric strain are recovered naturally from the element centroid (Andersen et al. 2017a). To get continuum strains, however, we must take a different approach. Starting with the matrix continuum, and comparing a volume averaged change in local porosity, Eq. (4), to the effective change in matrix porosity given by Eq. (41), such that dφ m = v m dϕ m , allows us to derive the following expression for the volumetric matrix strain We note the expression for E v m in Eq. (86) is only possible for an isotropic matrix as the inverse contraction map involving b m is otherwise ill-posed. With E v and E v m we can recover the fracture volumetric strain, E v f , for the DC model using Eq. (14).

Results and discussion
Here we present the results and analyses for the numerical test cases described in Section 4.1 under the modelling considerations described above. All results are given from observation points such as that shown in Fig. 4.

Undeformable isotropic
Fig . 5 shows the element averaged pressure evolutions from both fine-scale explicit and dual-continuum simulations for the undeformable isotropic material case. Both models show a rapid decrease in fracture pressure within the first millisecond followed by a delayed pressure response in the matrix. These general patterns can be attributed to the contrast in continuum permeabilities. Whilst both models show general decreasing trends, the FS fracture pressure decrease begins to smooth out at lower pressures. Further, the onset of fine-scale matrix pressure diffusion happens earlier. When matrix pressure diffusion does occur in the DC model, the process occurs more rapidly (indicated by a steeper gradient) than in the fine-scale case. FSm DCm FS f DC f Figure 5. Matrix and fracture continuum element averaged pressure evolutions for the undeformable isotropic test case. 'FS α ' and 'DC α ' denote quantities related to fine-scale and dual-continuum models for continuum α respectively.
The disparities in matrix and fracture pressure diffusion between the two modelling approaches arise from the first-order transfer term, Eq. (46), used by the dual-continuum model. In using a linear mass transfer model one implicitly places a pseudosteady state diffusion assumption on the communication between matrix and fracture continua. As a result, transient matrix drainage effects are neglected by the DC approach. Neglecting transient effects leads to the delay in DC matrix pressure diffusion we see in Fig. 5, and the loss of pressure support in the fractures.
Shortcomings of using simplified transfer concepts have been well documented in literature (e.g. Berkowitz et al. 1988;Lemonnier and Bourbiaux 2010). Previous works have thus sought to improve on the linear inter-continuum flow coupling term by including transient effects (Zimmerman et al. 1993;Sarma et al. 2004;March et al. 2016;Zhou et al. 2017). However, in the current work, we acknowledge the shortcomings of the transfer term used herein, with the focus being on understanding the coupled poroelastic behaviour.

Deformable isotropic
Pressure and total element volumetric strain results for the deformable isotropic case are shown in Fig. 6. In Fig. 6a both modelling approaches predict higher induced initial pressures in the fracture than in the matrix. Further, both approaches show rapid decreases in fracture pressure and gradual decreases in matrix pressure. In Fig. 6b the FS and DC models show increasing volumetric strain evolution behaviours. However, for both pressure and strain, specific differences of the variable fields between the two modelling approaches may be observed at both early and late times. The disparity in late-time matrix pressure evolution (Fig. 6a), is again due to the first-order inter-continuum transfer model (see Section 4.3.1). Over the same late-time period, we also observe a difference in volumetric strain (Fig. 6b). Of more interest in Fig. 6 are the early-time results for continuum pressures. In Fig. 6a both fine-scale matrix and fracture pressures exhibit non-monotonic behaviour, known as the Mandel-Cryer effect. These pressure rises are not seen in the dual-continuum pressure responses. A similar observation was also made in the work of Zhang et al. (2019), albeit for a different problem. Fig. 7 shows individual continuum volumetric strain evolutions. In Fig. 7b both modelling approaches show similar increasing strain behaviour with time. However in Fig. 7a, we observe the FS matrix strain shows early-time non-monotonic behaviour, contrary to DC matrix strain. The early non-monotonic differences in pressure and matrix strain between the two modelling approaches result from the underlying pressure assumption made for the DC model. In the derivation of the constitutive model in Section 2.2 we assumed a local equilibrium of pressure within each continuum over an REV. The induced response predicted in by the DC model is thus for a system in mechanical and hydrostatic equilibrium. Instead, the fine-scale model makes no such pressure assumption. To understand the specific impacts of the latter it is interesting to look at the local flow and deformation behaviours shown by the fine-scale. Fig. 8 shows the FS pressure and volumetric strain responses within the first 100 microseconds. At t 0+ in Fig. 8a, we observe pressure in the horizontal fracture is higher than the vertical fracture. This disequilibrium is concurrent with the negative and positive fracture strains for the horizontal and vertical fractures respectively (Fig. 8b). From t 0+ to t 1 , Fig. 8a shows, away from the fracture intersection, horizontal fracture pressure drops slightly. However, vertical fracture pressure increases. These pressure changes occur with further contraction and expansion respectively (Fig. 8b). Over the same time period matrix strain increases ( Fig. 8c). From t 1 to t 5 the pressure in both fractures is increasing (Fig. 8a), with matrix and fracture deformations following the same evolution paths described previously. Finally, at t 10 the fractures reach a pressure equilibrium (Fig. 8a).
(b) Volumetric strain (colour scale for fracture strain highlighted) Figure 8. Pressure (a) and volumetric strain (highlighted for fractures and matrix, (b) and (c) respectively) fields at different time levels, t i , for the FS representation of the deformable isotropic material. Each field plot is 5 mm×5 mm and is located at the observation point. Subscript 0+ denotes the time level corresponding to the undrained, loaded configuration.
We can now explain the early-time non-monotonic behaviours in Fig. 6a and Fig. 7a with the description of the local processes shown in Fig. 8. Following t 0+ , intra-fracture flow is driven by the pressure disequilibrium between the horizontal and vertical fractures. Between t 0+ and t 1 horizontal fracture contraction occurs primarily due to the dissipation of the fluid pressure support. Vertical fracture expansion follows due to poroelastic coupling to accommodate incoming fluid from the horizontal fracture. As the vertical fracture expands it forces the contraction of the matrix, and thus the increasing matrix strain shown in Fig. 7a and Fig. 8. After t 1 , deformation drives the horizontal fracture pressure increase due to fluid compressibility. The overall fracture continuum pressure increases (Fig. 6a), with strain generating pressure in the horizontal fracture, whilst the pressure change associated with vertical fracture expansion slows. The latter occurs due to the low matrix permeability which prohibits dissipation of excess matrix pressure, until later times. As a result, the undrained matrix stiffness increases with its progressive contraction, slowing vertical fracture expansion until a mechanical equilibrium is reached. The overall fracture continuum pressure rise finally stops when the fractures have reached mechanical equilibrium with the matrix and an internal hydrostatic equilibrium.
The local processes shown by the FS model are not captured by the DC model due to the underlying homogenisation assumptions made in the latter. However, Figs. 6 to 7 do show that, aside from the local equilibration processes, the DC model can capture the global poroelastic behaviours of the FS model.

Geometry-induced anisotropy: explicit computation of C *
Here we show the results for the geometry-induced (single-fracture set) anisotropy case. Pressure and total volumetric strain are given in Fig. 9, whilst Fig. 10 shows the individual continuum volumetric strains. Fig. 9a now shows a smaller disparity between the initial matrix and fracture pressures, with the fracture pressure being only slightly higher. Further, we do not observe the Mandel-Cryer effect in the FS model. However, away from the initial pressures, the general trends we see in Fig. 6a can still be observed in Fig. 9a. Specifically, a rapid decrease in fracture pressure followed by a smoother matrix pressure decrease. As expected, the late-time differences in matrix pressure observed previously are present in the current test. For both modelling approaches there is a good agreement in matrix and fracture pressure evolutions. The total element volumetric strain evolutions are also similar between the two modelling approaches, with an overall increase in strain as the material compacts.
The similarity in total volumetric strain between the two approaches is reflected in the individual continuum strains (Fig. 10). The matrix shows early-time expansion behaviour followed by contraction. Fracture deformation is coupled to matrix deformation (and vice versa). Fracture contraction is therefore followed by a period of expansion as the matrix drains and contracts.
The small difference in initial continuum pressures observed in Fig. 9a can be . Matrix and fracture continuum element averaged pressure (a) and total element volumetric strain (b) evolutions for the (deformable) anisotropic test case with one (vertical) fracture set. 'FS α ' and 'DC α ' denote quantities related to fine-scale and dual-continuum models for continuum α respectively. explained by considering the geometric anisotropy induced by the fractures. With the fracture set being aligned with the direction of loading (Fig. 3a), the stiffer matrix acts like a series of columns, supporting a significant portion of the applied load. Through the coupling between stress and pressure, the low portion of stress 'seen' by the fracture phase leads to the low induced fracture pressure shown in Fig. 9a. Finally, the absence of the Mandel-Cryer effect in the current case is due to the pressure being at equilibrium within the single fracture set. As a result local processes do not drive early-time poroelastic intra-and inter-continuum pressure generation.

4.3.4
Material-induced anisotropy: numerical computation of C * Fig. 11 and Fig. 12 show pressure and total strain, and individual continuum strain results respectively, for the material-induced anisotropy case. Both models in Fig. 11a show a strong difference in the early-time magnitudes of matrix and fracture pressures. The FS model shows similar early-time Mandel-Cryer fracture behaviour to what we observed in Fig. 6a. In contrast, the FS matrix nonmonotinicity is negligible in Fig. 11a compared to the isotropic case. At later times we see a significant non-montonic evolution in matrix pressure that is shown by both modelling approaches. This non-monotonic matrix pressure rise starts earlier in the FS model than the DC model. Finally, we observe that matrix and fracture diffusion starts at similar times, indicating a single-continuum response. Coupled to the delayed fracture diffusion response is the delayed increase in total volumetric strain (Fig. 11b). In Fig. 12 both modelling approaches give similar continuum strain evolutions. Similar to Fig. 7a, Fig. 12a shows the DC approach misses the early-time matrix strain non-monotinicity displayed by the FS approach. However, contrast to Fig. 7a, Fig. 12a shows a smoother early-time FS matrix strain non-monotinicity, whilst the late-time matrix strain non-monotinicity for both approaches is much sharper.
Results in Figs. 11 to 12 can be explained by considering the material anisotropy in the fracture continuum. The smoother early-time non-monotinicity in FS matrix strain occurs because the vertical fractures are stiffer. These fractures then expand less with incoming fluid, reducing poroelastic coupling (and thus deformation) with the matrix compared to the isotropic case. As a result, since FS matrix pressure does not change significantly, the initial matrix pressures for the two modelling approaches are similar. This result suggests mechanical anisotropy can noticeably affect the degree of inter-continuum coupling. The delay in fracture pressure diffusion occurs due to the low vertical fracture permeability. Accordingly, we see the non-monotonic rise in matrix pressure with local inter-continuum equilibration . Matrix and fracture continuum element averaged pressure (a) and total element volumetric strain (b) evolutions for the (deformable) anisotropic test case with two orthogonal fracture sets. 'FS α ' and 'DC α ' denote quantities related to fine-scale and dual-continuum models for continuum α respectively.
processes occuring at similar timescales to macroscopic fracture flow. Further, the pseudosteady state mass transfer assumption leads to the delayed response of this non-monotinicity in the DC model. The influx of fluid from the fractures into the matrix is accompanied by expansion of the matrix material, followed by contraction as fluid drains out (Fig. 12a).
The results in the current test show again how the DC model misses earlytime effects due to local equilibration processes. Neglecting these local processes is implicit due to the steady state pressure assumption made during homogenisation. However, once local equilibration is reached, the DC model does predict the general poroleastic behaviours of the FS model.

Conclusions
Dual-continuum models are an implicit approach to modelling multiscale materials. Further, with the appropriate extensions, they can be used to model complex multiphysics problems such as the coupled mechanics and flow phenomena studied in this work.
In this paper we derived a dual-continuum poroelastic constitutive model that makes no assumptions on the material symmetry and mechanical properties of the dual-material and its constituents. We termed the resulting model the anisotropic, dual-stiffness constitutive model. Further, we discussed how under isotropy of the continua and bulk material, and void space assumptions of the high permeability transport phase, previously introduced constitutive models can be recovered from the constitutive model developed herein. Secondly, using numerical modelling, we investigated whether the dual-continuum approach with the derived constitutive model, is able to capture the global poroelastic behaviours of fine-scale explicit models. We introduced the computational framework used to carry out our investigation and described the resulting numerical tests therein. We observed that anisotropy can have measurable impacts on flow and deformation behaviour. However, we showed that the DC approach is capable of capturing the global poroelastic behaviours for both isotropic and anisotropic FS equivalents. Discrepancies between the two model representations arise when local equilibration processes not accounted for in the homogenisation approach, are significant.
Finally, interesting extensions to the current work involve the study of nonlinear poromechanical effects, and measurement methodologies for the material parameters used herein. In the former, it is well known that fracture (and soil aggregate) deformation is geometrically non-linear leading to coupled material nonlinearities at the macroscopic scale. Modelling these scale effects requires compre-hensive multiscale modelling approaches, and is an active area of research Wang and Sun 2018;Castelletto et al. 2019). Lastly, in analogy to the work of Biot and Willis (1957), it highly desirable to develop methods of measurement for the parameters introduced in this work. In particular, the challenge remains how to map individual fracture characteristics to continuum properties. In this context, a microporomechanics framework could provide useful insights into experimental, and theoretical methodologies (e.g. Lemarchand et al. 2009).