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

TL;DR

This paper reviews and compares various models for constitutive coefficients in dual-continuum poromechanical systems, highlighting the effects of assumptions and bounds on system behavior and providing practical recommendations.

Contribution

It provides a comprehensive comparison of existing coefficient models, analyzes the implications of decoupling assumptions, and offers guidance on their appropriate application.

Findings

Explicit decouplings align with bounds on effective parameters.

Decoupling assumptions relate to end-member states like isostress or isostrain.

Using bounds and assumptions significantly impacts poromechanical behavior.

Abstract

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

Tables

Table 1. Table 1: poromechanical material cofficients by author. + superscript \text{}^{+} Assuming \pderiv ​ ψ α ​ t ≈ 0 \pderiv subscript 𝜓 𝛼 𝑡 0 \pderiv{\psi_{\alpha}}{t}\approx 0 . + + superscript absent \text{}^{++} Coefficients A 22 subscript 𝐴 22 A_{22} , A 23 subscript 𝐴 23 A_{23} and A 33 subscript 𝐴 33 A_{33} are given in Table   2 . 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 subscript 𝑏 𝑓 b_{f} and 1 M f 1 subscript 𝑀 𝑓 \frac{1}{M_{f}} and single porosity equivalents, b 𝑏 b and 1 M 1 𝑀 \frac{1}{M} (see for example Coussy ( 2004 ) ).

	$b_f$	$b_m$	$\frac{1}{M_f}$	$\frac{1}{Q}$	$\frac{1}{M_m}$
Khalili and Valliappan (1996)	$1-\frac{K_{dr}}{K_m}$	$\frac{K_{dr}}{K_m}-\frac{K_{dr}}{K_s}$	$\frac{{\varphi }_f^0}{K_l}+\frac{b_f-{\varphi }_f^0}{K_m}$	$\frac{b_f-{\varphi }_f^0}{K_s}-\frac{b_f-{\varphi }_f^0}{K_m}$	$\frac{{\varphi }_m^0}{K_l}+\frac{b_m-{\varphi }_m^0}{K_s}-\frac{1}{Q}$
Borja and Koliji (2009)${}^{+}$	${\psi }_f^{0}b$	${\psi }_m^{0}b$	$\frac{{\varphi }_f^0}{K_l}$	$0$	$\frac{{\varphi }_m^0}{K_l}$
Berryman (2002)${}^{++}$	$b_f^{\ast }\frac{K_{dr}-K_m}{K_f-K_m}$	$b_m^{\ast }\frac{K_{dr}-K_f}{K_m-K_f}$	$A_{33}-\frac{{\left(b_f\right)}^2}{K_{dr}}$	$A_{23}-\frac{b_m{b}_f}{K_{dr}}$	$A_{22}-\frac{{\left(b_m\right)}^2}{K_{dr}}$

Table 2. Table 2: Berryman ( 2002 ) material coefficient formulations for the pure stiffness constitutive model.

Coefficient	Berryman (2002) Formulation
$A_{11}$	$\frac{1}{K_{dr}}$
$A_{12}$	$\frac{b_m^{\ast }}{K_m}\frac{1-K_f/K_{dr}}{1-K_f/K_m}$
$A_{13}$	$\frac{b_f^{\ast }}{K_f}\frac{1-K_m/K_{dr}}{1-K_m/K_f}$
$A_{22}$	$\frac{v_m{b}_m^{\ast }}{B_m^{\ast }{K}_m}-{\left(\frac{b_m^{\ast }}{1-K_m/K_f}\right)}^2\left\{\frac{v_m}{K_m}+\frac{v_f}{K_f}-\frac{1}{K_{dr}}\right\}$
$A_{23}$	$\frac{K_m{K}_f{b}_m^{\ast }{b}_f^{\ast }}{{\left(K_f-K_m\right)}^2}\left\{\frac{v_m}{K_m}+\frac{v_f}{K_f}-\frac{1}{K_{dr}}\right\}$
$A_{33}$	$\frac{v_f{b}_f^{\ast }}{B_f^{\ast }{K}_f}-{\left(\frac{b_f^{\ast }}{1-K_f/K_m}\right)}^2\left\{\frac{v_m}{K_m}+\frac{v_f}{K_f}-\frac{1}{K_{dr}}\right\}$

Equations

$$\nabla \cdot \mathbf{\sigma }+\rho \mathbf{g}=\overline{\gamma },$$

$$\mathbf{ϵ}={\nabla }^{sym}\mathbf{u}=\frac{1}{2}(\nabla \mathbf{u}+{\nabla }^{\top }\mathbf{u}),$$

$$\overline{\gamma }=\sum _{\alpha =m,f}{\gamma }_{\alpha }{\mathbf{v}}_{l,\alpha },$$

$$\sum _{\alpha =m,f}{\gamma }_{\alpha }=0.$$

$${\mathbf{q}}_{\alpha }={\varphi }_{\alpha }\widetilde{{\mathbf{v}}_{l,\alpha }}={\varphi }_{\alpha }({\mathbf{v}}_{l,\alpha }-{\mathbf{v}}_s),$$

$${\varphi }_{\alpha }=\frac{{\Omega }_{p,\alpha }}{{\Omega }^0},\text{ such that }\varphi =\sum _{\alpha =m,f}{\varphi }_{\alpha },$$

$${\mathbf{q}}_{\alpha }=\frac{-{\mathbf{k}}_{\alpha }}{{\mu }_l}(\nabla p_{\alpha }-{\rho }_l^0\mathbf{g}),$$

$$\text{pderiv}{m}_{l,\alpha }t+{\rho }_l^0\nabla \cdot {\mathbf{q}}_{\alpha }={\gamma }_{\alpha },$$

$$\text{d}{W}_s=\mathbf{\sigma }\text{d}\mathbf{ϵ}+p_m\text{d}{\varphi }_m+p_f\text{d}{\varphi }_f.$$

$$\text{d}{W}_s=\text{d}{\Psi }_s,$$

$$\mathbf{\sigma }\text{d}\mathbf{ϵ}+p_m\text{d}{\varphi }_m+p_f\text{d}{\varphi }_f-\text{d}{\Psi }_s=0.$$

$$F_s={\Psi }_s-p_m{\varphi }_m-p_f{\varphi }_f,$$

$$\mathbf{\sigma }\text{d}\mathbf{ϵ}-{\varphi }_m\text{d}{p}_f-{\varphi }_f\text{d}{p}_f-\text{d}{F}_s=0.$$

$$\mathbf{\sigma }=\sigma 1+{\mathbf{\sigma }}_d,$$

$$\mathbf{ϵ}=\frac{1}{3}ϵ1+{\mathbf{ϵ}}_d,$$

$$\sigma =\frac{\partial F_s}{\partial ϵ}\text{;}{\mathbf{\sigma }}_d=\frac{\partial F_s}{\partial {\mathbf{ϵ}}_d}\text{;}{\varphi }_m=\frac{\partial F_s}{\partial p_m}\text{;}{\varphi }_f=\frac{\partial F_s}{\partial p_f}.$$

$$\text{d}\sigma =K_{dr}ϵ-b_m\text{d}{p}_m-b_f\text{d}{p}_f,$$

$$\text{d}{\varphi }_m=b_mϵ+\frac{1}{N_m}\text{d}{p}_m+\frac{1}{Q}\text{d}{p}_f,$$

$$\text{d}{\varphi }_f=b_fϵ+\frac{1}{Q}\text{d}{p}_m+\frac{1}{N_f}\text{d}{p}_f.$$

$$\text{d}{\mathbf{\sigma }}_d=2G{\mathbf{ϵ}}_d$$

$$\frac{\text{d}{m}_{l,\alpha }}{{\rho }_l}=\text{d}{\varphi }_{\alpha }+{\varphi }_{\alpha }\frac{\text{d}{\rho }_l}{{\rho }_l}\equiv \text{d}{\varphi }_{\alpha }+{\varphi }_{\alpha }\frac{\text{d}{p}_{\alpha }}{K_l},$$

$$\frac{1}{K_l}=\frac{1}{{\rho }_l}\frac{\text{d}{\rho }_l}{\text{d}p}.$$

$$\text{d}{\xi }_m=b_mϵ+\frac{1}{M_m}\text{d}{p}_m+\frac{1}{Q}\text{d}{p}_f,$$

$$\text{d}{\xi }_f=b_fϵ+\frac{1}{Q}\text{d}{p}_m+\frac{1}{M_f}\text{d}{p}_f,$$

$$\text{d}{\xi }_{\alpha }=\frac{\text{d}{m}_{l,\alpha }}{{\rho }_l^0},$$

$$\frac{1}{M_{\alpha }}=\frac{1}{N_{\alpha }}+\frac{{\varphi }_{\alpha }^0}{K_l}.$$

$$b=b_m+b_f=1-\frac{K_{dr}}{K_s},$$

$$\frac{1}{M_{\alpha }}+\frac{1}{Q}=\frac{b_{\alpha }-{\varphi }_{\alpha }^0}{K_s}+\frac{{\varphi }_{\alpha }^0}{K_l}.$$

$${\gamma }_m=\eta \frac{{\rho }_l^0\overline{k}}{{\mu }_l}(p_f-p_m),{\gamma }_f=\eta \frac{{\rho }_l^0\overline{k}}{{\mu }_l}(p_m-p_f),$$

$$\eta =\frac{2{\pi }^2}{d^2},$$

$${\sigma }^{\prime }=\text{d}\sigma +\sum _{\alpha =m,f}{\psi }_{\alpha }b\text{d}{p}_{\alpha },$$

$${\psi }_{\alpha }=\frac{{\phi }_{\alpha }}{1-{\phi }_s}.$$

$$b_{\alpha }\text{pderiv}ϵt+(1-{\varphi }_s)\text{pderiv}{\psi }_{\alpha }t+\frac{{\varphi }_{\alpha }^0}{K_l}\text{pderiv}{p}_{\alpha }t+\frac{{\mathbf{q}}_{\alpha }}{K_l}\nabla p_{\alpha }+\nabla \cdot {\mathbf{q}}_{\alpha }=\frac{1}{{\rho }_l^0}{\gamma }_{\alpha }.$$

Full text

1 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 et al. 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 phenomenon is known as poromechanical coupling (Rutqvist and Stephansson 2003), and is described by the well established poromechanical theory (see for example Biot 1941, 1977; Detournay and Cheng 1995; Coussy 1995, 2004; Wang 2000; De Boer 2012; Cheng 2016).

Aifantis (1977, 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 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 (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 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 (Fig. 1).

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

2 Double-porosity mathematical model

We present the balance equations and constitutive laws for the dual-continuum system within the macroscopic framework of Coussy (1995, 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

[TABLE]

where $\bm{\sigma}$ is the Cauchy stress tensor, $\bm{g}$ is the gravity vector, $\rho=\rho_{s}(1-\phi)+\rho_{l}\phi$ is the density of the bulk medium, $\rho_{s}$ is the intrinsic density of the solid matrix, $\rho_{l}$ is the intrinsic density of the fluid, and $\phi$ is the Lagrangian porosity. This property is defined as the ratio of the current pore volume, $\Omega_{p}$, to the bulk volume of the undeformed configuration, $\Omega^{0}$, where superscript [math] 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, $\overline{\gamma}$ represents the momentum transfer arising as a result of the mass transfer between the two pore continua. Often $\overline{\gamma}$ 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 strain-displacement compatibility relation

[TABLE]

where $\bm{u}$ denotes the displacement vector.

The inter-continuum momentum transfer is given by

[TABLE]

where $\gamma_{\alpha}$ $[\alpha=m,f]$ is the rate of mass transfer from pore continuum $\alpha$ to pore continuum $\beta$, and $\bm{v}_{l,\alpha}$ is the absolute fluid velocity within each pore continuum. The rate of mass transfer is conservative between the matrix and the fractures thus

[TABLE]

The absolute fluid velocity, $\bm{v}_{l,\alpha}$, is related to the volume flux, $\bm{q}_{\alpha}$, within each continuum by

[TABLE]

where $\tilde{\bm{v}_{l,\alpha}}$ is the relative fluid velocity, $\bm{v}_{s}$ is the velocity of the solid matrix, and $\phi_{\alpha}$ is the Lagrangian porosity associated to each continuum. This is defined as

[TABLE]

where $\Omega_{p,\alpha}$ is the current pore volume of continuum $\alpha$. The volume flux for each pore continuum is then given by Darcy’s law

[TABLE]

where $\bm{k}_{\alpha}$ and $p_{\alpha}$ denote the permeability tensor and fluid pressures associated with pore continuum $\alpha$.

The balance of fluid mass for each continuum is then given as

[TABLE]

where $m_{l,\alpha}=\rho_{l}\phi_{\alpha}$ is the fluid mass content of continuum $\alpha$.

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, 1979; Wilson and Aifantis 1982; Elsworth and Bai 1992, Lewis and Ghafouri 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 2004). It can be shown from a purely macroscopic approach (Coussy 2004) or via micromechanical considerations (Dormieux et al. 2006), that the increment in strain work density, $\text{d}W_{s}$, on the skeleton due to the loading triplet ($\text{d}\bm{\epsilon},\text{ d}p_{m},\text{ d}p_{f}$) can be expressed as

[TABLE]

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

[TABLE]

where $\Psi_{s}$ denotes the Helmholtz free energy of the skeleton, and from which it follows

[TABLE]

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

[TABLE]

into Eq. 11 results in

[TABLE]

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

[TABLE]

where $\sigma=\frac{1}{3}\text{tr}(\bm{\sigma})$ is the mean stress, $\bm{\sigma}_{d}$ is the deviatoric component of the stress tensor, $\epsilon=\text{tr}(\bm{\epsilon})$ is the volumetric strain, and $\bm{\epsilon}_{d}$ is the deviatoric component of the total strain tensor. Making use of the stress and strain decompositions from Eqs. 14 to 15 in Eq. 13, the state equations for double-porosity poroelasticity are then given as

[TABLE]

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

[TABLE]

where parameters $K_{dr}$ and $G$ are the drained bulk and shear moduli of the dual-medium respectively (Coussy 2004). Coefficients $b_{\alpha}$ can be thought of as effective Biot coefficients, and relate changes in effective Lagrangian porosity to skeletal straining under drained conditions. Coefficients $\frac{1}{N_{\alpha}}$ relate changes in the Lagrangian porosity of continuum $\alpha$ to changes in fluid pressure of the same medium, whilst the skeleton remains constrained and fluid pressure in continuum $\beta$ remains constant. Finally, $\frac{1}{Q}$ is a coupling coefficient that relates changes in the Lagrangian porosity of continuum $\alpha$ and pressure changes in continuum $\beta$.

In poromechanics it is common to formulate the constitutive equations in terms of the fluid mass content such that

[TABLE]

where $K_{l}$ is the fluid compressibility given by

[TABLE]

With Eq. 21 we can express Eqs. 18 to 19 as

[TABLE]

where

[TABLE]

Comparison between Eqs. 18 to 19 and Eqs. 23 to 24 gives the additional relation

[TABLE]

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

[TABLE]

where $b$ is the single-porosity Biot coefficient (Berryman and Wang 1995).

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 et al. (1963) the simplest model for mass transfer between the two continua is given by

[TABLE]

where $\overline{k}$ is the interface permeability, which here is taken as the matrix permeability (Barenblatt et al. 1960; Choo and Borja 2015), and $\eta$ is the shape factor. The first order nature of Eq. 29 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 Lim and Aziz (1995). For an isotropic material in two dimensions $\eta$ is defined as

[TABLE]

where $d$ denotes the average spacing between the fractures.

3 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 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:

1. (i)

Khalili and Valliappan (1996) - Constituent mechanical properties, assuming the high permeability, low storage continuum is all void space (no intrinsic fracture properties), 2. (ii)

Borja and Koliji (2009) - Constituent pore fractions, assuming the high permeability, low storage continuum is all void space, 3. (iii)

Berryman (2002) - 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 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), Choo et al. (2016) and Wang and Sun (2018). 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 Aifantis (1977, 1979), Wilson and Aifantis (1982), Elsworth and Bai (1992), Berryman and Wang (1995), Tuncay and Corapcioglu (1995, 1996), Loret and Rizzi (1999) and Dormieux et al. (2006) to name but a few. These works should thus by borne in mind in the section to follow.

3.1 Model approaches and assumptions

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}\phi^{0}_{f}=b_{f}\phi^{0}_{m}$. This relation 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 Berryman and Wang (1995), Loret and Rizzi (1999), and Dormieux et al. (2006).

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

[TABLE]

where $\sigma^{\prime}=K_{dr}\epsilon$, and $\psi_{\alpha}$ denotes the pore fraction of continuum $\alpha$ such that

[TABLE]

Notation $\varphi_{\alpha}$ is the Eulerian porosity of continuum $\alpha$. This is defined as the ratio of the current pore volume, $\Omega_{p,\alpha}$, to the bulk volume of the current (deformed) configuration, $\Omega$. Notation $\varphi_{s}=\frac{\Omega_{s}}{\Omega}$ is the current volume fraction of the solid phase. In the limit of infinitesimal transformations, $\varphi_{\alpha}\approx\phi_{\alpha}$, from which, together with small perturbations in Lagrangian porosity, follows $\psi_{\alpha}\approx\psi^{0}_{\alpha}$.

Comparison of Eq. 31 with Eq. 17 leads to the following relations for effective Biot coefficients, $b_{m}=\psi^{0}_{m}b$ and $b_{f}=\psi^{0}_{f}b$.

Borja and Koliji (2009) identify the requirement for a constitutive expression for $\psi_{\alpha}$ based on energy conjugacy with $p_{\alpha}$ (their Eq. (76)). However, explicit constitutive equations for $\psi_{\alpha}$ 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 $\psi_{\alpha}$ we first consider the following mass balance equation given in Borja and Koliji (2009),

[TABLE]

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. 8 such that

[TABLE]

where we have used Eq. 26 to decompose $\frac{1}{M_{\alpha}}$.

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

[TABLE]

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 $\pderiv{\psi_{\alpha}}{t}\approx 0$. This assumption is equivalent to

[TABLE]

in Eq. 35. We discard the latter relation in Eq. 36 as it is overly restrictive with respect to fluid exchange between matrix and fractures. Hence we identify the closure condition $\pderiv{\psi_{\alpha}}{t}\approx 0$ with values for material coefficients $\frac{1}{N_{\alpha}}$ and $\frac{1}{Q}$ of zero.

A summary of the coefficient models from Borja and Koliji (2009), under the explicit assumption $\pderiv{\psi_{\alpha}}{t}\approx 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 $\phi^{*}_{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 ($\sigma$ and $p_{\alpha}$ as primary variables; contrary to conventional poromechanical modelling),

[TABLE]

where the goal is then to find expressions for coefficients $A_{11}-A_{33}$. Symmetry of the coefficient matrix in Eq. 37 is implied. A comparison of Eq. 37 with Eq. 17 and Eqs. 23 to 24 reveals the following relations

[TABLE]

To identify expressions for the parameters in Eq. 37, Berryman (2002) considers scenarios of uniform expansion (or contraction). These scenearios are equivalent to asking whether we can find variations in uniform stress, $\text{d}\sigma=\text{d}\sigma_{m}=d\text{d}\sigma_{f}$, and variations in pore pressures, $\text{d}p_{m}$ and $\text{d}p_{f}$, such that $\epsilon=\epsilon_{m}=\epsilon_{f}$ (Berryman 2002) As a result the authors are able to relate Eq. 37 to intrinsic constituent equations. The microscale equations used in Berryman (2002) 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). In the case of a Gassmann-material, the constitutive equations for a material $\alpha$ can be written as

[TABLE]

where $v_{\alpha}$ denotes the volume fraction of material $\alpha$, and poroelastic coefficients intrinsic to a material $\alpha$ have been denoted by superscript ∗. Volume fractions for matrix and fracture materials are given as

[TABLE]

where $\Omega^{0}_{m}$ and $\Omega^{0}_{f}$ are the volumes of matrix and fracture continua at reference conditions respectively. Under the void space assumption $v_{f}=\phi^{0}_{f}$. The intrinsic Biot and Skempton coefficients, $b^{*}_{\alpha}$ and $B^{*}_{\alpha}$ for material $\alpha$ respectively, are

[TABLE]

where $K^{\alpha}_{s}$ is the solid grain modulus of the intact material, and $M^{*}_{\alpha}$ and $\phi^{*}_{\alpha}$ are the intrinsic Biot modulus and intrinsic porosity (measurement at reference conditions is implied) respectively.

Expressions for $A_{11}-A_{33}$ are finally recovered (Table 2), using the uniform expansion/contraction thought experiments described above. With the relations in Eq. 38 we get material coefficient formulations pertaining to the conventional mixed compliance setting ($\epsilon$ and $p_{\alpha}$ as primary variables) (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 dual-continuum 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) and Borja and Koliji (2009) both make an underlying void space assumption for the high permeability, low storage continuum. 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 the starting point (or the requirement for constitutive expressions are at least identified in the case of Borja and Koliji (2009)). 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.

4 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 1993). 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.

4.1 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 Bear and Bachmat (2012). 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

[TABLE]

where $\bm{f}$ is an arbitrary field and $\Omega$ 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 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, $\bm{s}(\bm{x})$ and $\bm{e}(\bm{x})$ respectively, are related through the drained fourth order stiffness and compliance tensors, $\bm{C}^{*}_{dr}(\bm{x})$ and $\bm{S}_{dr}^{*}(\bm{x})$ respectively, by

[TABLE]

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,

[TABLE]

where $\overline{\bm{e}}_{\alpha}=\Omega^{-1}_{\alpha}\int_{\Omega_{\alpha}}e_{\alpha}(x)\text{ }dV$ (resp. $\overline{\bm{s}}_{\alpha}$). Substitution of Eqs. 43 to 44 in Eqs. 45 to 46, and making use of the decompositions in Eqs. 14 to 15 due to the isotropy of the dual-medium leads to

[TABLE]

where $\overline{e}_{\alpha}$ and $\overline{s}_{\alpha}$ are the average microscopic volumetric strain and mean stress fields of continuum $\alpha$ respectively, and $S_{dr}$ is the drained compressibility of continuum $\alpha$.

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)

[TABLE]

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

Additionally, $A_{\alpha}$ (resp. $B_{\alpha}$) admit the following compatibility relations

[TABLE]

Use of Eqs. 49 to 51 in Eq. 47 and Eq. 48 allows us to develop the following fundamental relationships for effective moduli of two phase composites,

[TABLE]

Finding effective moduli then amounts to determining $A_{\alpha}$ or $B_{\alpha}$.

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}$ (Table 1). This dependence is further explored in the following section.

4.2 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 (Torquato 1991).

We use arguments from micromechanics to show heuristically that the inverse to explicit decoupling is to assume that Eq. 52 or Eq. 53 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: $\frac{1}{Q}=0$

Whilst this explicit assumption has not been used within literature its consideration remains instructive. We define assumption $\frac{1}{Q}=0$ in terms of $\text{d}\phi_{m}$ (although converse arguments may be used for $\text{d}\phi_{f}$). The term $\frac{1}{Q}=0$ is then equivalent to

[TABLE]

Next we define the following local constitutive model for a continuum $\alpha$ (Dormieux et al. 2006),

[TABLE]

Consider the local model given by Eqs. 55 to 56 for the case $\alpha=m$. Under a drained matrix ($\text{d}p_{m}=0$), Eq. 56 shows that a zero matrix porosity variation ($\text{d}\phi_{m}=0$) can only hold if the average matrix strain is zero ($\overline{e}_{m}=0$). Further, we can see from Eq. 55 that if the matrix is drained and does not deform, then the variation in average matrix stress is also zero ($\text{d}\overline{s}_{m}=0$). This will be a useful result for discussions on the explicit decoupling assumption $A_{23}=0$.

Proceeding with Eq. 54 and the macroscopic strain constraint, $\epsilon=0$, the strain partition in Eq. 46 shows that if $\epsilon=\overline{e}_{m}=0$ then $\overline{e}_{f}=0$. The condition $\frac{1}{Q}=0$ is therefore consistent with a condition of isostrain,

[TABLE]

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

[TABLE]

Hill (1963) 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 Eq. 52. 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 ($\text{d}p_{f}=0$), and from Eq. 56, effective fracture strain is equal to the effective fracture pore strain, $\overline{e}_{f}=\frac{\text{d}\phi_{f}}{\phi^{0}_{f}}$, and thus in analogy to Eq. 57 isostrain can be summarised in this case as

[TABLE]

Note that Eq. 59, and isostrain more generally, Eq. 57, are only equal to zero under the constrained macroscopic strain condition ($\epsilon=0$). However, in general Eq. 57 and Eq. 59 are non-zero due to macroscopic deformation ($\epsilon\not=0$).

From Eq. 19, with drained matrix and fracture continua, together with the void-space isotrain condition, Eq. 59, we get the following lower bound on $b_{f}$

[TABLE]

From Eq. 27 the lower bound on $b_{f}$ corresponds to an upper bound for $b_{m}$

[TABLE]

where $v_{m}=1-\phi^{0}_{f}$, and where we have used $K_{dr}=K^{V}_{dr}=v_{m}K_{m}$.

From the first equality in Eq. 54 we expect $\frac{1}{Q}\leq 0$ since matrix porosity must reduce in order to accommodate the pressure driven fracture expansion (see also similar arguments in Berryman and Wang (1995)). Thus, based on the arguments described in this section we can infer that the explicit decoupling assumption $\frac{1}{Q}=0$ is concurrent with an upper bound on $\frac{1}{Q}$.

Incompressible grain isostrain: $\frac{1}{N_{\alpha}}=\frac{1}{Q}=0$

We now consider the coefficient models from Borja and Koliji (2009) under the assumption $\pderiv{\psi_{\alpha}}{t}\approx 0$ made by Choo and Borja (2015) and Choo et al. (2016).

In Section 3.1.2, we indetified that $\pderiv{\psi_{\alpha}}{t}\approx 0$ amounts to $\frac{1}{N_{\alpha}}=\frac{1}{Q}=0$ when mapping to the constitutive model shown in Eqs. 17 to 19. It is of interest to see under what conditions the result $\frac{1}{N_{\alpha}}=\frac{1}{Q}=0$ arises when starting from the void space coefficient models built from constituent mechanical properties (i.e. those from Khalili and Valliappan (1996)).

The set of explicit assumptions: $\frac{1}{N_{\alpha}}=\frac{1}{Q}=0$, can easily be derived from the microscale by first considering isostrain (and thus $\frac{1}{Q}=0$). With the resulting bounds arising from isostrain, Eqs. 60 to 61, along with the assumption $K_{s}=\infty$, we obtain $\frac{1}{N_{\alpha}}=0$ using the coefficient models of Khalili and Valliappan (1996) (Table 1, with Eq. 26 to decompose $\frac{1}{M_{\alpha}}$). We therefore refer to conditions resulting in $\frac{1}{N_{\alpha}}=\frac{1}{Q}=0$ as incompressible grain isostrain

As far as parameters in the balance of mass are concerned, Eq. 33 and Eq. 34 are identical when assuming $\pderiv{\psi_{\alpha}}{t}\approx 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_{\alpha}$ 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 Eqs. 60 to 61.

Under the incompressible grain assumption the upper bound for $b_{m}$ now reads $b_{m}=v_{m}$. When using $K^{V}_{dr}=v_{m}K_{m}$ in the coefficient models of Khalili and Valliappan (1996) we can see that the bounds in Eqs. 60 to 61 are naturally recovered (see Table 1). In contrast, from Table 1 the effective Biot coefficients calculated using the models of Borja and Koliji (2009) 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 Nguyen and Abousleiman (2010), Kim et al. (2012), Mehrabian and Abousleiman (2014) and Mehrabian (2018). We continue to work in terms of $\text{d}\phi_{m}$. Accordingly, and using similar arguments to those used in the derivations of Eqs. 23 to 24,

[TABLE]

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. 62, $\text{d}\sigma=0$, and the stress partition in Eq. 45, if $\text{d}\sigma=\text{d}\overline{s}_{m}=0$ then $\text{d}\overline{s}_{f}=0$. The condition $A_{23}$ is therefore consistent with a condition of isostress,

[TABLE]

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

[TABLE]

Hill (1963) showed that the Reuss average is a lower bound on $K_{dr}$, and is naturally obtained by substitution of $B_{f}=1$ into Eq. 53, or from $A_{23}=0$ in Table 2. 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 Eq. 64 that $K_{dr}$ is also zero. Using the void space assumption and the isostress condition, Eq. 63, with Eq. 17 for $d\sigma$ and Eq. 55 for $d\sigma_{m}$ and $d\sigma_{f}$ gives us the following

[TABLE]

From the required isostress equality, $\sigma=\overline{s}_{f}$, Eq. 65 followed by Eq. 27 allows us to establish the following bounds

[TABLE]

Eq. 66 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. $\phi^{0}_{f}\leq b_{f}\leq 1$ are very similar to those established on $b$ for the single-porosity model (see Dormieux et al. (2006)).

4.3 On moduli upscaling

In Nguyen and Abousleiman (2010), Kim et al. (2012), Kim and Moridis (2013), Mehrabian and Abousleiman (2014, 2015) and more recently Mehrabian (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 1963). 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

[TABLE]

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

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 (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. Comprehensive summaries of such approaches can be found in the works of Aboudi (1992), Nemat-Nasser and Hori (1993) and Torquato (2002) for example.

5 Qualitative analysis using the Mandel problem

We use solutions to the double-porosity Mandel-problem (Nguyen and Abousleiman 2010), 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 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 Khalili and Valliappan (1996) and Berryman (2002) respectively. Second, we study the effects of decoupling assumptions.

5.1 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. 4). The domain is sandwiched between two impermeable, rigid plates, and is free to displace both laterally and vertically. A constant compressive force, $\int^{a}_{-a}\sigma_{zz}\text{ }dx=-2Fa$, is applied at the rigid plate boundaries, $\Gamma_{N}$ and $\Gamma_{S}$ (north and south boundaries respectively). The east and west boundaries, $\Gamma_{E}$ and $\Gamma_{W}$ respectively, are then free to drain such that $p_{m}=p_{f}=0$ at these boundaries.

In summary the boundary conditions are,

[TABLE]

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

For an isotropic double-porosity material Eq. 17 and Eq. 20 can be re-written and extended to

[TABLE]

where $\lambda=\frac{2Gv}{1-2v}$ is the Lamé constant, in which $v$ is Poisson’s ratio. Substitution of Eq. 69 and Eq. 2 in Eq. 1 yields

[TABLE]

where we assume $\overline{\gamma}\approx 0$. It can be shown that $\nabla^{2}\bm{u}=\nabla\epsilon$, with which, and in the absence of body forces, Eq. 70 reduces to

[TABLE]

where $c_{\alpha}=\frac{b_{\alpha}(1-2v)}{2G(1-v)}$ is the consolidation coefficient belonging to material $\alpha$.

Integration of Eq. 71 leads to

[TABLE]

where $f(t)$ is an integration function. Use of Eq. 72 in Eq. 34 leads to a set of diffusion equations written entirely in terms of continuum pressures. Solutions to the resulting system of equations are presented in Nguyen and Abousleiman (2010), and Mehrabian and Abousleiman (2014). Further, these solutions allow for calculation of vertical stress and strain (see Nguyen and Abousleiman 2010 or Mehrabian and Abousleiman 2014).

5.2 Data for analysis

For the qualitative analysis we use a quarter of a $2\text{ m}\times 2\text{ m}$ deformable porous domain. The studied domain is subjected to a constant top boundary force, $\int^{1}_{0}-2\times 10^{6}\text{ Pa }dx=-2\text{ 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\text{ GPa}$, $K_{s}=70\text{ 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 $\dagger$). Fluid properties are for that of water; $\rho^{0}_{l}=1000\text{ kg m}^{\text{-3}}$, $\mu_{l}=1\text{ cp}$ and $K_{l}=2.3\text{ GPa}$. Rock properties related to fluid storage and flow are $\phi^{0}_{m}=0.05$, $k_{m}=0.01\text{ md}$ and $\phi^{0}_{f}=0.01$ (Nelson 2001). Effective permeability of the fracture is assumed to be $k_{f}\approx 1000\text{ md}$. With the cubic law (Witherspoon et al. 1980), this corresponds to an aperture, $a_{f}=7\times 10^{-5}\text{ m}$, and fracture spacing, $d=2.8\times 10^{-2}\text{ m}$.

5.3 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 Khalili and Valliappan (1996) could still be used even when a fracture has an associated phase stiffness.

From Tables 1 to 2, we hypothesise that provided that $\phi^{*}_{f}\approx 1$ and the fracture phase stiffness is orders of magnitude lower than the grain stiffness ($K_{f}\ll K^{f}_{s}$), effects arising due to deviations of intrinsic fracture coefficients under the void space assumption (i.e. deviations from: $b^{*}_{f}=1$, $\frac{1}{M^{*}_{\alpha}}=\frac{1}{K_{l}}$ and $B^{*}_{\alpha}=1$) are negligible. Void space coefficient models could then be used in place of coefficient models that include intrinsic fracture properties. When $\phi^{*}_{f}=1$ we use coefficient models from Khalili and Valliappan (1996), and when $\phi^{*}_{f}<1$ we use coefficient models from Berryman (2002) 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 $\phi^{*}_{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^{\dagger}_{f}=\frac{K^{f}_{s}}{1750}$), versus ones that use a fracture modulus only an order of magnitude lower than the solid modulus, ($K^{\dagger}_{f}=\frac{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 Alberto et al. (2019) and Hajiabadi and Khoei (2019). To mimic implicit assumptions we consider $\frac{1}{Q}=0$ and $A_{23}=0$ and make no acknowledgement of these assumptions with respect to relations between mechanical properties. When considering $\frac{1}{Q}=0$ and $A_{23}=0$ we use coefficient models from Khalili and Valliappan (1996) and Berryman (2002), with $K_{f}=0$, 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}$ (a fundamental difference between the implicit assumptions considered here and explicit assumptions) we take an arbitrary value of $K^{\dagger}_{dr}=10\text{ 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 Khalili and Valliappan (1996). Under isostrain the composite and constituent bulk moduli are linked by the Voigt average ($K^{V}_{dr}=19.8\text{ GPa}$). This leads naturally to bounds on the effective constitutive coefficients, with $\frac{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\text{ GPa}$ and $K^{\dagger}_{dr}=10\text{ 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 $\pderiv{\psi_{\alpha}}{t}=0$, and Khalili and Valliappan (1996) under the assumption of incompressible grain isostrain. We use $K^{V}_{dr}=19.8\text{ GPa}$ and $K_{s}=\infty$ 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 ($A_{23}=0$) 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: $\phi^{*}_{f}=0.7$ with $K^{f}_{s}=K_{s}$ and $K^{\dagger}_{f}=\frac{K_{m}}{500}$. Coefficient models from Berryman (2002) are then used.

We compare results arising when calculating the composite bulk modulus with the Reuss average ($K^{R}_{dr}=3.3\text{ GPa}$), the HS lower bound ($K^{HS-}_{dr}=5.7\text{ GPa}$), and an arithmetic average of the HS bounds ($K^{AHS}_{dr}=12.7\text{ 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.

6 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 Eq. 70. In the instantaneous problem, fluid pressure for continuum $\alpha$ can be shown to be a state function of total stress and fluid pressure $\beta$. In the time-dependent problem, continuum pressures are governed by way of the diffusion equation, Eq. 34.

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^{t(0^{+})}_{l,\alpha}=m^{t(0)}_{l,\alpha}$. From Eq. 23 and Eq. 24, together with the undrained condition ($\text{d}m_{l,\alpha}=0$) we can recover

[TABLE]

We note that $\epsilon_{zz}\propto\epsilon$. Substitution of Eq. 17 into Eq. 73 leads to

[TABLE]

6.1 Case 1: Intrinsic fracture properties

Fig. 5 and Fig. 5 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 ($\phi^{*}_{f}\approx 1$) and the fracture phase stiffness is orders of magnitude lower than the solid grain stiffness ($K_{f}\ll K^{f}_{s}$) (Fig. 5). Even when the intrinsic fracture porosity is diminishing ($\phi^{*}_{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. 5). However, when the fracture bulk modulus is only an order of magnitude lower than the solid modulus ($K_{f}\not\ll K^{f}_{s}$), early time fracture pressure differences become measurable as intrinsic porosity decreases (Fig. 5). This phenomenon can be explained by considering the intrinsic Skempton coefficient, $B^{*}_{\alpha}$, written in the following form for the fracture phase (Cheng 2016),

[TABLE]

Use of Eq. 75 allows us to cast our observations as a bounding problem such that

[TABLE]

The lower bound in Eq. 76 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, Eq. 75 shows that as $\phi^{*}_{f}$ approaches zero, $B^{*}_{f}$ must asymptotically approach one.

If the lower bound in Eq. 76 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}\ll 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 $\phi^{*}_{f}$ decreases (Fig. 5).

Fig. 5 and Fig. 5 show the variation in vertical strain for the softer and stiffer fracture phases respectively. In Fig. 5, 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, Fig. 5 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. 73 for the fracture phase, which scales proportionally with the variation in fracture pressure.

6.2 Case 2: Implicit decoupling

Pressure and vertical strain results for the implicit decoupling assumptions test case are presented in Fig. 6. When assuming $A_{23}=0$, the matrix and fracture pressure evolutions are almost identical to the reference case (Fig. 6), with vertical strains also being correspondingly very similar (Fig. 6). When assuming $\frac{1}{Q}=0$, the early time matrix and fracture pressures are measurably lower than the reference case (Fig. 6). The early time vertical strain is greater when $\frac{1}{Q}=0$ than the strain in the reference case (Fig. 6).

The results in Fig. 6 suggest that the assumption $\frac{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 $\frac{1}{Q}<0$. From Eq. 74, setting $\frac{1}{Q}=0$ thus has the effect of removing a pressure source from continuum $\alpha$. Removal of this poromechanical coupling explains the lower than expected induced matrix and fracture pressures. From Eq. 73 we can see $\text{d}p_{\alpha}\propto\epsilon^{-1}$. Underestimated pressures therefore explain the over estimated strain when taking $\frac{1}{Q}=0$.

In contrast, Eq. 38 shows that although $A_{23}=0$, $\frac{1}{Q}\not=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.

6.3 Case 3: Explicit decoupling - isostrain

Fig. 7 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. 74.

An alternative, heuristic approach to explaining the pressure distributions in Fig. 7 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 Fig. 7 can be explained by way of Eq. 73 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 Eqs. 60 to 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, Eq. 70. This explains the growth in vertical strain separation at later times shown in Fig. 7.

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. 8 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 $\pderiv{\psi_{\alpha}}{t}\approx 0$, induced matrix and fracture pressures are equal. This can be explained by considering Eq. 73 which allows us to equate induced variations in matrix and fracture pressures such that

[TABLE]

Assuming $\pderiv{\psi_{\alpha}}{t}\approx 0$ together with the effective Biot coefficient expressions from Borja and Koliji (2009) (Table 1), Eq. 77 reduces to

[TABLE]

from which it is easy to see that $\text{d}p_{m}=\text{d}p_{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 $\text{d}p_{m}=\text{d}p_{f}$ therefore suggests that the closure condition $\pderiv{\psi_{\alpha}}{t}\approx 0$ may be an even stronger assumption than incompressible grain isostrain alone.

Fig. 8 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. 73, which is affected by differences in $b_{\alpha}$ arising from each set of coefficient models.

In light of the discussions in Section 4.2.2 and the results presented herein, assuming $\pderiv{\psi_{\alpha}}{t}\approx 0$ appears to be a strong closure assumption to make. Thus, we suggest further development of a constitutive model for $\psi_{\alpha}$, along with its relationship to the constitutive model shown in Eqs. 17 to 19.

6.4 Case 4: Explicit decoupling - isostress

Figure 9 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 (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. 9 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, Eq. 70.

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.

7 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:

1. (i)

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

Constituent pore fractions; assuming the high permeability, low storage continuum is all void space, 3. (iii)

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 $\phi^{*}_{f}<1$ and $K_{f}\not\ll 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 $\phi^{*}_{f}<1$ and $K_{f}\not\ll 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 $\phi^{*}_{f}<1$, provided $K_{f}\ll 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.