Foundations and their practical implications for the constitutive coefficients of poromechanical dual-continuum models
Mark Ashworth Florian Doster

Self-Supervised Learning via Conditional Motion Propagation
    [
    December 30, 2023
===========================================================

[figure]labelfont=bf,labelformat=default,labelsep=period,name=Figure




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, a holistic overview of the different modelling concepts 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. Within the framework of micromechanics we show heuristically that these explicit decouplings are in fact coincident with bounds on the effective parameters themselves. Depending on the formulation, these bounds correspond to end-member states of isostress or isostrain. We show the impacts of using constitutive coefficient models, decoupling assumptions and parameter bounds on poromehcanical behaviours using analytical solutions for a 2D model problem. Based on the findings herein we offer recommendations for how and when to use different coefficient modelling concepts.





§ 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 (; ; ; ; ; ). 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) ().  
   
With respect to fractured systems, using explicit methods can be computationally prohibitive at large scales (; ; ). Additionally, explicit methods may require data (e.g. spatial data) that is not obtainable without direct access (; ). 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 <cit.>. 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 phenomenon is known as poromechanical coupling (), and is described by the well established poromechanical theory (see for example ; ; ; ; ; ).  

<cit.> and <cit.> were the first to introduce the generalised notion of deformation within the dual-continuum setting. Further offerings then came from <cit.>, <cit.> and <cit.>. 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 <cit.>. 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 <cit.> still remains an open question. 

Additional models in which the constitutive equations included inter-continuum pressure coupling were introduced by <cit.>, <cit.>, <cit.>, <cit.>, <cit.>,  <cit.> and <cit.>. 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. ), whilst others allow for an intrinsic phase stiffness (e.g. ). 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 ; ; ; ). 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 aim of this paper is to formulate a set of recommendations for how and when to use different constitutive modelling concepts. In doing we show the impacts of making implicit and explicit decoupling assumptions. In the case of the latter we use heuristic arguments from micromechanics to show that these assumptions are coincident with bounds on the effective parameters themselves.  

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 (). 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 and when to use coefficient models in light of (a) intrinsic fracture stiffness effects and (b) pressure decoupling assumptions between pore domains. 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 (<ref>).     

We note that work has been done on the determination of effective properties of multiple-porosity materials via homogenisation methods (e.g.  and ). 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 (). 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 (), 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 <cit.>. 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


    ∇·σ + ρg = γ,


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 property 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 (; ; ; ). 

Next we introduce the linearised strain tensor given according to the strain-displacement compatibility relation


    ϵ = ∇^symu = 1/2(∇u + ∇^⊤u),


where u denotes the displacement vector.

The inter-continuum momentum transfer is given by


    γ = ∑_α = m,fγ_αv_l,α,


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 


    ∑_α = m,fγ_α = 0.


The absolute fluid velocity, v_l,α, is related to the volume flux, q_α, within each continuum by


    q_α = ϕ_αṽ_̃l̃,̃α̃ = ϕ_α(v_l,α - v_s),


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


    ϕ_α= Ω_p,α/Ω^0,  such that ϕ = ∑_α = m,fϕ_α,
   

where Ω_p,α is the current pore volume of continuum α. The volume flux for each pore continuum is then given by Darcy's law


    q_α = -k_α/μ_l(∇ p_α-ρ^0_lg),


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 


    m_l,αt + ρ^0_l∇·q_α = γ_α,


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 (; ; , ; ). 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 (). It can be shown from a purely macroscopic approach () or via micromechanical considerations (), 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


    dW_s = σdϵ + p_mdϕ_m + p_fdϕ_f.


Due to elasticity, our system is non-dissipative and thus the skeletal strain energy is stored entirely as an elastic potential


    dW_s = dΨ_s,


where Ψ_s denotes the Helmholtz free energy of the skeleton, and from which it follows 


    σdϵ + p_mdϕ_m + p_fdϕ_f - dΨ_s = 0.


<ref> 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 ). Now introducing the following Legendre transform 


    F_s = Ψ_s - p_mϕ_m - p_fϕ_f,


into <ref> results in


    σdϵ - ϕ_m dp_f - ϕ_f dp_f - dF_s = 0.


Next, it is useful to decompose the stress and strain tensors by way of their volumetric and deviatoric parts, 


    σ = σ1 + σ_d, 
    ϵ = 1/3ϵ1 + ϵ_d,


where σ = 1/3tr(σ) 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 from eqn:14eqn:15 in <ref>, the state equations for double-porosity poroelasticity are then given as


    σ = ∂ F_s/∂ϵ;  σ_d = ∂ F_s/∂ϵ_d;  ϕ_m = ∂ F_s/∂ p_m;  ϕ_f = ∂ F_s/∂ p_f.


Applying <ref> to <ref>, and making use of the Maxwell symmetry relations which arise naturally from <ref> whilst also assuming isotropy of the material, we arrive at the constitutive equations for a linear isotropic poroelastic dual-continuum 


    dσ = K_drϵ - b_mdp_m - b_fdp_f, 
    dϕ_m = b_mϵ + 1/N_mdp_m + 1/Qdp_f, 
    dϕ_f = b_fϵ + 1/Qdp_m + 1/N_fdp_f. 
    dσ_d = 2Gϵ_d


where parameters K_dr and G are the drained bulk and shear moduli of the dual-medium respectively (). Coefficients b_α can be thought of as effective Biot coefficients, and relate changes in effective Lagrangian porosity 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 the Lagrangian porosity 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


    dm_l,α/ρ_l = dϕ_α + ϕ_αdρ_l/ρ_l≡dϕ_α + ϕ_αdp_α/K_l,
 

where K_l is the fluid compressibility given by


    1/K_l = 1/ρ_ldρ_l/dp.
  

With <ref> we can express eqn:18eqn:19 as


    dξ_m =  b_mϵ + 1/M_mdp_m + 1/Qdp_f, 
    dξ_f = b_fϵ + 1/Qdp_m + 1/M_fdp_f,


where 


    dξ_α = dm_l,α/ρ^0_l,


Comparison between eqn:18eqn:19 and eqn:23eqn:24 gives the additional relation 


    1/M_α = 1/N_α + ϕ^0_α/K_l.
  

Under long-term drainage conditions the double-porosity model must be able to reduce to the well known single-porosity model. This provides us with the following compatibility relations


    b = b_m + b_f = 1 - K_dr/K_s, 
    1/M_α + 1/Q = b_α-ϕ^0_α/K_s + ϕ^0_α/K_l.


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 <cit.> and <cit.> the simplest model for mass transfer between the two continua is given by


    γ_m = ηρ^0_lk/μ_l(p_f - p_m),     γ_f = ηρ^0_lk/μ_l(p_m - p_f),


where k is the interface permeability, which here is taken as the matrix permeability (; ), and η is the shape factor. The first order nature of <ref> may in some cases over-simplify the physics of inter-continuum 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 <cit.>. For an isotropic material in two dimensions η is defined as     


    η = 2π^2/d^2,


where d denotes the average spacing between the fractures. 



§ MODELS OF CONSTITUTIVE COEFFICIENTS

We require substantiation of the constitutive coefficients in <ref> and eqn:23eqn:24. One option is direct measurement of these effective parameters. However, this approach 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. In the following we compare three modelling concepts that make use of different properties for the calculation of the constitutive parameters: 



  * <cit.> - Constituent mechanical properties, assuming the high permeability, low storage continuum is all void space (no intrinsic fracture properties),

  * <cit.> - Constituent pore fractions, assuming the high permeability, low storage continuum is all void space,

  * <cit.> - Constituent mechanical properties, including intrinsic fracture properties.


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 <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.> and <cit.>. It should be stressed that in most cases the modelling concepts introduced in this section build on, or are aligned with, previous works and concepts introduced by <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.> and <cit.> to name but a few. These works should thus by borne in mind in the section to follow. 



 §.§ Model approaches and assumptions



  §.§.§ <cit.>

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.  

<cit.> 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 relation is restrictive, and is later removed in <cit.> and <cit.>, due to the resulting compatibility enforced between bulk moduli (). We present coefficient models derived by the authors without this assumption (<ref>). The coefficient models from <cit.> (<ref>), are then consistent with results from <cit.>, <cit.>, and <cit.>. 


   


  §.§.§ <cit.>
 
In <cit.>, 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 <cit.> for an isotropic single phase dual-continuum system as follows 


    σ' = dσ + ∑_α=m,fψ_α bdp_α,
  

where σ'=K_drϵ, and ψ_α denotes the pore fraction of continuum α such that


    ψ_α=φ_α/1-φ_s.
  

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 <ref> with <ref> leads to the following relations for effective Biot coefficients, b_m = ψ^0_m b and b_f = ψ^0_f b. 

<cit.> 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 <cit.> 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 <cit.>,


    b_αϵt + (1-ϕ_s)ψ_αt + ϕ^0_α/K_lp_αt+q_α/K_l∇ p_α + ∇·q_α  = 1/ρ^0_lγ_α.


We can derive an alternative form of the mass balance here by substitution of <ref> or <ref>, together with Darcy's law, into <ref> such that  


    b_αϵt + (1/N_α+ϕ^0_α/K_l)p_αt + 1/Qp_βt + ∇·q_α = 1/ρ^0_lγ_α,


where we have used <ref> to decompose 1/M_α. 

Under the assumption of small perturbations in fluid density, the fourth term on the left-hand side of <ref> can be neglected. Comparing <ref> and <ref> leads to the following identity


    (1-ϕ_s)ψ_αt = 1/N_αp_αt + 1/Qp_βt


In the works of <cit.> and <cit.>, the authors achieve non-algebraic closure of the mass balance equations for each continuum by assuming ψ_αt≈0. This assumption is equivalent to


    1/N_α = 1/Q = 0   or  p_αt = -N_α/Qp_βt,


in <ref>. We discard the latter relation in <ref> 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 <cit.>, under the explicit assumption ψ_αt≈0, mapped to the constitutive model of <ref> and eqn:23eqn:24 is shown in <ref>.



  §.§.§ <cit.>
 
The motivation behind the approach  by <cit.> 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. Inclusion of intrinsic fracture properties is concurrent with the fracture continuum having an associated stiffness. 

The authors use a top down approach whose starting point is the macroscopic constitutive model written within a pure stiffness setting (σ and p_α as primary variables; contrary to conventional poromechanical modelling),


    [    ϵ; dξ_m; dξ_f ]
    =[ A_11 A_12 A_13; A_21 A_22 A_23; A_31 A_32 A_33 ][   dσ; dp_m; dp_f ],


where the goal is then to find expressions for coefficients A_11-A_33. Symmetry of the coefficient matrix in <ref> is implied. A comparison of <ref> with <ref> and eqn:23eqn:24 reveals the following relations


    K_dr = 1/A_11,      b_m = A_12/A_11,      b_f = A_13/A_11, 
    1/M_m = A_22-(A_12)^2/A_11,     1/Q = A_23-A_12A_13/A_11,     1/M_f = A_33-(A_13)^2/A_11.
  

To identify expressions for the parameters in <ref>, <cit.> considers scenarios of uniform expansion (or contraction). These scenearios are 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 () As a result the authors are able to relate <ref> to intrinsic constituent equations. The microscale equations used in <cit.> 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 (). In the case of a Gassmann-material, the constitutive equations for a material α can be written as   


    [  ϵ_α; dξ_α ] 
    = 1/K_α[               1           b^*_α;       v_α b^*_α v_α b^*_α/B^*_α ][ dσ_α; dp_α ],
       

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 


    v_m = Ω^0_m/Ω^0,     v_f =  Ω^0_f/Ω^0,


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


    b^*_α = 1 - K_α/K^α_s, 
    
    B^*_α = b^*_α M^*_α/K_α + (b^*_α)^2M^*_α  where  1/M^*_α = ϕ^*_α/K_l + b^*_α - ϕ^*_α/K^α_s.


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.  

Expressions for A_11-A_33 are finally recovered (<ref>), using the uniform expansion/contraction thought experiments described above. With the relations in <ref> we get material coefficient formulations pertaining to the conventional mixed compliance setting (ϵ and p_α as primary variables) (<ref>).

As a final note on the <cit.> coefficient models, recent users have explicitly assumed  A_23=A_32=0 as a closure condition to generalise the dual-continuum system to a multi-continuum one (; ; ), 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 <cit.> and <cit.> both make an underlying void space assumption for the high permeability, low storage continuum. Models from <cit.> make use of continuum pore fractions, but still require a final closure relationship for the evolution of each pore fraction. Finally, models from <cit.> 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 the starting point (or the requirement for constitutive expressions are at least identified in the case of <cit.>). The explicit decoupling assumptions are then made so as to provide relations between material properties and to simplify coefficient models (e.g. due to non-algebraic closure). However, the physical justifications and/or implications of these assumptions still remain an open question. 

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) (). 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.



 §.§ Effective elastic properties

We define effective properties over a representative elementary volume (REV) taken from a macroscopic body that is at least an order of magnitude greater in size than the REV itself. The scale requirements for the identification of such an REV can be found in <cit.>. In addition to the material isotropy assumption, we assume the dual-material to be statistically homogeneous. In its most simple interpretation, statistical homogeneity implies that averages taken over any REV from a large composite body, approach averages taken over the whole body itself. Averages are thus positionally invariant and we can recover the classical volume averaging expression accordingly   


    f = 1/Ω∫_Ωf(x) dV,
    

where f is an arbitrary field and Ω is the volume of any REV within a large macroscopic body.

With the averaging operation defined we proceed to formulate the effective elastic property problem. In doing we make use of the works of <cit.> and <cit.>. 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 


    s(x) = C_dr^*(x)e(x), 
    e(x) = S_dr^*(x)s(x).


It can be shown that macroscopic stress and strain tensors are related to their underlying fields through the volume averaging operator (see for example ) such that,


    σ = s = v_ms_m + v_fs_f, 
    ϵ = e = v_me_m + v_fe_f.


where e_α=Ω^-1_α∫_Ω_αe_α(x)dV (resp. s_α). Substitution of eqn:43eqn:44 in eqn:45eqn:46, and making use of the decompositions in eqn:14eqn:15 due to the isotropy of the dual-medium leads to


    σ = K_drϵ = v_mK_me_m + v_fK_fe_f  
    ϵ = S_drσ = v_mS_ms_m + v_fS_fs_f,


where e_α and s_α are the average microscopic volumetric strain and mean stress fields of continuum α respectively, and S_dr 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 () 


    e_m = A_mϵ,   e_f = A_fϵ, 
    s_m = B_mσ,   s_f = B_fσ,


where A_α and B_α are concentration factors mapping macroscopic volumetric strain and mean stress to the equivalent averaged microscopic fields. Importantly A_α and B_α encode the geometrical information of the problem.

Additionally, A_α (resp. B_α) admit the following compatibility relations 


    v_mA_m + v_fA_f = 1,    v_mB_m + v_fB_f = 1.


Use of eqn:49eqn:51 in <ref> and <ref> allows us to develop the following fundamental relationships for effective moduli of two phase composites,


    K_dr = K_m + v_f(K_f - K_m)A_f, 
    
    S_dr = S_m + v_f(S_f - S_m)B_f.


Finding effective moduli then amounts to determining A_α or B_α. 

For given underlying constituent properties, K_dr (resp. S_dr) can take on a range of values, existing between certain well-defined lower and upper bounds, due to the geometrical variability of the problem. As a result, the effective coefficients also exhibit a bounded range of values due to their dependence on K_dr (<ref>). This dependence is further explored in the following section.   



 §.§ Physical implications of explicit decoupling

In the following we show that explicit decoupling assumptions are coincident with effective coefficient bounds. Bounds are attractive as they provide useful estimates of effective properties of interest, as well as a means to verify the values of these properties ().    

We use arguments from micromechanics to show heuristically that the inverse to explicit decoupling is to assume that <ref> or <ref> can be calculated directly by considering certain limiting behaviours and thus bounds on K_dr. Bounds, from which explicit decouplings naturally arise, on effective properties then follow from the bounds on K_dr.     



  §.§.§ 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


    1/Q = ϕ_mp_f|_ϵ=0
    dp_m=0 = 0.


Next we define the following local constitutive model for a continuum α (),


    ds_α = K_αe_α - b^*_αdp_α, 
    dϕ_α = v_α(b^*_αe_α + 1/N^*_αdp_α).


Consider the local model given by eqn:55eqn:56 for the case α=m. Under a drained matrix (dp_m=0), <ref> 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 <ref> 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 <ref> and the macroscopic strain constraint, ϵ=0, the strain partition in <ref> shows that if ϵ=e_m=0 then e_f=0. The condition 1/Q=0 is therefore consistent with a condition of isostrain, 


    ϵ = e_m = e_f,    A_m = A_f = 1.


This distribution of strain can be obtained when elements are set in parallel to the direction of loading (<ref>). Under isostrain, the bulk modulus of the composite is then the Voigt average of the constituent moduli (), 


    K^V_dr = v_mK_m + v_fK_f.


<cit.> used variational principles to show that the Voigt average represents an upper bound on K_dr, and is naturally obtained by substitution of A_f=1 into <ref>. As a result, effective coefficients calculated with K^V_dr correspond to bounds on these parameters. To explore this further we consider the case for the effective Biot coefficients under the the void space assumption.    



When assuming a drained void fracture continuum (dp_f=0), and from <ref>, effective fracture strain is equal to the effective fracture pore strain, e_f=dϕ_f/ϕ^0_f, and thus in analogy to <ref> isostrain can be summarised in this case as


    ϵ = e_m = dϕ_f/ϕ^0_f.


Note that <ref>, and isostrain more generally, <ref>, are only equal to zero under the constrained macroscopic strain condition (ϵ=0). However, in general <ref> and <ref> are non-zero due to macroscopic deformation (ϵ≠0).

From <ref>, with drained matrix and fracture continua, together with the void-space isotrain condition, <ref>, we get the following lower bound on b_f 


    b_f = ϕ_fϵ|_dp_m=0
    dp_f=0 = ϕ^0_f.
   

From <ref> the lower bound on b_f corresponds to an upper bound for b_m 


    b_m = v_m (1 - K_m/K_s),


where v_m = 1 - ϕ^0_f, and where we have used K_dr = K^V_dr = v_mK_m. 

From the first equality in <ref> we expect 1/Q≤ 0 since matrix porosity must reduce in order to accommodate the pressure driven fracture expansion (see also similar arguments in <cit.>). Thus, based on the arguments described in this section we can infer that the explicit decoupling assumption 1/Q = 0 is concurrent with an upper bound on 1/Q. 



  §.§.§ Incompressible grain isostrain: 1/N_α=1/Q=0

We now consider the coefficient models from <cit.> under the assumption ψ_αt≈0 made by <cit.> and <cit.>. 

In Section 3.1.2, we indetified that ψ_αt≈0 amounts to 1/N_α=1/Q=0 when mapping to the constitutive model shown in eqn:17eqn:19. It is of interest to see under what conditions the result 1/N_α=1/Q=0 arises when starting from the void space coefficient models built from constituent mechanical properties (i.e. those from <cit.>). 

The set of explicit assumptions: 1/N_α=1/Q=0, can easily be derived from the microscale by first considering isostrain (and thus 1/Q=0). With the resulting bounds arising from isostrain, eqn:60eqn:61, along with the assumption K_s=∞, we obtain 1/N_α=0 using the coefficient models of <cit.> (<ref>, with <ref> to decompose 1/M_α). We therefore refer to conditions resulting in 1/N_α=1/Q=0 as incompressible grain isostrain

As far as parameters in the balance of mass are concerned, <ref> and <ref> are identical when assuming ψ_αt≈0 in the former and incompressible grain isostrain using the constituent mechanical property void space coefficient models in the latter. However, differences in mass balance behaviour may be introduced through the way in which b_α is modelled. It is therefore of interest see how the effective Biot coefficients calculated using the respective void space coefficient models under incompressible grain isostrain compare to the bounds established in eqn:60eqn:61. 

Under the incompressible grain assumption the upper bound for b_m now reads b_m=v_m. When using K^V_dr=v_mK_m in the coefficient models of <cit.> we can see that the bounds in eqn:60eqn:61 are naturally recovered (see <ref>). In contrast, from <ref> the effective Biot coefficients calculated using the models of <cit.> are equivalent to pore fractions, since b=1 for incompressible grains. 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: A_23(=A_32)=0
 
Finally, we study the pure stiffness setting with the condition A_23=0 in light of the assumptions made in <cit.>, <cit.>, <cit.> and <cit.>. We continue to work in terms of dϕ_m. Accordingly, and using similar arguments to those used in the derivations of eqn:23eqn:24, 


    A_23 = ξ_mp_f|_dσ=0
    dp_m=0≡ϕ_mp_f|_dσ=0
    dp_m=0 = 0.


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 <ref>, dσ=0, and the stress partition in <ref>, if dσ=ds_m=0 then ds_f=0. The condition A_23 is therefore consistent with a condition of isostress,


    dσ = ds_m = ds_f,    B_m = B_f = 1.


The classical configuration under which an isostress distribution is observed, are elements that are arranged transversely to the direction of applied load (<ref>). For isotropic fracture networks, a more frequent configuration that shows isostress behaviour is within a solid-fluid suspension (<ref>). 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 (), 


    1/K^R_dr = v_m/K_m + v_f/K_f.
 

<cit.> showed that the Reuss average is a lower bound on K_dr, and is naturally obtained by substitution of B_f=1 into <ref>, or from A_23=0 in <ref>. In analogy to K^V_dr, use of K^R_dr for the calculation of the effective constitutive parameters will result in bounds on these parameters.    

Consider the case when K_f is zero (void space fracture phase) we can see from <ref> that K_dr is also zero. Using the void space assumption and the isostress condition, <ref>, with <ref> for dσ and <ref> for dσ_m and dσ_f gives us the following 


    b_mdp_m - b_fdp_f = K_me_m - b^*_mdp_m = -1dp_f.


From the required isostress equality, σ = s_f, <ref> followed by <ref> allows us to establish the following bounds 


    b_f = 1,    b_m = 0.
   

<ref> shows that use of the Reuss bound corresponds to an upper bound on b_f and a lower bound on b_m. Interestingly the bounds on b_f from this section and Section 4.2.1, i.e. ϕ^0_f ≤ b_f ≤ 1 are very similar to those established on b for the single-porosity model (see <cit.>). 



 §.§ On moduli upscaling
    
In <cit.>, <cit.>, <cit.>, <cit.> and more recently <cit.>, 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 (). Moduli bounds represent a minimum and maximum limit, between, or at which, effective moduli arise. For the bulk modulus, bounds by <cit.> 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


    K^HS-_dr = K_f + v_m/(K_m - K_f)^-1 + v_f(K_f + 4/3G_f)^-1,
     

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 <ref>.

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^R_dr=K^HS+_dr=0. However, in situations when the fracture phase has an intrinsic stiffness, the Reuss and HS lower bounds may be significantly different (). 

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 (). The effective bulk modulus of the composite will then coincide with the HS lower bound (; ). In the general case fractures are likely to have an associated stiffness (). 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 (). 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. 
Comprehensive summaries of such approaches can be found in the works of <cit.>, <cit.> and <cit.> for example.
               


§ 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, <ref> and eqn:23eqn:24, is the starting point but pressure coupling is neglected leading to bounds being used for the calculation of K_dr (and thus the calculation of the effective constitutive coefficients).

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 <cit.> and <cit.> 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) (<ref>). 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= -2Fa, 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,


    σ_xx = σ_xz = p_m = p_f = 0      on Γ_E ∪Γ_W, 
    σ_xz = q_m,z = q_f,z = 0      on Γ_N ∪Γ_S, 
    ∫^a_-aσ_zz dx = -2Fa      on Γ_N ∪Γ_S.




Due to the symmetry of the problem only a quarter of the domain need be considered.

For an isotropic double-porosity material <ref> and <ref> can be re-written and extended to


    dσ = 2Gϵ + λϵ1 -∑_α=m,fb_αdp_α1


where λ=2Gv/1-2v is the Lamé constant, in which v is Poisson's ratio. Substitution of <ref> and <ref> in <ref> yields  


    G∇^2u + (λ + G)∇ϵ = ∑_α=m,fb_α∇ p_α1 - ρg,
   

where we assume γ≈0. It can be shown that ∇^2u=∇ϵ, with which, and in the absence of body forces, <ref> reduces to  


    ∇ϵ = ∑_α=m,fc_α∇ p_α,


where c_α=b_α(1-2v)/2G(1-v) is the consolidation coefficient belonging to material α.

Integration of <ref> leads to 


    ϵ = ∑_i=m,f c_α p_α + f(t),


where f(t) is an integration function. Use of <ref> in <ref> leads to a set of diffusion equations written entirely in terms of continuum pressures. Solutions to the resulting system of equations are presented in <cit.>, and <cit.>. Further, these solutions allow for calculation of vertical stress and strain (see  or ).    



 §.§ 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 (). 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 (). Effective permeability of the fracture is assumed to be k_f≈1000  md. With the cubic law (), 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. In each case, analytical solutions to the double-porosity Mandel problem are compared. Differences in solutions for each case then arise due the parameter permutations described in the case descriptions that follow.  



  §.§.§ 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 <cit.> could still be used even when a fracture has an associated phase stiffness.  

From tab:1tab:2, we hypothesise that provided that ϕ^*_f≈1 and the fracture phase stiffness is orders of magnitude lower than the grain stiffness (K_f≪ K^f_s), effects arising due to deviations of intrinsic fracture coefficients under the void space assumption (i.e. deviations from: b^*_f=1, 1/M^*_α = 1/K_l and B^*_α = 1) are  negligible. Void space coefficient models could then be used in place of coefficient models that include intrinsic fracture properties. When ϕ^*_f=1 we use coefficient models from <cit.>, and when ϕ^*_f<1 we use coefficient models from <cit.> with K^f_s=K_s. 

Strictly speaking, use of void space coefficient models implies K_f=0. However, the aim of this test case is to highlight the effect of missing physics by not including intrinsic poromechanical parameters within coefficient model formulations. To do this, and to ensure non-zero values of K_dr, we use the same upscaled bulk modulus for both sets of coefficient models, which is calculated with a non-zero K_f. 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. 

We compare results from coefficient models calculated using a fracture modulus several orders of magnitude lower than the solid grain modulus, (K^†_f=K^f_s/1750), versus ones that use a fracture modulus only an order of magnitude lower than the solid modulus, (K^†_f=K^f_s/35). 



  §.§.§ 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 <cit.> and <cit.>. 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 <cit.> and <cit.>, with K_f=0, respectively. We reference these results against ones for which no decoupling is made with coefficient models from <cit.>.

Use of the void space models implies K_f=0. Since we do not enforce relations on K_dr (a fundamental difference between the implicit assumptions considered here and explicit assumptions) we take an arbitrary value of K^†_dr=10  GPa. As a result bounds on the effective constitutive coefficients are not enforced.  



  §.§.§ Case 3: Explicit decoupling assumption - isostrain

We investigate the effect of assuming isostrain at the microscale whilst making use of coefficient models from <cit.>. Under isostrain the composite and constituent bulk moduli are linked by the Voigt average (K^V_dr= 19.8  GPa).  This leads naturally to bounds on the effective constitutive coefficients, with 1/Q=0 representing an upper bound. 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^HS+_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 <cit.> under the assumption of ψ_αt=0, and <cit.> under the assumption of incompressible grain isostrain. We use K^V_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 <cit.> have been used with an explicit decoupling assumption (A_23=0) that implies isostress (; ; ). 

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=K_m/500. Coefficient models from <cit.> are then used. 

We compare results arising when calculating the composite bulk modulus with the Reuss average (K^R_dr=3.3  GPa), the HS lower bound (K^HS-_dr=5.7  GPa), and an arithmetic average of the HS bounds (K^AHS_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 double-porosity 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 <ref>. 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, <ref>. 

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 (). Consequently, each continuum is undrained, m^t(0^+)_l,α=m^t(0)_l,α. From <ref> and <ref>, together with the undrained condition (dm_l,α=0) we can recover


    dp_α = (1/M_α - M_β/(Q)^2)^-1(M_β b_β/Q-b_α)ϵ.
 

We note that ϵ_zz∝ϵ. Substitution of <ref> into <ref> leads to   


    dp_α = -M_α b_α/{K_dr + M_α(b_α)^2}{dσ + ( b_β +K_dr/b_α Q) dp_β}.
   



 §.§ Case 1: Intrinsic fracture properties

<ref> and <ref> 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 <cit.> and <cit.> when the fracture is almost all void space (ϕ^*_f ≈ 1) and the fracture phase stiffness is orders of magnitude lower than the solid grain stiffness (K_f≪ K^f_s) (<ref>). 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 (<ref>). 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 (<ref>). This phenomenon can be explained by considering the intrinsic Skempton coefficient, B^*_α, written in the following form for the fracture phase (),


    B^*_f = 1 - ϕ^*_fK_f(K^f_s-K_l)/K_l(K^f_s-K_f)+ϕ^*_fK_f(K^f_s-K_l).
 

Use of <ref> allows us to cast our observations as a bounding problem such that


    B^*_f(ϕ^*_f=1,K_f) < B^*_f(ϕ^*_f<1,K_f) ≤ 1.
 

The lower bound in <ref> is a fictitious one in that technically materials with an intrinsic porosity of one should also have a zero stiffness. However, accepting this contradiction is useful in approaching the bounding problem from a purely quantitative point of view. For the upper bound, <ref> shows that as ϕ^*_f approaches zero, B^*_f must asymptotically approach one.

If the lower bound in <ref> 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 (<ref>).  

<ref> and <ref> show the variation in vertical strain for the softer and stiffer fracture phases respectively. In <ref>, the strain evolutions are almost identical when using the more compliant fracture phase for the range of intrinsic fracture porosities. This observation is coupled to the similarity in pressure evolutions. Whilst induced fracture pressures are significantly larger when fracture stiffness and intrinsic porosity are non-negligible, <ref> 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 <ref> 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 <ref>. When assuming A_23=0, the matrix and fracture pressure evolutions are almost identical to the reference case (<ref>), with vertical strains also being correspondingly very similar (<ref>). When assuming 1/Q=0, the early time matrix and fracture pressures are measurably lower than the reference case (<ref>). The early time vertical strain is greater when 1/Q=0 than the strain in the reference case (<ref>).    

The results in <ref> suggest that the assumption 1/Q=0 has the most noticeable effect on the poromechanical behaviour of the dual-medium. As discussed in Section 4.2.1 we would expect 1/Q<0. From <ref>, 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 <ref> we can see dp_α∝ϵ^-1. Underestimated pressures therefore explain the over estimated strain when taking 1/Q=0.

In contrast, <ref> 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 bounds on bulk moduli. The remainder of this results section aims at qualitatively supporting these findings.  



 §.§ Case 3: Explicit decoupling - isostrain

<ref> 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 <ref>. 



An alternative, heuristic approach to explaining the pressure distributions in <ref> 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. 

The early time vertical strains in <ref> can be explained by way of <ref> 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 eqn:60eqn:61. 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, <ref>. This explains the growth in vertical strain separation at later times shown in <ref>.     

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. 

<ref> displays the results for the incompressible grain isostrain investigation. For the case of coefficient models from <cit.>, the induced matrix pressure is significantly larger than the induced fracture pressure. Conversely, when using coefficient models from <cit.> with the assumption ψ_αt≈ 0, induced matrix and fracture pressures are equal. This can be explained by considering <ref> which allows us to equate induced variations in matrix and fracture pressures such that


    (1/M_m - M_f/(Q)^2)    (M_fb_f/Q-b_m)^-1dp_m 
       = (1/M_f - M_m/(Q)^2) (M_mb_m/Q-b_f)^-1dp_f.


Assuming ψ_αt≈ 0 together with the effective Biot coefficient expressions from <cit.> (<ref>), <ref> reduces to


    (ϕ^0_m/K_l1-ϕ_s/ϕ^0_m)dp_m =
    (ϕ^0_f/K_l1-ϕ_s/ϕ^0_f)dp_f,


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. 

<ref> shows vertical strain is lower at early times when using coefficient models from <cit.> under incompressible grain isostrain. This can be explained by <ref>, 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 eqn:17eqn:19.
  
  


 §.§ Case 4: Explicit decoupling - isostress

<ref> shows further impacts 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.  

Non-monotonic pressure rises are often referred to as the Mandel-Cryer effect within literature (; ). 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. 

 

<ref> 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.       
 
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, <ref>.   

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 goal of this paper was to formulate a set of recommendations for how and when to use different constitutive modelling concepts. In doing we revisited and reviewed three main modelling approaches, which differ in the parameters, and/or assumptions, used for the calculation of the effective constitutive coefficients. The parameters and assumptions used in these approaches are summarised as follows: 



  * Constituent mechanical properties; assuming the high permeability, low storage continuum is all void space (no intrinsic fracture properties),

  * Constituent pore fractions; assuming the high permeability, low storage continuum is all void space,

  * Constituent mechanical properties, including intrinsic fracture properties.


Based on theoretical and qualitative findings in this paper we recommend further work on algebraic closure conditions for models built using continuum pore fractions. Comparing coefficient models (i) and (iii) we found that the effects of intrinsic fracture properties become measurable when there are significant deviations from the intrinsic poromechanical constitutive coefficients of a void space fracture continuum. In this case ϕ^*_f<1 and K_f≪̸K^f_s, and it is advisable to use coefficient models where intrinsic fracture properties are naturally incorporated. We envisage the aforementioned conditions, and thus real benefit of using models (iii), to be observed when considering non-linear poromechanics, where the internal structure of the high permeability, low storage continuum is evolving (e.g. fracture closure) thus leading to cases where ϕ^*_f<1 and K_f≪̸K^f_s. However, in the linear case, our results show models of type (i) to give very good matches to models of type (iii) even when ϕ^*_f<1, provided K_f≪ K^f_s. Therefore as a first approach we recommend the use of models (i) for poromechanical dual-continuum modelling given that we expect the high permeability, low storage continuum to be mechanically weaker than solid grains.     

The second set of recommendations is formed on the basis of our investigations into implicit and explicit decoupling assumptions. In both cases mechanical coupling between continuum pressures is neglected. In the former we showed that implicit assumptions can lead to the removal of pressure sources leading to physically unreliable results. We therefore recommend the use of a full constitutive system where possible.  

Even with a full constitutive system, explicit assumptions have been made as a passage to simplifying relations between composite and constituent moduli without considering the physical implications of their use. In this case we showed that explicit decoupling assumptions are coincident with bounds on composite moduli that arise naturally under end-member states of isostrain and isostress. However, for isotropic composite materials it is well known that the bounds obtained under isotrain and isostress can be loose, and that tighter bounds using similar quantities are readily available within literature. Our qualitative investigations showed clear differences in poromechanical behaviour when using these different bounds. 

To conclude, bounds arising from isostrain and isostress states, which are concurrent with explicit decoupling assumptions, can provide a useful means for guiding our intuition into multiscale poromechanical behaviour, given their ease of computation. However, for practical subsurface applications, we recommend against the use of explicit decoupling assumptions, as they have physical and geometrical implications that are unlikely to be justified within isotropic multiscale materials.                     



 §.§ Acknowledgements

The authors are very grateful for the funding provided to them by the National Environmental Research Council to carry out this work, as well as to the reviewers who helped to improve this work.  

agsm
xx

Aboudi1992Aboudi1992
Aboudi, J.  1992, Mechanics of
  composite materials: a unified micromechanical approach, Elsevier.

Aifantis1977Aifantis1977
Aifantis, E. C.  1977, `Introducing a
  multi-porous medium', Dev. Mech 8(3), 209–211.

Aifantis1979Aifantis1979
Aifantis, E. C.  1979, `On the response of
  fissured rocks', Dev. Mech 10, 249–253.

[Alberto et al.]Alberto, Cordero, Mejia   Roehl2019Alberto2019
Alberto, J., Cordero, R., Mejia, E. C.  Roehl, D.  2019, `Dual permeability models for fluid flow in
  deformable fractured media', Journal of Petroleum Science and
  Engineering 175(November 2018), 644–653.

Arbogast  Douglas1990Arbogast1990
Arbogast, T.  Douglas, J. I. M.  1990, `Derivation of the double porosity model of single
  phase flow via homogenization theory',  21(4), 823–836.

[Bai et al.]Bai, Meng, Elsworth, Abousleiman   Roegiers1999Bai1999
Bai, M., Meng, F., Elsworth, D., Abousleiman, Y.  Roegiers, J.
  1999, `Numerical modelling of coupled
  flow and deformation in fractured rock specimens', International
  journal for numerical and analytical methods in geomechanics 23(2), 141–160.

[Bandis et al.]Bandis, Lumsden   Barton1983Bandis1983
Bandis, S., Lumsden, A.  Barton, N.  1983, Fundamentals of rock joint deformation, in
  `International Journal of Rock Mechanics and Mining Sciences & Geomechanics
  Abstracts', Vol. 20, Elsevier, pp. 249–268.

[Barenblatt et al.]Barenblatt, Zheltov   Kochina1960Barenblatt1960
Barenblatt, G., Zheltov, I.  Kochina, I.  1960, `Basic concepts in the theory of seepage of
  homogeneous liquids in fissured rocks [strata]', Journal of Applied
  Mathematics and Mechanics 24(5), 1286–1303.

Bear  Bachmat2012Bear2012
Bear, J.  Bachmat, Y.  2012,
  Introduction to modeling of transport phenomena in porous media,
  Vol. 4, Springer Science & Business Media.

Beran1971Beran1971
Beran, M.  1971, `Application of statistical
  theories to heterogeneous materials', physica status solidi (a) 6(2), 365–384.

Berkowitz2002Berkowitz2002
Berkowitz, B.  2002, `Characterizing flow
  and transport in fractured geological media: A review', Advances in
  Water Resources 25(8-12), 861–884.

[Berre et al.]Berre, Doster   Keilegavlen2018Berre2018
Berre, I., Doster, F.  Keilegavlen, E.  2018, `Flow in fractured porous media: A review of
  conceptual models and discretization approaches', Transport in Porous
  Media .

Berryman2002Berryman2002a
Berryman, J. G.  2002, `Extension of
  Poroelastic Analysis to Double-Porosity Materials: New Technique in
  Microgeomechanics', Journal of Engineering Mechanics 128(8), 840–847.

Berryman2006Berryman2006
Berryman, J. G.  2006, `Effective medium
  theories for multicomponent poroelastic composites', Journal of
  Engineering Mechanics 132(5), 519–531.

Berryman  Pride2002Berryman2002
Berryman, J. G.  Pride, S. R.  2002, `Models for computing geomechanical constants of
  double‐porosity materials from the constituents' properties', Journal
  of Geophysical Research: Solid Earth 107(B3).

Berryman  Wang1995Berryman1995
Berryman, J. G.  Wang, H. F.  1995, `The elastic coefficients of double‐porosity
  models for fluid transport in jointed rock', Journal of Geophysical
  Research: Solid Earth 100(B12), 24611–24627.

Biot1941Biot1941
Biot, M. A.  1941, `General theory of
  three‐dimensional consolidation', Journal of Applied Physics 12(2), 155–164.

Biot1977Biot1977
Biot, M. A.  1977, `Variational
  Lagrangian-thermodynamics of nonisothermal finite strain mechanics of porous
  solids and thermomolecular diffusion', International Journal of Solids
  and Structures 13(6), 579–597.

[Blessent et al.]Blessent, Jørgensen   Therrien2014Blessent2014
Blessent, D., Jørgensen, P. R.  Therrien, R.  2014, `Comparing discrete fracture and continuum models to
  predict contaminant transport in fractured porous media', Groundwater
  52(1), 84–95.

Borja  Choo2016Borja2016
Borja, R. I.  Choo, J.  2016,
  `Cam-clay plasticity, part viii: A constitutive framework for porous
  materials with evolving internal structure', Computer Methods in Applied
  Mechanics and Engineering 309, 653–679.

Borja  Koliji2009Borja2009
Borja, R. I.  Koliji, A.  2009,
  `On the effective stress in unsaturated porous continua with double
  porosity', Journal of the Mechanics and Physics of Solids 57(8), 1182–1193.

Boucher1974Boucher1974
Boucher, S.  1974, `On the effective moduli
  of isotropic two-phase elastic composites', Journal of Composite
  Materials 8(1), 82–89.

Budiansky  O'connell1976Budiansky1976
Budiansky, B.  O'connell, R. J.  1976, `Elastic moduli of a cracked solid', International journal of Solids and structures 12(2), 81–97.

Callari  Federico2000Callari2000
Callari, C.  Federico, F.  2000,
  `FEM validation of a double porosity elastic model for consolidation of
  structurally complex clayey soils', International Journal for Numerical
  and Analytical Methods in Geomechanics 24(4), 367–402.

Cheng2016Cheng2016
Cheng, A. H.-D.  2016, Poroelasticity, Vol. 27, Springer.

Choo  Borja2015Choo2015
Choo, J.  Borja, R. I.  2015,
  `Stabilized mixed finite elements for deformable porous media with double
  porosity', Computer Methods in Applied Mechanics and Engineering 293, 131–154.

[Choo et al.]Choo, White  Borja2016Choo2016
Choo, J., White, J. A.  Borja, R. I.  2016, `Hydromechanical modeling of unsaturated flow in
  double porosity media', International Journal of Geomechanics 16(6), D4016002.

Coussy1995Coussy1995
Coussy, O.  1995, Mechanics of porous
  continua, Wiley.

Coussy2004Coussy2004
Coussy, O.  2004, Poromechanics,
  John Wiley & Sons.

Coussy2005Coussy2005
Coussy, O.  2005, `Poromechanics of
  freezing materials', Journal of the Mechanics and Physics of Solids
  53(8), 1689–1718.

Coussy2007Coussy2007
Coussy, O.  2007, `Revisiting the
  constitutive equations of unsaturated porous solids using a Lagrangian
  saturation concept', International Journal for Numerical and Analytical
  Methods in Geomechanics 31(15), 1675–1694.

Coussy2011Coussy2011
Coussy, O.  2011, Mechanics and
  physics of porous solids, John Wiley & Sons.

De Boer2012DeBoer2012
De Boer, R.  2012, Theory of porous
  media: highlights in historical development and current state, Springer
  Science & Business Media.

Detournay  Cheng1995Detournay1995
Detournay, E.  Cheng, A. H.-D.  1995, Fundamentals of poroelasticity, in `Analysis
  and design methods', Elsevier, pp. 113–171.

[Dormieux et al.]Dormieux, Kondo   Ulm2006Dormieux2006
Dormieux, L., Kondo, D.  Ulm, F.-J.  2006, Microporomechanics, John Wiley & Sons.

Elsworth  Bai1992Elsworth1992
Elsworth, D.  Bai, M.  1992,
  `Flow-deformation response of dual-porosity media', Journal of
  Geotechnical Engineering 118(1), 107–124.

Eshelby1957Eshelby1957
Eshelby, J. D.  1957, `The determination of
  the elastic field of an ellipsoidal inclusion, and related problems', Proc. R. Soc. Lond. A 241(1226), 376–396.

[Fornells et al.]Fornells, García-Aznar   Doblaré2007Fornells2007
Fornells, P., García-Aznar, J. M.  Doblaré, M.
  2007, `A finite element dual porosity
  approach to model deformation-induced fluid flow in cortical bone', Annals of Biomedical Engineering 35(10), 1687–1698.

[Garipov et al.]Garipov, Karimi-Fard   Tchelepi2016Garipov2016
Garipov, T. T., Karimi-Fard, M.  Tchelepi, H. A.  2016, `Discrete fracture model for coupled flow and
  geomechanics', Computational Geosciences 20(1), 149–160.

Gerke2006Gerke2006
Gerke, H. H.  2006, `Preferential flow
  descriptions for structured soils', Journal of Plant Nutrition and Soil
  Science 169(3), 382–400.

Gong2007Gong2007
Gong, B.  2007, Effective models of
  fractured systems, PhD thesis, Stanford University.

Hajiabadi  Khoei2019Hajiabadi2019
Hajiabadi, M. R.  Khoei, A. R.  2019, `A bridge between dual porosity and multiscale
  models of heterogeneous deformable porous media',  (July 2018), 212–238.

Hashin1972Hashin1972
Hashin, Z.  1972, Theory of fiber
  reinforced materials, NASA CR-1974.

Hashin1983Hashin1983
Hashin, Z.  1983, `Analysis of composite
  materials—a survey', Journal of Applied Mechanics 50(3), 481–505.

Hashin  Shtrikman1963Hashin1963
Hashin, Z.  Shtrikman, S.  1963,
  `A variational approach to the theory of the elastic behaviour of multiphase
  materials', Journal of the Mechanics and Physics of Solids 11(2), 127–140.

Hill1952Hill1952
Hill, R.  1952, `The elastic behaviour of a
  crystalline aggregate', Proceedings of the Physical Society. Section A
  65(5), 349.

Hill1963Hill1963
Hill, R.  1963, `Elastic properties of
  reinforced solids: some theoretical principles', Journal of the
  Mechanics and Physics of Solids 11(5), 357–372.

Hill1965Hill1965
Hill, R.  1965, `A self-consistent mechanics
  of composite materials', Journal of the Mechanics and Physics of Solids
  13(4), 213–222.

[Karimi-Fard et al.]Karimi-Fard, Gong   Durlofsky2006Karimi-Fard2006
Karimi-Fard, M., Gong, B.  Durlofsky, L. J.  2006, `Generation of coarse-scale continuum flow models
  from detailed fracture characterizations', Water Resources Research
  42(10), 1–13.

[Kazemi et al.]Kazemi, Merrill Jr, Porterfield   Zeman1976Kazemi1976
Kazemi, H., Merrill Jr, L. S., Porterfield, K. L.  Zeman, P. R.
  1976, `Numerical simulation of water-oil
  flow in naturally fractured reservoirs', Society of Petroleum Engineers
  Journal 16(06), 317–326.

Khalili2003Khalili2003
Khalili, N.  2003, `Coupling effects in
  double porosity media with deformable matrix', Geophysical Research
  Letters 30(22).

Khalili2008Khalili2008
Khalili, N.  2008, `Two‐phase fluid flow
  through fractured porous media with deformable matrix', Water Resources
  Research 44(5).

[Khalili et al.]Khalili, Khabbaz   Valliappan2000Khalili2000
Khalili, N., Khabbaz, M. H.  Valliappan, S.  2000, `An effective stress based numerical model for
  hydro-mechanical analysis in unsaturated porous media', Computational
  Mechanics 26(2), 174–184.

Khalili  Selvadurai2003Khalili2003a
Khalili, N.  Selvadurai, A. P. S.  2003, `A fully coupled constitutive model for
  thermo-hydro-mechanical analysis in elastic media with double porosity',
  Geophysical Research Letters 30(24), 1–5.

Khalili  Valliappan1996Khalili1996
Khalili, N.  Valliappan, S.  1996, `Unified theory of flow and deformantion in double porous media', European journal of mechanics. A. Solids 15(2), 321–336.

Kim  Moridis2013Kim2013
Kim, J.  Moridis, G. J.  2013,
  `Development of the t+ m coupled flow–geomechanical simulator to describe
  fracture propagation and coupled flow–thermal–geomechanical processes in
  tight/shale gas systems', Computers & geosciences 60, 184–198.

[Kim et al.]Kim, Sonnenthal  Rutqvist2012Kim2012
Kim, J., Sonnenthal, E. L.  Rutqvist, J.  2012, `Formulation and sequential numerical algorithms of
  coupled fluid/heat flow and geomechanics for multiple porosity materials',
  International Journal for Numerical Methods in Engineering 92(5), 425–456.

Koliji2008Koliji2008
Koliji, A.  2008, Mechanical behaviour
  of unsaturated aggregated soils, PhD thesis, EPFL.

[Levin et al.]Levin, Kanaun  Markov2012Levin2012
Levin, V., Kanaun, S.  Markov, M.  2012, `Generalized Maxwell's scheme for homogenization of
  poroelastic composites', International Journal of Engineering Science
  61, 75–86.

Lewis  Ghafouri1997Lewis1997
Lewis, R. W.  Ghafouri, H. R.  1997, `A novel finite element double porosity model for
  multiphase flow through deformable fractured porous media', International Journal for Numerical and Analytical Methods in Geomechanics
  21(11), 789–816.

Lim  Aziz1995Lim1995
Lim, K. T.  Aziz, K.  1995,
  `Matrix-fracture transfer shape factors for dual-porosity simulators', Journal of Petroleum Science and Engineering 13(3-4), 169–178.

Loret  Rizzi1999Loret1999
Loret, B.  Rizzi, E.  1999,
  `Strain localization in fluid-saturated anisotropic elastic-plastic porous
  media with double porosity', Journal of the Mechanics and Physics of
  Solids 47(3), 503–530.

Mehrabian2018Mehrabian2018
Mehrabian, A.  2018, `The poroelastic
  constants of multiple-porosity solids', International Journal of
  Engineering Science 132, 97–104.

Mehrabian  Abousleiman2014Mehrabian2014
Mehrabian, A.  Abousleiman, Y. N.  2014, `Generalized biot's theory and mandel's problem of
  multiple-porosity and multiple-permeability poroelasticity', Journal of
  Geophysical Research: Solid Earth 119(4), 2745–2763.

Mehrabian  Abousleiman2015Mehrabian2015
Mehrabian, A.  Abousleiman, Y. N.  2015, `Gassmann equations and the constitutive relations
  for multiple-porosity and multiple-permeability poroelasticity with
  applications to oil and gas shale', International Journal for Numerical
  and Analytical Methods in Geomechanics 39(14), 1547–1569.

Nelson2001Nelson2001
Nelson, R.  2001, Geologic analysis of
  naturally fractured reservoirs, Gulf Professional Publishing.

Nemat-Nasser  Hori1993Nemat1993
Nemat-Nasser, S.  Hori, M.  1993, Micromechanics: overall properties of heterogeneous materials,
  Elsevier.

Nguyen  Abousleiman2010Nguyen2010
Nguyen, V. X.  Abousleiman, Y. N.  2010, `Poromechanics solutions to plane strain and
  axisymmetric mandel-type problems in dual-porosity and dual-permeability
  medium', Journal of Applied Mechanics 77(1), 011002.

Norris1985Norris1985
Norris, A.  1985, `A differential scheme for
  the effective moduli of composites', Mechanics of materials 4(1), 1–16.

Pao  Lewis2002Pao2002
Pao, W. K.  Lewis, R. W.  2002,
  `Three-dimensional finite element simulation of three-phase flow in a
  deforming fissured reservoir', Computer Methods in Applied Mechanics
  and Engineering 191(23-24), 2631–2659.

Reuss1929Reuss1929
Reuss, A.  1929, `Berechnung der
  fließgrenze von mischkristallen auf grund der plastizitätsbedingung
  für einkristalle.', ZAMM-Journal of Applied Mathematics and
  Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 9(1), 49–58.

Rice1997Rice1997
Rice, R.  1997, `Limitations of pore-stress
  concentrations on the mechanical properties of porous materials', Journal of Materials Science 32(17), 4731–4736.

[Romero et al.]Romero, Della Vecchia   Jommi2011Romero2011
Romero, E., Della Vecchia, G.  Jommi, C.  2011, `An insight into the water retention properties of
  compacted clayey soils', Géotechnique 61(4), 313.

Rutqvist  Stephansson2003Rutqvist2003
Rutqvist, J.  Stephansson, O.  2003, `The role of hydromechanical coupling in fractured
  rock engineering', Hydrogeology Journal 11(1), 7–40.

[Taron et al.]Taron, Elsworth  Min2009Taron2009
Taron, J., Elsworth, D.  Min, K. B.  2009, `Numerical simulation of
  thermal-hydrologic-mechanical-chemical processes in deformable, fractured
  porous media', International Journal of Rock Mechanics and Mining
  Sciences 46(5), 842–854.

Torquato1991Torquato1991
Torquato, S.  1991, `Random heterogeneous
  media: microstructure and improved bounds on effective properties', Applied mechanics reviews 44(2), 37–76.

Torquato2002Torquato2002
Torquato, S.  2002, Random
  heterogeneous materials: microstructure and macroscopic properties, Springer
  Science & Business Media.

Tuncay  Corapcioglu1995Tuncay1995
Tuncay, K.  Corapcioglu, M. Y.  1995, `Effective stress principle for saturated fractured
  porous media', Water Resources Research 31(12), 3103–3106.

Tuncay  Corapcioglu1996Tuncay1996a
Tuncay, K.  Corapcioglu, M. Y.  1996, `Wave propagation in fractured porous media', Transport in Porous Media 23(3), 237–258.

Voigt1928Voigt1928
Voigt, W.  1928, Lehrbuch der
  Kristallphysik (mit Ausschluss der Kristalloptik), Springer-Verlag.

Wang2000Wang2000
Wang, H. F.  2000, Theory of linear
  poroelasticity with applications to geomechanics and hydrogeology, Princeton
  University Press.

Wang  Sun2018Wang2018
Wang, K.  Sun, W.  2018, `A
  multiscale multi-permeability poroplasticity model linked by recursive
  homogenizations and deep learning', Computer Methods in Applied
  Mechanics and Engineering 334, 337–380.

[Warren et al.]Warren, Root et al.1963Warren1963
Warren, J., Root, P. J. et al.  1963, `The
  behavior of naturally fractured reservoirs', Society of Petroleum
  Engineers Journal 3(03), 245–255.

[Watt et al.]Watt, Davies  O'Connell1976Watt1976
Watt, J. P., Davies, G. F.  O'Connell, R. J.  1976, `The elastic properties of composite materials', Reviews of Geophysics 14(4), 541–563.

Wilson  Aifantis1982Wilson1982
Wilson, R. K.  Aifantis, E. C.  1982, `On the theory of consolidation with double
  porosity', International Journal of Engineering Science 20(9), 1009–1035.

[Witherspoon et al.]Witherspoon, Wang, Iwai   Gale1980Witherspoon1980
Witherspoon, P. A., Wang, J. S., Iwai, K.  Gale, J. E.
  1980, `Validity of cubic law for fluid
  flow in a deformable rock fracture', Water Resources Research 16(6), 1016–1024.


 


