Foundations and their practical implications for the constitutive coefficients of poromechanical dual-continuum models

A dual-continuum model can offer a practical approach to understanding first order behaviours of poromechanically coupled multiscale systems. To close the governing equations, constitutive equations with models to calculate effective constitutive coefficients are required. Several coefficient models have been proposed within the literature. However, an assessment of how these models relate to one another, and the underlying assumptions used in their construction, is still missing. To address this we first compare and contrast the dominant models existing within the literature. In terms of the constitutive relations themselves, early relations were indirectly postulated that implicitly neglected the effect of the mechanical interaction arising between continuum pressures. Further, recent users of complete constitutive systems that include inter-continuum pressure coupling have explicitly neglected these couplings as a means of providing direct relations between composite and constituent properties, and to simplify coefficient models. We analyse the impact of implicit and explicit decoupling assumptions using both theoretical and qualitative means. The former can be equated to the removal of a pressure source term in the governing equations. Analysis of the latter is conducted within the framework of micromechanics. Depending on the formulation, explicit decoupling assumptions are equivalent to conditions of isostress or isostrain. These conditions represent end-member states that correspond to the Reuss and Voigt bounds for composite moduli respectively. Use of these bounds goes on to directly affect coefficient models themselves. We conclude by offering recommendations for when to use different coefficient models, in addition to upscaling recommendations for when composite moduli are unavailable.


Introduction
Many natural and manufactured geomaterials exhibit strong heterogeneities in their material properties owing to the existence of porous constituents at various length scales. Examples of multiscale systems that are commonly encountered in subsurface operations include fissured or fractured rock and soil aggregates (Warren and Root 1963;Kazemi et al. 1976;Nelson 2001;Gerke 2006;Koliji (2008); Romero et al. 2011). Modelling of such materials is invaluable in understanding how these systems behave in response to extraneous activities. In general, modelling can be done using either explicit (e.g. discrete fracture matrix models) or implicit methods (e.g. continuum approaches) (Berre et al. 2018).
With respect to fractured systems, using explicit methods can be computationally prohibitive at large scales (Karimi-Fard et al. 2006;Gong 2007;Garipov et al. 2016). Additionally, explicit methods may require data (e.g. spatial data) that is not obtainable without direct access (Berkowitz 2002;Blessent et al. 2014). In cases where field-scale modelling of multiscale systems is required, implicit representations are then often preferred. The most common type of implicit model is the dual-continuum (or double-porosity) model, originally attributed to Barenblatt et al. (1960). In this, the dual-material is considered as the superposition of two overlapping continua, which communicate through a mass transfer term. Continua are defined on the basis of their material properties. For example, fractures (or inter-aggregate pores) generally have high permeabilities and poor storage capabilities, vice-versa the matrix. Although less detailed than their explicit counterparts, dual-continuum models can provide practical and valuable insight into the first order behaviours of multiscale systems. Further, use of these models is desirable due to the low number of fitting parameters that allow for efficient calibration to historical data.
Multiscale systems can also exhibit strong coupling between deformation and fluid flow, and vice-versa. This is known as poromechanical coupling (Rutqvist and Stephansson 2003), and is described by the well established poromechanical theory (see for example Biot 1941Biot , 1977Detournay and Cheng 1995;Coussy 1995Coussy , 2004Wang 2000;De Boer 2012;Cheng 2016). Aifantis (1977Aifantis ( , 1979 and Wilson and Aifantis (1982) were the first to introduce the generalised notion of deformation within the dual-continuum setting. Further offerings then came from Elsworth and Bai (1992), Lewis and Ghafouri (1997) and Bai et al. (1999). However, all of these models implicitly neglected the effects of coupling between pressures of different pore domains due to their postulation of the form of the constitutive equations. The absence of these pressure couplings was shown to give unphysical responses by Khalili (2003). Specifically, the author's results showed discontinuous pressure jumps in the matrix and fracture continua that were incompatible with the prescribed boundary conditions. The cause of the observations made by Khalili (2003) still remains an open question.
Additional models in which the constitutive equations included inter-continuum pressure coupling were introduced by Berryman and Wang (1995), Tuncay and Corapcioglu (1996), Loret and Rizzi (1999), Berryman and Pride (2002), Berryman (2002), Khalili and Valliappan (1996) and Khalili (2008). The difference between these presentations comes in the way that the authors choose to calculate the constitutive coefficients that govern a dual-continuum's poromechanical behaviour. For example, some authors implicitly assume the high permeability continuum to be all void space (e.g. Khalili and Valliappan 1996), whilst others allow for an intrinsic phase stiffness (e.g. Berryman 2002). Models (referred to as coefficient models herein) used to calculate the constitutive coefficients are required due to the potential difficulty in measuring these properties experimentally. Whilst various coefficient models exist within literature there is still no general guideline for how and when to use them.
More recent users of these later constitutive/coefficient models have explicitly decoupled pore domain pressures when expressing the constitutive relations in terms of stress and continuum pressures (pure stiffness setting) (see for example Nguyen and Abousleiman 2010;Kim et al. 2012;Mehrabian and Abousleiman 2014;Mehrabian 2018). This has been done as a form of non-algebraic closure and to provide explicit relations between composite and constituent properties resulting in simplified coefficient models. However, such decoupling assumptions have been made without discussing the origin and sensitivities that may arise as a result.
The first of aim this paper is to compare and contrast the various coefficient models existing within the literature. Our second aim is to provide a thorough investigation into decoupling assumptions, and how they may go on to affect a poromechanical system. Fundamental to our approach is treating the dual-medium as a composite. This allows us to analyse coefficient models, and certain intercontinuum coupling assumptions made on these models, within the context of micromechanics. Throughout this paper our reference multiscale material is that of a naturally fractured system. Such systems can be considered as void space inclusion composites or stiff inclusion composites (Fig. 1).
We structure the paper as follows. In Section 2 we introduce the governing and constitutive equations pertinent to double-porosity materials. For the latter set of equations we support their form using arguments from the energy approach to poromechanics (Coussy 1995(Coussy , 2004. Section 3 presents the most prevalent modelling approaches for calculating the effective poromechanical coefficients. Section 4 details the origins of explicit assumptions made on constitutive/coefficient models within the framework of micromechanics. From here we offer upscaling recommendations for constituent moduli when composite moduli may not be available. In Section 5 we use analytical solutions to the double-porosity Mandel problem to explore the physical implications, and relevance, of different coefficient models and decoupling assumptions. We conclude by way of offering recommendations for how to use coefficient models in light of (a) pressure coupling between pore domains, and (b) when composite moduli are not available and/or (c) intrinsic fracture stiffness effects need to be accounted for.
We note that work has been done on the determination of effective properties of multiple-porosity materials via homogenisation methods (e.g. Berryman 2006 andLevin et al. 2012). However, equivalent continuum models can fail to provide insight into processes occurring at the different porosity scales due to use of an averaged flow field (Berre et al. 2018). In contrast, this work is concerned with double-porosity materials for which two distinct flow fields exist. Upscaling of such flow fields for inelastic materials has been addressed by periodic homogenisation (Arbogast and Douglas 1990), but such a treatment for deformable materials is, to the best of the authors knowledge, still missing. Given this context, the introduction of the phenomological approaches described herein for the determination of constitutive coefficients is desirable due to their ease of use, and resulting explicit relations to underling properties.

Double-porosity mathematical model
We present the balance equations and constitutive laws for the dual-continuum system within the macroscopic framework of Coussy (1995Coussy ( , 2004). The dual system is considered as the superposition of two overlapping poroelastic continua. Elastic deformation of each continuum is thus implied. Quantities denoted by m and f refer to matrix and fracture continua respectively. It is assumed that the poroelastic double-porosity material is isotropic, and is saturated by a slightly compressible fluid which can undergo isothermal flow. Under the assumptions of quasi-static deformations and infinitesimal transformations the momentum balance for the dual medium recovered as where σ is the Cauchy stress tensor, g is the gravity vector, ρ = ρ s (1 − φ) + ρ l φ is the density of the bulk medium, ρ s is the intrinsic density of the solid matrix, ρ l is the intrinsic density of the fluid, and φ is the Lagrangian porosity. This is defined as the ratio of the current pore volume, Ω p , to the bulk volume of the undeformed configuration, Ω 0 , where superscript 0 denotes measurement at reference conditions. Assuming small perturbations in Lagrangian porosity, and solid and fluid densities, allows us to take these quantities at reference conditions where necessary. In keeping with convention, stress is taken as positive in the tensile direction. Finally, γ represents the momentum transfer arising as a result of the mass transfer between the two pore continua. Often γ is assumed to be negligible with respect to the other force density terms (Elsworth and Bai 1992;Pao and Lewis 2002;Fornells et al. 2007;Kim et al. 2012). Next we introduce the linearised strain tensor given according to the straindisplacement compatibility relation where u denotes the displacement vector. The inter-continuum momentum transfer is given by where γ α [α = m, f ] is the rate of mass transfer from pore continuum α to pore continuum β, and v l,α is the absolute fluid velocity within each pore continuum. The rate of mass transfer is conservative between the matrix and the fractures thus The absolute fluid velocity, v l,α , is related to the volume flux, q α , within each continuum by whereṽ l,α is the relative fluid velocity, v s is the velocity of the solid matrix, and φ α is the Lagrangian porosity associated to each continuum. This is defined as where Ω p,α is the current pore volume of continuum i. The volume flux for each pore continuum is then given by Darcy's law where k α and p α denote the permeability tensor and fluid pressures associated with pore continuum, α. The balance of fluid mass for each continuum is then given as where m l,α = ρ l φ α is the fluid mass content of continuum α. We require constitutive laws to provide closure to the model. In early presentations, constitutive relations were indirectly postulated in which inter-continuum pressures were implicitly decoupled (Aifantis 1977(Aifantis , 1979Wilson and Aifantis 1982;Bai 1992, Lewis andGhafouri 1997;Bai et al. 1999). To provide more rigour to the form of the constitutive equations we make use of the energy approach to poromechanics under the assumption of infinitesimal strain theory (Coussy 1995(Coussy , 2004. It can be shown from a purely macroscopic approach (Coussy 2011) or via micromechanical considerations (Dormieux et al. 2006), that the increment in strain work density, dW s , on the skeleton due to the loading triplet (d , dp m , dp f ) can be expressed as Due to elasticity, our system is non-dissipative and thus the skeletal strain energy is stored entirely as an elastic potential where Ψ s denotes the Helmholtz free energy of the skeleton, and from which it follows Eq. (11) is a trivial extension of the skeleton free energy expression for single porosity materials, and is indeed analogous (and identical) to the expression for the multiphase fluid single porosity poromechanical problem (see for example Coussy 2004). Now introducing the following Legendre transform into Eq. (11) results in Next, it is useful to decompose the stress and strain tensors by way of their volumetric and deviatoric parts, where σ = 1 3 tr(σ) is the mean stress, σ d is the deviatoric component of the stress tensor, = tr( ) is the volumetric strain, and d is the deviatoric component of the total strain tensor. Making use of the stress and strain decompositions in Eq. (13), the state equations for double-porosity poroelasticity are then given as Applying Eq. (16) to Eq. (13), and making use of the Maxwell symmetry relations which arise naturally from Eq. (16) whilst also assuming isotropy of the material, we arrive at the constitutive equations for a linear isotropic poroelastic dual-continuum where parameters K dr and G are the drained bulk and shear moduli of the dualmedium respectively. Coefficients b α can be thought of as effective Biot coefficients, and relate changes in fluid volume content to skeletal straining under drained conditions. Coefficients 1 Nα relate changes in the Lagrangian porosity of continuum α to changes in fluid pressure of the same medium, whilst the skeleton remains constrained and fluid pressure in continuum, β, remains constant. Finally, 1 Q is a coupling coefficient that relates changes in fluid content of continuum α and pressure changes in continuum β.
In poromechanics it is common to formulate the constitutive equations in terms of the fluid mass content such that where K l is the fluid compressibility given by where Comparison between Eqs. (18) to (19) and Eqs. (23) to (24) gives the additional relation Under long-term drainage conditions the double-porosity model must be able to reduce to the well known single-porosity model presented (Berryman and Wang 1995;Coussy 2007). This provides us with the following compatibility relations where b is the single-porosity Biot coefficient. A final constitutive equation is required to describe the mass transfer rate between the two pore continua. In accordance with Barenblatt et al. (1960) and Warren and Root (1963) the simplest model for mass transfer between the two continua is given by where k is the interface permeability, which here is taken as the matrix permeability (Barenblatt et al. 1960;Choo and Borja 2015), and η is the shape factor. The first order nature of Eq. (29) may in some cases over-simplify the physics of intercontinuum mass transfer. This should be taken into consideration when working with the dual-continuum paradigm. We make use of an analytical based shape factor as introduced in Lim and Aziz (1995). For an isotropic material in two dimensions this is given as where d denotes the average spacing between the fractures.

Models of constitutive coefficients
We require substantiation of the constitutive coefficients in Eq. (17) and Eqs. (23) to (24). One option is direct measurement of these effective parameters. However, this is predicted to be non-trivial for dual-continua. For example, isolating matrix and fracture contributions would be challenging. An alternative option would be to calculate the effective parameters with models that use more accessible quantities. Examples of more accesible quantities include the intrinsic bulk moduli and effective porosities of each continuum. Such models exist within the literature, and in the following we compare and contrast those from Berryman and Wang (1995), Khalili and Valliappan (1996), Berryman (2002), and Borja and Koliji (2009). We recognise these models, to the best of our ability, as the most dominant within the literature. They have been used in the works of Khalili et al. (2000), Callari and Federico (2000), Pao and Lewis (2002), Fornells et al. (2007), Taron et al. (2009), Kim et al. (2012, Mehrabian and Abousleiman (2014), Choo and Borja (2015) and Choo et al. (2016) to name but a few.

Model approaches and assumptions
3.1.1 Khalili and Valliappan (1996) The authors take a top down approach, postulating macroscopic balance laws and an effective stress expression. Closure for the model equations therein is then sought through thought experiments that isolate volumetric changes of the constituents. Superposition due to linearity, and Betti's reciprocal work theorem finally allow for recovery of the macroscopic behaviour in terms of constituent responses. In doing, expressions for the constitutive coefficients are identified. Khalili and Valliappan (1996) implicitly assume that the fracture phase is all void space. Additionally, the following assumption is also made: b m φ 0 f = b f φ 0 m . This is restrictive, and is later removed in Khalili (2003) and Khalili and Selvadurai (2003), due to the resulting compatibility enforced between bulk moduli (Loret and Rizzi 1999). We present coefficient models derived by the authors without this assumption ( Table 1). The coefficient models from Khalili and Valliappan (1996) (Table 1), are then consistent with results from Loret and Rizzi (1999), Coussy (2005) and Dormieux et al. (2006).

Berryman and Wang (1995)
The model formulation approach taken by Berryman and Wang (1995) shares the same fracture continuum void space assumption as Khalili and Valliappan (1996). Further, the authors also make use of a top down approach. The starting point is the macroscopic constitutive model written within a pure stiffness setting (σ and p α as primary variables; contrary to conventional poromechanical modelling), where the goal is then to find expressions for coefficients A 11 − A 33 . Symmetry of the coefficient matrix in Eq. (31) is implied. A comparison of Eq. (31) with Eq. (17) and Eqs. (23) to (24) reveals the following relations To proceed, the authors formulate expressions for the material coefficients through a series of thought experiments that allow for the isolation of matrix and fracture behaviours. As a result the authors are able to relate Eq. (31) to intrinsic constituent equations. The microscale equations used in Berryman and Wang (1995) are postulated based on the assumption that each constituent belonging to the dual-medium behaves as a Gassmann-material. That is, constituent solid phases are homogeneous and isotropic (Cheng 2016). This assumption is implicit within the approach of Khalili and Valliappan (1996). In the case of a Gassmann-material, the constitutive equations for a material α can be written as where v α denotes the volume fraction of material α, and poroelastic coefficients intrinsic to a material α have been denoted by superscript * . Volume fractions for matrix and fracture materials are given as where Ω 0 m and Ω 0 f are the volumes of matrix and fracture continua at reference conditions respectively. Under the void space assumption v f = φ 0 f . The intrinsic Biot and Skempton coefficients, b * α and B * α for material α respectively, are where K α s is the solid grain modulus of the intact material, and M * α and φ * α are the intrinsic Biot modulus and intrinsic porosity (measurement at reference conditions is implied) respectively.
Khalili and Valliappan  Table 3. Note that our expressions for coefficient models from Khalili and Valliappan (1996) are slightly different to their original presentation. This is done to highlight the similarities between expressions for b f and 1 Expressions for A 11 − A 33 are finally recovered (Table 2), using the prescribed thought experiments. With the relations in Eq. (32) we get material coefficient formulations pertaining to the conventional mixed compliance setting ( and p α as primary variables) ( Table 1). The void space assumption leads to K m s = K s , and thus we recover the same coefficient models as in Khalili and Valliappan (1996).

Borja and Koliji (2009)
In Borja and Koliji (2009), the authors consider the evolution of internal energy density and derive a thermodynamically consistent effective stress law. Aggregate material is used as their reference material. The void space assumption is thus implicit. We present the effective stress expression from Borja and Koliji (2009) for an isotropic single phase dual-continuum system as follows where σ = K dr , and ψ α denotes the pore fraction of continuum α such that Notation ϕ α is the Eulerian porosity of continuum α. This is defined as the ratio of the current pore volume, Ω p,α , to the bulk volume of the current (deformed) configuration, Ω. Notation ϕ s = Ωs Ω is the current volume fraction of the solid phase. In the limit of infinitesimal transformations, ϕ α ≈ φ α , from which, together with small perturbations in Lagrangian porosity, follows ψ α ≈ ψ 0 α . Comparison of Eq. (36) with Eq. (17) leads to the following relations for effective Biot coefficients Borja and Koliji (2009) identify the requirement for a constitutive expression for ψ α based on energy conjugacy with p α (their Eq. (76)). However, explicit constitutive equations for ψ α remain, to the best of the current authors' knowledge, an open question. We do note, however, that Borja and Choo (2016) develop a framework that allows for the tracking of pore fraction evolutions numerically.
As an initial approach to deriving an algebraic expression for variations in ψ α we first consider the following mass balance equation given in Borja and Koliji (2009) We can derive an alternative form of the mass balance here by substitution of Eq. (23) or Eq. (24), together with Darcy's law, into Eq.
where we have used Eq. (26) to decompose 1 Mα . Under the assumption of small perturbations in fluid density, the fourth term on the left-hand side of Eq. (38) can be neglected. Comparing Eq. (38) and Eq. (39) leads to the following identity In the works of Choo and Borja (2015) and Choo et al. (2016), the authors achieve non-algebraic closure of the mass balance equations for each continuum by assuming ∂ψα ∂t ≈ 0. This is equivalent to in Eq. (40). We discard the latter relation in Eq. (41) as it is overly restrictive with respect to fluid exchange between matrix and fractures. Hence we identify the closure condition ∂ψα ∂t ≈ 0 with values for material coefficients 1 Nα and 1 Q of zero. A summary of the coefficient models from Borja and Koliji (2009), under the explicit assumption ∂ψα ∂t ≈ 0, mapped to the constitutive model of Eq. (17) and Eqs. (23) to (24) is shown in Table 1.

Berryman (2002)
The motivation behind the approach by Berryman (2002) is to formulate coefficient models using intrinsic fracture properties. Therefore, contrary to the previous models, no assumption is made on the values of φ * f and K f . This may be useful when assuming that the fracture continuum has an associated stiffness. Homogeneity and isotropy of the constituent solid phases are still assumed.
As per Berryman and Wang (1995), the authors use a top down approach starting with the pure stiffness formulation of the constitutive model (Eq. (31)). To derive expressions for A 11 −A 13 , Berryman (2002) considers scenarios of uniform expansion (or contraction). This is equivalent to asking whether we can find variations in uniform stress, dσ = dσ m = ddσ f , and variations in pore pressures, dp m and dp f , such that = m = f (Berryman 2002). With this thought experiment, models for A 11 − A 13 are recovered accordingly (Table 3). Relations in Eq. (32) allow for inversion of the constitutive coefficients in Table 3 into the standard mixed compliance setting (Table 1).
As a final note on the Berryman (2002) coefficient models, recent users have explicitly assumed A 23 = A 32 = 0 as a closure condition to generalise the dualcontinuum system to a multi-continuum one (Kim et al. 2012;Mehrabian and Abousleiman 2014;Mehrabian 2018), providing explicit relations between material properties and simplifying the coefficient models. It is still unclear how this assumption may affect the system.

In sum
Coefficient models from Khalili and Valliappan (1996), Berryman and Wang (1995) and Borja and Koliji (2009) all make an underlying void space assumption for the high permeability, low storage continuum. We demonstrate expressions for Khalili and Valliappan (1996) and Berryman and Wang (1995) type models are identical under the void space assumption. Models from Borja and Koliji (2009) make use of continuum pore fractions, but still require a final closure relationship for the evolution of each pore fraction. Finally, models from Berryman (2002) make no underlying assumption on the intrinsic porosity, and thus the stiffness of the high permeability, low storage continuum.
In terms of pressure decoupling assumptions we have two types. The first are implicit assumptions for which the constitutive relations are postulated without inter-continuum pressure coupling. The second are explicit assumptions for which the full constitutive model is described (or the requirement for constitutive expressions are at least identified in the case of Borja and Koliji (2009)). These decoupling assumptions are made explicitly so as to provide relations between material properties and to simplify coefficient models (e.g. due to non-algebraic closure).
It is of interest to investigate how differences in coefficient models, in addition to deceoupling assumptions, may impact poromechanical behaviour. This is pursued in the remaining sections.

Micromechanics of dual-continua
To establish the physical implications of explicit decoupling assumptions such as taking A 23 = A 32 = 0 we make use of the theoretical framework of micromechancis. Micromechanics is used as a tool to relate the macroscopic behaviour and properties of a composite to those of its underlying constituents (microscale) (Nemat-Nasser and Hori 2013). In the following we consider the dual-medium as a composite of two (poro-) elastic materials with each material having its own intrinsic constitutive model. Fundamental results from micromechancis can then be used to look at the stress and strain distributions within the composite, as well as relations between composite and constituent moduli.

Volume averaging
The distinction of scales is central to a micromechanical approach. In general, a dual-continuum approach can be justified if a representative elementary volume (REV) can be identified from a large macroscopic structure e.g. a fracture outcrop (Berkowitz 2002). The scale requirements for an REV of size D are that it must be at least an order of magnitude larger than the length scale of local heterogeneities (d), and at least an order of magnitude smaller than the length scale of the macroscopic body (L) (Bear and Bachmat 2012). This is depicted in Fig. 2 and can be summarised as The REV thus represents the scale at which averaged quantities are defined, and at which dual-continuum relationships between these quantities are established (Fig. 2). Typical approaches for averaging include the moving average and the ensemble average (Hashin 1983). The latter is a statistical treatment that represents the Figure 2. Identification of an REV over the microscopic scale from a large macroscopic structure (MS). The REV is used to define the macroscopic dual-continuum (DC) model for which averaged quantities are used. The DC model is the superposition of matrix and fracture continua (or particles) in space and time, whose mass exchange is described by the transfer term γ α . Notations d, D and L denote characteristic lengths of local heterogeneities, the REV and the macroscopic structure respectively. average taken over a composite ensemble whose realisations differ in the spatial distribution of its properties (Hashin 1972). The former can be defined over a volume defined according to Eq. (42) such that where f is an arbitrary field, y is a position vector locating an REV within the composite body, and y is a local vector that locates points within the REV. Volume Ω o is a fixed volume of an REV whose centroid is located at the origin o of a body of size L (Nemat-Nasser and Hori 2013). For various positions of y, Eq. (43) gives a moving average within the composite body.
A fundamental assumption often made in micromechanics is that of statistical homogeneity. This has several implications. First, it implies that the average taken over any REV of a large composite body approaches the average taken over the whole body itself. In this case, the moving average is independent of y, and the shape and size of Ω o , so long as the volume is sufficiently larger than than local heterogeneities. Assuming statistical homogeneity for our dual-material then lets us recover the classical volume averaging expression, where Ω is the volume of any REV within the body, and x = y + y . Secondly, statistical homogeneity allows us to define a link between the moving average and the ensemble average in the form of an ergodic hypothesis. That is, we are able to replace the ensemble average with the body average, or the REV average due to implicit statistical homogeneity, in the limit that the volume of the body goes to infinity (is very large) (Beran 1971;Torquato 2013).

Effective elastic properties
With the averaging operation defined we proceed to formulate the effective elastic property problem. In doing we make use of the works of Hill (1963) and Hashin (1972). Additionally we assume that the composite is drained and thus make no distinction between effective and total stress fields.
Starting at the microscale, for which we have assumed linear elasticity, the micro stress and strain fields, s(x) and e(x) respectively, are related through the drained fourth order stiffness and compliance tensors, C * dr (x) and S * dr (x) respectively, by It can be shown that macroscopic stress and strain tensors are related to their underlying fields through the volume averaging operator (see for example Hashin 1972) such that, where e α = V −1 α Vα e α (x) dV (resp. s α ). Substitution of Eqs. (45) to (46) in Eqs. (47) to (48), and making use of the decompositions in Eqs. (13) to (14) due to the isotropy of the dual-medium leads to where e α and s α are the average microscopic volumetric strain and mean stress fields of continuum α respectively, and S α is the drained compressibility of continuum α.
Fundamental to the micromechanics approach is understanding how macroscopic stress and strain distribute between individual constituents. Due to the linearity of the material behaviour, the following relations hold (Hashin and Shtrikman 1963) where A α and B α are concentration factors mapping macroscopic volumetric strain and mean stress to the equivalent averaged microscopic fields. Additionally, A α (resp. B α ) admit the following compatibility relations Use of Eqs. (51) to (53) in Eq. (49) and Eq. (50) allows us to develop the following fundamental relationships for effective moduli of two phase composites, Finding effective moduli then amounts to determining A α or B α . This is equivalent to determining the ratio of macroscopic stress or strain to the average constituent stress or strain fields respectively. One can do this exactly for dilute concentrations of spherical or ellipsoidal particles of continuum β hosted in continuum α using results from Eshelby (1957). By dilute, we refer to the condition in which particle β behaves as an isolated body embedded in an infinite medium α. There are therefore no pore-stress (mechanical) interactions with other embedded particles of the same continuum (Rice 1997).

Physical implications of explicit decoupling
To relate the micromechanics problem described above with the explicit decoupling assumptions mentioned in Section 3, we consider the form of the constitutive equations given by Eqs. (18) to (19). We show that under certain distributions of strain and stress that decoupling between continuum pressures naturally arises. Further, these distributions lead to explicit relations between composite and constituent moduli, as well as compatibility relations for other constitutive coefficients.
4.2.1 Isostrain: 1 Q = 0 Whilst this explicit assumption has not been used within literature its consideration remains instructive. We define assumption 1 Q = 0 in terms of dφ m (although converse arguments may be used for dφ f ). The term 1 Q = 0 is then equivalent to Next we define the following local constitutive model for the matrix (Dormieux et al. 2006), Under a drained matrix (dp m = 0), Eq. (58) shows that a zero matrix porosity variation (dφ m = 0) can only hold if the average matrix strain is zero (e m = 0). Further, we can see from Eq. (57) that if the matrix is drained and does not deform, then the variation in average matrix stress is also zero (ds m = 0). This will be a useful result for discussions on the explicit decoupling assumption A 23 = 0. Proceeding with Eq. (56) and the macroscopic strain constraint, = 0, the strain partition in Eq. (48) shows that if = e m = 0 then e f = 0. The condition 1 Q = 0 is therefore consistent with a condition of isostrain, This distribution of strain can be obtained when elements are set in parallel to the direction of loading (Fig. 3). Under isostrain, the bulk modulus of the composite is then the Voigt average of the constituent moduli (Voigt 1928), The Voigt average is naturally obtained by substitution of A f = 1 into Eq. (54). Under the void space assumption, isostrain can be summarised as From Table 1, Eq. (62) is implicitly satisfied by the coefficient models from Khalili and Valliappan (1996) when using the Voigt average of the constituent moduli of a void space inclusion composite (i.e. K dr = v m K m ). Further, Table 1 shows compatibilities enforced on effective Biot coefficients will also impact the remaining constitutive coefficients.

Incompressible grain isostrain:
We now consider the coefficient models from Borja and Koliji (2009) under the assumption ∂ψα ∂t ≈ 0 made by Choo and Borja (2015) and Choo et al. (2016). We structure the following in two parts: In Section 3.1.3, we indetified that ∂ψα ∂t ≈ 0 amounts to 1 Nα = 1 Q = 0 when mapping to the constitutive model shown in Eqs. (17) to (19). Our first goal is to see under what conditions the result 1 Nα = 1 Q = 0 arises when using the alternative void space coefficient models. To do so we make use of the coefficient models of Khalili and Valliappan (1996).
Our second goal is to see if there is any correspondence between the effective Biot coefficients, b α , from the coefficient models of Borja and Koliji (2009) and the other void space coefficient models. This is of interest because differences in mass balance behaviour could be introduced through the way in which b α is modelled. To investigate this we use the conditions required to make 1 Nα = 1 Q = 0 with the coefficient models of Khalili and Valliappan (1996) and calculate b α for the various void space coefficient models accordingly.
The set of explicit assumptions: 1 Nα = 1 Q = 0, can easily be derived from the microscale by first considering isostrain ( 1 Q = 0). With the resulting compatibility from isostrain, Eq. (62), along with the assumption K s = ∞, we can see from the coefficient models of Khalili and Valliappan (1996) in Table 1 that 1 Nα = 0 (using Eq. (26) to decompose 1 Mα ). We therefore refer to conditions resulting in 1 Nα = 1 Q = 0 as incompressible grain isostrain, which, in analogy to Eq. (61), can be summarised as As far as parameters in the balance of mass are concerned, Eq. (38) and Eq. (39) are identical when assuming ∂ψα ∂t ≈ 0 in the former and incompressible grain isostrain using the other void space coefficient models in the latter.
With equivalent mass balance expressions we now wish to see if there is correspondence between the effective Biot coefficients.
From Eq. (61) or Eq. (62) the compatibilities for the effective Biot coefficients under incompressible grain isostrain read Table 1, the coefficient models of Borja and Koliji (2009) are unable to satisfy these compatibilities. Instead the effective Biot coefficients given by this set of models are equivalent to pore fractions, since b = 1 for incompressible grains. In contrast, the coefficient models of Khalili and Valliappan (1996) do satisfy the compatibility requirements. Due to the differences in effective Biot coefficients (and thus other constitutive parameters), we expect disparity in poroelastic behaviour when using the two sets of void space coefficient models.

Isostress:
Finally, we study the pure stiffness setting with the condition A 23 = 0 in light of the assumptions made in Nguyen and Abousleiman (2010), Kim et al. (2012), Mehrabian and Abousleiman (2014) and Mehrabian (2018). We continue to work in terms of dφ m . Accordingly In Section 4.2.1 we established that a drained matrix experiencing no deformation is concurrent with a condition of zero (average) matrix stress. From the macroscopic stress constraint in Eq. (64), dσ = 0, and the stress partition in Eq. (47), if dσ = ds m = 0 then ds f = 0. The condition A 23 is therefore consistent with a condition of isostress, The classical configuration under which an isostress distribution is observed, are elements that are arranged transversely to the direction of applied load (Fig. 4a). For isotropic fracture networks, a more frequent configuration that shows isostress behaviour is within a solid-fluid suspension (Fig. 4b). Such a situation may occur if a network of open fractures totally permeates the solid, thus completely dissociating the matrix material. Under isostress, the bulk modulus of the composite is then the Reuss average of the constituent bulk moduli (Reuss 1929), The Reuss average is naturally obtained by substitution of B f = 1 into Eq. (55), or from A 23 = 0 in Table 3.  When K f is zero, as is the case for a void space fracture phase, we can see from Eq. (66) that K dr is also zero. Using this result and the isostress condition with Eq. (17) and Eq. (57) for the matrix and fracture phases gives us the following b m dp m − b f dp f = K m e m − b * m dp m = −1dp f .
From Eq. (67) we can establish the following compatibility relations for the effective Biot coefficients assuming a void space fracture phase From Table 1, Eq. (68) is implicitly satisfied by the coefficient models from Khalili and Valliappan (1996) when using the Reuss average of the constituent moduli of a void space inclusion composite (i.e. K dr = 0).  (2018), isostress is implicitly assumed. Upscaling of constituent moduli is then admitted through the Reuss average. This raises the question as to whether this is a reasonable approach to upscaling or not? The Reuss and Voigt average represent lower and upper bounds for effective moduli respectively (Hill 1952). It is instructive to therefore first discuss the upscaling problem, and the Reuss average, within the context of bounds.

On moduli upscaling
Moduli bounds represent a minimum and maximum limit, between, or at which, effective moduli arise. For the bulk modulus, bounds by Hashin and Shtrikman (1963) have been shown in the same paper to be the best possible for isotropic composites, given only constituent moduli and volume fractions. This result is obtained by solving exactly for the bulk modulus of a specific geometry composite (composite sphere assemblage), and observing that the resultant solution coincides with either bound depending on the stiffness of the inclusion relative to the host material. The lower Hashin-Shtrikman (HS) effective bulk modulus bound is given as where G f denotes the shear modulus of the fracture material. The upper bound, K HS+ dr can be determined by swapping subscripts f and m in Eq. (69). When G f = 0, such as in a fluid suspension geometry, the HS lower bound and the Reuss bound coincide. It follows that in this geometry, with the stiffness of a void space fracture phase (K f = 0), both the Reuss and HS lower bound result in K dr = 0. However, in situations when the fracture phase has an intrinsic stiffness, the Reuss and HS lower bounds may be significantly different (Watt et al. 1976).
Bounds can be used as a first approach to upscaling under certain geometries. If a fracture network completely percolates a matrix then it will have the maximum effect of weakening the rock (Watt et al. 1976). The effective bulk modulus of the composite will then coincide with the HS lower bound (Boucher 1974;Watt et al. 1976). In the general case fractures are likely to have an associated stiffness (Bandis et al. 1983). We thus recommend using the HS lower bound over the Reuss average as a first approach to upscaling for such geometries. This procedure is also in line with the assumptions built into the continuum approach: The continuum assumption is linked to ones ability to define an REV over which properties can be averaged. For a fractured system, such an REV cannot be justified if the system is poorly connected (Berkowitz 2002). Use of the HS lower bound, as a first approach to moduli averaging, thus supports the notion of a well-connected isotropic dense fracture network, over which an REV could be defined.
When the underlying composite geometry precludes the use of bounds as methods for upscaling, one must use other methods of averaging. In addition, upscaling procedures need to be selected on whether to incorporate interaction effects of the high permeability, low storage porosity with stress. Under the dilute assumption such interactions are negligible (Rice 1997). In this case direct approaches can be used to calculate A α (or B α ) (see for example Hashin (1983) and Dormieux et al. 2006). At higher values of φ f , mechanical pore-stress interactions need to be accounted for in the averaging procedure. Such procedures include, but are not limited to, the self-consistent method (e.g. Hill 1965) and the differential scheme (e.g. Norris 1985). Upscaling with these schemes generally requires iterative procedures due to coupling between upscaled bulk and shear moduli. Additionally, we note that for fractured media, it has been shown that the fracture density is the controlling factor on effective elastic moduli (Budiansky and O'connell 1976). This is because fractures can interact even at vanishingly small fractions of fracture material (Rice 1997). Thus schemes that account for mechanical interactions are also recommended for more general fracture arrangements (e.g. when the fracture phase does not completely dissociate the rock).

Qualitative analysis using the Mandel problem
We use solutions to the double-porosity Mandel-problem to investigate the physical impacts of different coefficient modelling concepts and assumptions on the poromechanical response of a dual-continuum material. We consider implicit assumptions, where the constitutive model starts with no pressure coupling, and explicit assumptions, where the full constitutive model, Eq. (17) and Eqs. (23) to (24), is the starting point but pressure coupling is neglected in order to introduce explicit relations between bulk and constituent moduli, and to simplify coefficient models.
We first investigate the effects of considering the fractured dual-medium as a void space inclusion composite or stiff inclusion composite as assumed by the coefficient models of Khalili and Valliappan (1996) and Berryman (2002) respectively. Second, we study the effects of decoupling assumptions.

Double-porosity Mandel problem
The problem geometry is described as an infinitely long (rectangular) cuboid domain such that the plane-strain condition holds (i.e. u y = 0 and q m,y = q f,y = 0) (Fig. 5). The domain is sandwiched between two impermeable, rigid plates, and is free to displace both laterally and vertically. A constant compressive force, a −a σ zz dx = −2F a, is applied at the rigid plate boundaries, Γ N and Γ S (north and south boundaries respectively). The east and west boundaries, Γ E and Γ W respectively, are then free to drain such that p m = p f = 0 at these boundaries.
In summary the boundary conditions are, Due to the symmetry of the problem only a quarter of the domain need be considered.
where we assume γ ≈ 0. It can be shown that ∇ 2 u = ∇ . With this and in the absence of body forces, Eq. (72) reduces to where c α = bα(1−2v) 2G(1−v) is the consolidation coefficient belonging to material α.

Data for analysis
For the qualitative analysis we use a quarter of a 2 m × 2 m deformable porous domain. The studied domain is subjected to a constant top boundary force, 1 0 −2× 10 6 Pa dx = −2 MPa.m. Where possible we use values for material properties that are typically encountered in naturally fractured carbonate reservoirs. The mechanical properties are then K m = 20 GPa, K s = 70 GPa, and v = 0.2 (Wang 2000). Values for K dr and K f are problem dependent. In cases where no rigorous justification is used, we choose values for K dr and K f arbitrarily (denoted by a superscript †). Fluid properties are for that of water; ρ 0 l = 1000 kg m -3 , µ l = 1 cp and K l = 2.3 GPa. Rock properties related to fluid storage and flow are φ 0 m = 0.05, k m = 0.01 md and φ 0 f = 0.01 (Nelson 2001). Effective permeability of the fracture is assumed to be k f ≈ 1000 md. With the cubic law (Witherspoon et al. 1980), this corresponds to an aperture, a f = 7 × 10 −5 m, and fracture spacing, d = 2.8 × 10 −2 m.

Test cases
We consider one test case to investigate the differences between void space and stiff inclusion composite coefficient models, and three test cases to investigate implicit and explicit decoupling assumptions.

Case 1: Intrinsic fracture properties
In Case 1 we are interested in comparing the differences that arise when considering intrinsic fracture properties. In particular, it is of interest to investigate if coefficient models from Khalili and Valliappan (1996) and Berryman and Wang (1995) could still be used even when a fracture has an associated phase stiffness. We hypothesise that provided φ * f ≈ 1 and the fracture phase stiffness is orders of magnitude lower than the stiffness of the solid material (K f K f s ), effects arising due to deviations of intrinisc fracture coefficients from the void space assumption (b * f = 1, 1 and B * α = 1) are negligible. Void space coefficient models would then be usable in place of composite coefficient models We consider a composite medium with a network of fractures that completely dissociates the matrix. Upscaling is then done through the HS lower bound. To test our hypothesis we use combinations of various values of φ * f and K f . When φ * f = 1 we use coefficient models from Khalili and Valliappan (1996). When φ * f < 1 we use coefficient models from Berryman (2002) with K f s = K s . We compare results from coefficient models calculated using a fracture modulus several orders of magnitude lower than the solid modulus (K † f = Km 500 ), versus ones that use a fracture modulus only an order of magnitude lower than the solid modulus (K † f = Km 10 ).

Case 2: Implicit decoupling assumptions
For Case 2 we investigate the impact of implicit decoupling assumptions. This is particularly poignant considering the use of such decoupled constitutive models in the recent works of Alberto et al. (2019) and Hajiabadi and Khoei (2019). To mimic implicit assumptions we consider 1 Q = 0 and A 23 = 0 and make no acknowledgement of these assumptions with respect to relations between mechanical properties. When considering 1 Q = 0 and A 23 = 0 we use coefficient models from Khalili and Valliappan (1996) and Berryman and Wang (1995) respectively. We reference these results against ones for which no decoupling is made with coefficient models from Khalili and Valliappan (1996).
Use of the void space models implies K f = 0. Since we do not enforce relations on K dr (as per requirement of implicit assumptions) we take an arbitrary value of K † dr = 10 GPa. As a result compatibilities for b m and b f that arise when making explicit assumptions do not hold.

Case 3: Explicit decoupling assumption -isostrain
We investigate the effect of assuming isostrain at the microscale whilst making use of coefficient models from Khalili and Valliappan (1996). Under isostrain the composite and constituent bulk moduli are linked by the Voigt average (K dr = 19.8 GPa). This leads naturally to 1 Q = 0, as well as additional compatibility requirements on the effective Biot coefficients which go on to affect the remaining constitutive coefficients. We compare the isostrain results to those computed when using coefficient models with a composite bulk modulus coinciding with the HS upper bound and an arbitrary value (K dr = 19.5 GPa and K † dr = 10 GPa respectively). Further, we investigate the disparity between results when using the void space coefficient models of Borja and Koliji (2009) under the assumption of ∂ψα ∂t = 0, and Khalili and Valliappan (1996) under the assumption of incompressible grain isostrain. We use K dr = 19.8 GPa and K s = ∞ for both sets of coefficient models in this latter isostrain investigation.

Case 4: Explicit decoupling assumption -isostress
For Case 4 we study the effect of assuming isostress at the microscale. In previous works the coefficient models of Berryman (2002) have been used with an explicit decoupling assumption that implies isostress (Kim et al. 2012;Mehrabian and Abousleiman 2014;Mehrabian 2018).
To avoid cases where K dr = 0 we consider the fracture phase to have the following properties: φ * f = 0.7 with K f s = K s and K † f = Km 500 . Coefficient models from Berryman (2002) are then used.
We compare results arising when calculating the composite bulk modulus with the Reuss average (K dr = 3.3 GPa), the HS lower bound (K dr = 5.7 GPa), and an arithmetic average of the HS bounds (K dr = 12.7 GPa). The latter modulus is tested in analogy to a dual system with inclusions that do not have the maximum weakening effect on the host material. One example would be a fracture system composed of a network of open and closed fractures. Another would be aggregate material.

Results and discussions
In the following, we show results of the test cases described above for the doubleporosity Mandel problem. Results are given in terms in evolutions of matrix and fracture pressures, and vertical strain with time.
To aid in our analysis for pressure and vertical strain we introduce the following notions of the instantaneous problem and the time-dependent problem. In both cases, mechanical equilibrium is governed, for this system, by Eq. (72). In the instantaneous problem, fluid pressure for continuum α can be shown to be a state function of total stress and fluid pressure β. In the time-dependent problem, continuum pressures are governed by way of the diffusion equation, Eq. (39).
The instantaneous problem considers the change of the system from an unloaded state, t(0), to a loaded state, t(0 + ), upon application of instantaneous loading. Under such conditions the rate of loading is infinitely faster than the rate of inter-continuum fluid transfer (Coussy 2004). Consequently, each continuum is undrained, m l,α . From Eq. (23) and Eq. (24), together with the undrained condition (dm l,α = 0) we can recover We note that zz ∝ . Substitution of Eq. (17) into Eq. (75) leads to 6.1 Case 1: Intrinsic fracture properties Fig. 6a and Fig. 6c show matrix and fracture pressure evolutions whilst varying intrinsic fracture properties. We find a good match in both matrix and fracture pressure evolutions using coefficient models from Khalili and Valliappan (1996) and Berryman (2002) when the fracture is almost all void space (φ * f ≈ 1) and the fracture phase stiffness is orders of magnitude lower than the stiffness of the solid material (K f K f s ) (Fig. 6a). Even when the intrinsic fracture porosity is diminishing (φ * f = 0.2), provided the fracture phase stiffness is orders of magnitude lower than the solid stiffness, the difference between early time matrix pressures with the different coefficient model formulations is small (Fig. 6a). However, when the fracture bulk modulus is only an order of magnitude lower than the solid modulus (K f K f s ), early time fracture pressure differences become measurable as intrinsic porosity decreases (Fig. 6c). This phenomenon can be explained by considering the intrinsic Skempton coefficient, B * α , written in the following form for the fracture phase (Cheng 2016),      ).
Use of Eq. (77) allows us to cast our observations as a bounding problem such that The void space assumption is concurrent with the lower bound of B * f in Eq. (78). Additionally, Eq. (77) show that as φ * f approaches zero, B * f must asymptotically approach one.
If the lower bound in Eq. (78) for a given K f is close to one then changes in intrinsic fracture porosity are negligible due to the proximity of the lower and upper bounds. This is the case when K f K f s . If the fracture phase stiffness is not orders of magnitude lower than the solid stiffness, then the lower bound of B * f may be significantly less than one. In this case, changes in B * f cannot be captured when using void space coefficient models, and thus early time fracture pressure is underestimated as φ * f decreases (Fig. 6c). Fig. 6b and Fig. 6d show the variation in vertical strain for the softer and stiffer fracture phases respectively. In Fig. 6b, the strain evolutions are almost identical when using the more compliant fracture phase for the range of intrinsic fracture porosities. This is coupled to the similarity in pressure evolutions. Whilst induced fracture pressures are significantly larger when fracture stiffness and intrinsic porosity are non-negligible, Fig. 6d shows very little differences in vertical strain across the whole intrinsic porosity range. This is a direct consequence of the algebraic coefficient of strain in Eq. (75) for the fracture phase, which scales proportionally with the variation in fracture pressure.

Case 2: Implicit decoupling
Pressure and vertical strain results for the implicit decoupling assumptions test case are presented in Fig. 7. When assuming A 23 = 0, the matrix and fracture pressure evolutions are almost identical to the reference case (Fig. 7a), with vertical strains also being correspondingly very similar (Fig. 7b). When assuming 1 Q = 0, the early time matrix and fracture pressures are measurably lower than the reference case (Fig. 7a). The early time vertical strain is greater when 1 Q = 0 than the strain in the reference case (Fig. 7b).
The results in Fig. 7a suggest that the assumption 1 Q = 0 has the most noticeable effect on the poromechanical behaviour of the dual-medium. As fracture pressure increases we expect contraction of the matrix in order to accommodate fracture expansion under a constrained strain condition (vice-versa fracture contraction under matrix expansion). This requires 1 Q < 0 (Berryman and Wang 1995). From Eq. (76), setting 1 Q = 0 thus has the effect of removing a pressure source from continuum α. Removal of this poromechanical coupling explains the lower than expected induced matrix and fracture pressures. From Eq. (75) we  Khalili and Valliappan (1996) for which no assumptions have been made can see dp α ∝ −1 . Underestimated pressures therefore explain the over estimated strain when taking 1 Q = 0. In contrast, Eq. (32) shows that although A 23 = 0, 1 Q = 0. The pressures in each continuum therefore remain poromechanically coupled with respect to the mixed compliance setting, hence the similarities in pressures and vertical strain when A 23 = 0 versus the reference case.
The current results suggest that assuming A 23 = 0 is a reasonable implicit assumption to make. However, we advise caution when interpreting this result. Based on the results in this section it would be easy, and incorrect, to use them as a justification for assumptions made at the microscale. In Section 4.2 it was shown that explicit assumptions affect all of the constitutive coefficients due to compatibilities enforced between bulk moduli. The remainder of this results section aims at qualitatively supporting these findings. Fig. 8a gives the results for the comparison between coefficient models calculated with different values of K dr , pertaining to cases of isostrain (Voigt average), the HS upper bound and an arbitrary value. In the isostrain and HS upper bound cases, matrix pressure evolutions appear to be almost identical. However, for the arbitrary K dr case the fracture pressures are significantly higher. Further, for the arbitrary case the matrix pressures are slightly higher at early to middle times relative to the aforementioned upper bound cases. The contrast in fracture pressures is most pronounced between the isostrain and arbitrary composite bulk modulus cases. For the induced problem, different values of K dr , which go on to affect constitutive coefficient calculations, explain the differences in pressure by way of Eq. (76).

Case 3: Explicit decoupling -isostrain
An alternative, heuristic approach to explaining the pressure distributions in Fig. 8a can be achieved by considering the required geometry that would be necessary to give a K dr that is lower than the Voigt average. Departure from this upper bound occurs when inclusions are arranged so as to weaken the composite. They then take on a greater portion of the distributed stress. We therefore observe a greater induced fracture pressure when using a lower value of K dr compared to the Voigt bound, due to the proportionality between stress and pressure.  evolutions for the double-porosity Mandel problem considering isostrain. 'ISN' denotes coefficient models for which isostrain is assumed (K dr = 19.8 GPa). 'HS + ' are models calculated using the upper Hashin-Shtrikman bound for K dr (19.5 GPa). 'K † dr ' are coefficient models calculated with an arbitrary bulk modulus of 10 GPa.
The early time vertical strains in Fig. 8b can be explained by way of Eq. (75) or through the heuristic argument. With the latter, towards the upper bounds of K dr the matrix supports the majority of deformation. Since the matrix is stiff, deformation is low. In contrast, when fractures are arranged such that they have a greater weakening effect on the solid, deformation is high. At later times vertical strain between both the upper bound and arbitrary K dr cases diverge. Towards the upper bounds for K dr , b f < b m . In the case of isostrain this fact is easily seen from the relations in Eq. (62). The magnitude of b m means that deformation is more strongly coupled to differences in matrix pressure relative to fracture pressure by way of momentum balance, Eq. (72). This explains the growth in vertical strain separation at later times shown in Fig. 8b.
Of further interest, considering these both represent upper bounds on the composite bulk modulus, is the difference in early time fracture pressures between the isostrain and HS upper bound cases. The early time fracture pressure associated with the HS upper bound is over double that of the isostrain case. This highlights the need for caution before making assumptions on the distribution of strain between constituents. Fig. 9a displays the results for the incompressible grain isostrain investigation. For the case of coefficient models from Khalili and Valliappan (1996), the induced matrix pressure is significantly larger than the induced fracture pressure. Conversely, when using coefficient models from Borja and Koliji (2009) with the assumption ∂ψα ∂t ≈ 0, induced matrix and fracture pressures are equal. This can be explained by considering Eq. (75) which allows us to equate induced variations in matrix and fracture pressures such that Assuming ∂ψα ∂t ≈ 0 together with the effective Biot coefficient expressions from Borja and Koliji (2009) (Table 1), Eq. (79) reduces to from which it is easy to see that dp m = dp f . Under isostrain, we would expect the distribution of stress required to maintain strain uniformity between matrix and fractures to lead to disparate matrix and fracture pressures. The result dp m = dp f therefore suggests that the closure condition ∂ψα ∂t ≈ 0 may be an even stronger assumption than incompressible grain isostrain alone. Fig. 9b shows vertical strain is lower at early times when using coefficient models from Khalili and Valliappan (1996) under incompressible grain isostrain. This can be explained by Eq. (75), which is affected by differences in b α arising from each set of coefficient models.
In light of the discussions in Section 4.2.2 and the results presented herein, assuming ∂ψα ∂t ≈ 0 appears to be a strong closure assumption to make. Thus, we suggest further development of a constitutive model for ψ α , along with its relationship to the constitutive model shown in Eqs. (17) to (19).   Khalili and Valliappan (1996) and Borja and Koliji (2009) respectively. The latter includes the assumption ∂ψα ∂t ≈ 0. This allows us to map the effective stress model from Borja and Koliji (2009) onto the constitutive model Eqs. (17) to (19). Figure 10a shows the impact of the upscaling method on matrix and fracture pressure evolutions. A stiffer composite bulk modulus leads to a lower induced fracture pressure and earlier onset of pressure diffusion in the same continuum. This is the case when using the arithmetic mean of the HS bounds. Non-monotonic rises in matrix and fracture pressures are more pronounced for more compliant composites. This is the case when assuming isotress (Reuss average) and the HS lower bound. Additionally, more compliant composites exhibit a faster decrease in matrix pressure at later times when compared to the stiffer modulus case.

Case 4: Explicit decoupling -isostress
Non-monotonic pressure rises are often referred to as the Mandel-Cryer effect within literature (Wang 2000;Cheng 2016). We predict that such rises lead to greater pressure differentials at middle to late times between the matrix and fracture domains. This would explain the higher rate of matrix pressure diffusion in the compliant material cases. Fig. 10b shows pronounced distinctions in vertical strain for the three different upscaling cases. As expected, the stiffer arithmetic mean of the HS bounds shows the lowest deformation. Of more interest, considering that they both correspond to lower bounds, are the difference in strains between the cases of isostress and HS lower bound. Vertical strain at late times is approximately 75% larger when assuming isostress compared to when using the HS lower bound.  Notation 'HS − ' and 'AHS' are models calculated using the lower Hashin-Shtrikman bound and the arimthmetic mean of the Hashin-Shtrikman bounds for K dr respectively (5.7 GPa and 12.7 GPa respectively).
The cause of the difference in vertical strains between the isostress and HS lower bound cases can be explained using similar discussions as to those used in Sections 6.2 and 6.3. First, using the heuristic argument the HS lower bound is higher than the Reuss bound suggesting that the matrix is capable of supporting a greater distribution of strain. This explains the difference in vertical strain at early times. Late time differences can be explained by the differences in fracture pressure. In contrast to Section 6.3, towards the lower bounds for K dr , b f > b m . The magnitude of b f means that deformation is more strongly coupled to differences in fracture pressure relative to matrix pressure by way of momentum balance, Eq. (72).
With a view toward multi-continuum generalisations, based on the results in Section 4.2.3 and the qualitative results herein we recommend care before assuming isostress. This stress distribution has strong geometrical implications that, without experimental substantiation to prove otherwise, would seem unlikely to hold within a multi-continuum material.

Conclusion
The goals of this paper were twofold. The first was to compare and contrast several prevelant poromechanical dual-continuum constitutive coefficient models. We considered coefficient models from Khalili and Valliappan (1996), Berryman and Wang (1995), Koliji (2009) andBerryman (2002). The first three sets of models constitute a fracture continuum that is all void space. The latter allows for the addition of fracture stiffness effects.
The second goal was to investigate inter-continuum pressure decoupling assumptions made at various stages of modelling. At the first stage we investigate implicit assumptions in which constitutive equations are postulated without pore domain pressure coupling. To follow we investigated explicit assumptions where decoupling is explicitly done on (complete) constitutive equations which include pressure coupling.
Based on the findings herein we make the following conclusions.
On constitutive coefficient models: 1. Coefficient models from Khalili and Valliappan (1996) and Berryman and Wang (1995) are equivalent due to the void space assumption. Algebraic closure for models from Borja and Koliji (2009) is still an open question. Without such closure, poromechanical behaviour is measurably different versus using other void space coefficient models.
2. Measurable differences in poromechanical behaviour are observed between void space and stiff fracture coefficient models when the fracture phase has non-negligible intrinsic porosity, and stiffness with respect to the solid material (φ * f ≈ 1 and K f K f s respectively). On decoupling assumptions: 3. Mechanical decoupling of inter-continuum pressures can lead to misleading poromechanical behaviours. An implicit assumption which neglects intercontinuum pressure coupling within a mixed compliance setting is equivalent to removing a fracture and matrix pressure source for both matrix and fracture phases respectively. Neglecting inter-continuum pressure coupling within a pure stiffness setting still allows for inter-continuum pressure coupling. However, the latter result should not be used as a justification for explicit assumptions. These latter assumptions impact all of the constitutive coefficients, not just those coupling continuum pressures.
4. In a complete mixed compliance constitutive model, explicitly decoupling between continuum pressures is equivalent to enforcing equality between matrix and fracture strains (isostrain). Effective moduli are then linked to constituent moduli through the Voigt average. This leads to compatibility requirements between remaining constitutive coefficients. The Voigt average is an upper bound on composite moduli.
5. In a complete pure stiffness constitutive model, explicitly decoupling between continuum pressures is equivalent to enforcing equality between matrix and fracture stresses (isostress). Such a state may be observed in systems of solidfluid suspension. Under the isostress assumption, effective moduli are linked to constituent moduli through the Reuss average. This leads to compatibility requirements between remaining constitutive coefficients. The Reuss average is a lower bound on composite moduli.
6. If the effective moduli are not available then averaging of constituent moduli is required. This should be done with the geometry and type of constituents of the composite in mind. In special cases, bounds can be used in place of effective moduli. In these instances we advocate the use of the optimal Hashin-Shtrikman bounds in place of the Reuss-Voigt bounds. Use of the latter may result in measurable under or over estimations of effective moduli. In more general cases, averaging procedures, in addition to the aforementioned considerations, should also account for pore interaction effects when necessary.
Recommendation: For simplicity we recommend the use of either coefficient models from Khalili and Valliappan (1996) or Berryman and Wang (1995). When effective moduli are not available we still recommend the use of these models, provided φ * f ≈ 1 and K f << K f s , as all that is required is averaging of constituent moduli. If the above conditions cannot be satisfied then we recommend the use of coefficient models from Berryman (2002), along with the appropriate averaging procedures. Finally, we advise caution before assuming isostrain and isostress as these represent end-member states that are unlikely to be observed in isotropic multiscale media.