Multi-porous extension of anisotropic poroelasticity: consolidation and related coefficients

We propose the generalisation of the anisotropic poroelasticity theory. At a large scale, a medium is viewed as quasi-static, which is the original assumption of Biot. At a smaller scale, we distinguish different porosity clusters (sets of pores or fractures) that are characterized by various fluid pressures, which is the original poroelastic extension of Aifantis. In consequence, both instantaneous and time-dependent deformation lead to fluid content variations that are different in each cluster. We present the equations for such phenomena, where the anisotropic properties of both the solid matrix and pore sets are assumed. Novel poroelastic coefficients that relate solid and fluid phases in our extension are proposed, and their physical meaning is determined. To demonstrate the utility of our equations and emphasize the meaning of new coefficients, we perform numerical simulations of a triple-porosity consolidation. These simulations reveal positive pore pressure transients in the drained behaviour of weakly connected pore sets, and these may result in mechanical weakening of the material.


Introduction
The deformation of a porous medium containing fluids can be described using equations proposed by Maurice Anthony Biot; an applied physicist born in Belgium.His phenomenological theory relates strains of a solid phase with displacements of fluids.A medium is viewed at a macroscopic, bulk scale, where all pores are treated as interconnected.Constant fluid pressure throughout the medium and a unique measure of the fluid content change is adopted.
His theory of three-dimensional consolidation describes time-dependent deformation.It assumes quasi-static stress conditions, incompressible fluid, and flow obeying Darcy's law.In 1935, the isotropic version of consolidation was formulated by Biot in French.Six years later, the more rigorous treatment was written in English (Biot, 1941).
Finally, Biot discussed a more general process of anisotropic consolidation (Biot, 1955).The developments from the aforementioned papers formulate the core of the so-called (quasi-static) poroelasticity.The quasi-static approach will also be considered in this paper.
Since the publication of the original theory of poroelastic consolidation, some generalisations were proposed.Biot (1956) considered a viscoelastic extension.Further, Biot (1972) formulated the finite elastic description of porous structures.Other researchers introduced more modifications.Booker and Savvidou (1984) generalized the consolidation equations by including a thermal stress term.Chemical effects were considered by Sherwood (1993).Various impacts were combined to provide multi-physical formulations for porous materials by Taron et al. (2009).They linked the usual mechanical and hydraulic poroelastic coupling with thermal and chemical effects.Note that all the above extensions go beyond the theory of poroelasticity.These generalisations constitute foundations of the theories of finite poroelasticity, porothermoelasticity, porochemoelasticity, or porothermochemoelasticity, respectively (Chapter 12, Cheng, 2016).In this paper, however, we do not go beyond poroelasticity.Instead, we consider an extension within the frame of elastic, mechanical-hydraulic coupling.
Such an extension within Biot's poroelastic frame was first formulated by Aifantis (1977Aifantis ( , 1980a,b,c),b,c), a physicist born in Greece.His generalisation of consolidation equations was based on the idea of the so-called dual or multipleporosity.A "multi-porous" medium consists of two or more sets of pores that are characterized by different fluid pressures and diffusions.In other words, the assumption of uniform pore pressure-originally introduced by Biot-is removed.It is important to mention that the concept of spatially varying pressure has been known long before the aforementioned work.Barenblatt et al. (1960) were the first to introduce the dual-porosity approach and Warren and Root (1963) were the first to discuss its use in practice.However, these authors considered the problem of fluid diffusion only.As noticed by Elsworth and Bai (1992), until the work of Aifantis (1977), the hydrologic description was not coupled with mechanics since the stresses were assumed to be constant (Odeh, 1965;Duguid and Lee, 1977).In the classical dual-porosity analysis of fluid diffusion, the fractures and solid matrix contribute to the overall behaviour in distinctly different ways due to decidedly different characteristics; the fractures commonly have high permeability and low storativity, whereas the solid matrix tends to have the opposite characteristics.Therefore, it is reasonable to expect the poromechanical contributions of the fractures to vary from that of the solid matrix.The works by clusters that are either weakly connected or isolated from each other.Then, our multi-porous extension of anisotropic poroelasticity theory can be used to describe the mechanics of such media.Although the original Biot theory has a macroscopic nature (quasi-static assumption), in our extension, the mesoscale structures are also distinguished (different pore sets), as schematically depicted in Figure 1.Further, one can also search for the micromechanical link to the theory (pore geometry); such linkage is discussed in our parallel paper (Adamus et al., 2023).We subjectively distinguish three porosity types, where our extension is pertinent: hierarchical porosity, complex porosity, and porosity formed by spatially-distributed pore sets.Each type is presented in Figure 2 and discussed below.macro meso Figure 1: Schematic view of the multi-porous extension that considers a quasi-static medium at a macroscopic scale and distinguishes pore sets at a mesoscopic scale.Sets are allowed to be weakly connected (see upper part).
Hierarchical porosity describes a scenario, where there are multiple pore sizes in the medium (Halliwell et al., 2022) forming a nested structure, like in a Russian matryoshka doll (Cowin et al., 2009).Such porosity can be found in geomaterials, extended framework materials, foams, fibres, vascular plants, or body tissues.Examples of hierarchically porous rocks involve carbonates during the dolomitization process (Fig. 2a) or microporous sandstones having macropores at grain boundaries (Fig. 2b).The former consist of microporosity, abundant intercrystalline pores, and some irregular vuggy pores that form during the advanced stage of the dolomitization process (Wang et al., 2015).
The dolostones are interesting not only from a purely geological perspective (dolomitization) but also from a resource exploration point of view-approximately 50% of the world's gas and oil reservoirs are in carbonate rocks (Xu et al., 2020), and many ore deposits are hosted in dolostones (Warren, 1999).The porosity of clean sandstones is simpler; macroporosity forms a relatively uniform intergranular network, and microporosity forms from detrital and authigenic clays (Thomson et al., 2019).Sandstones are also important for petroleum geology since they commonly host oil and gas (e.g., Bjørlykke and Jahren, 2015).Biomechanical examples involve bone tissue, tendon tissue, or meniscus tissue; extension of the poroelasticity theory can be used to model the mechanical and blood pressure load-driven movements (Cowin et al., 2009).Also, synthesis and applications of hierarchically structured porous materials has become a rapidly evolving field of current interest (Wu et al., 2020).The theory of dual porosity was applied to, for instance, perforated concrete (Carbajo et al., 2017) and sound-absorbing materials (Bécot et al., 2008), like foams and fibres (Liu et al., 2022).
We refer to a complex porosity where microcracks or larger fractures intersect a medium having background porosity.Such a scenario commonly occurs in geomaterials (e.g., coals, tight-gas sands).Among many fields, the dual-porosity theory is widely applied in petroleum science (e.g., Lemonnier and Bourbiaux, 2010;Karimi-Fard et al., 2006;Luo et al., 2022).Originally, the dual-porosity concept was studied to describe the mechanics of conventional fractured reservoirs (Warren and Root, 1963), represented by e.g., fractured porous sandstone (Fig. 2c).Also, the extended theory of Biot was used to describe unconventional reservoirs, such as coalbed methane (Huang et al., 2022).Coals (Fig. 2d) contain natural fractures with different permeability than the porous background (Espinoza et al., 2016).
Alternatively, dual porosity can be applied to a fractured geothermal reservoir (Mahzari et al., 2021), such as hydrothermally altered granite (Genter and Traineau, 1992).Note that fractured rocks are usually anisotropic due to the preferred orientation of the fractures and, at a large scale, fracture networks can demonstrate strongly different fluid flow behaviours (e.g., Berkowitz, 2002).Hence, two or more pore sets might need to be considered to describe such phenomena.
We use the term "spatially-distributed pore sets" to describe pores that form non-overlapping concentrations.In other words, pores (e.g.cracks) having fixed poroelastic constants are not dispersed throughout the medium, but they constitute a set that is spatially separated from pores having different characteristics (see Fig. 1).For instance, Heidsiek et al. (2020) examine a reservoir sandstone with pore sets that vary spatially when viewed at both the samplescale (Fig. 3a) and the grid-cell scale relevant to reservoir studies (Fig. 3b).Brantut and Aben (2021) measured local fluid-pressure variations within laboratory-scale samples of sandstone and granite.At a large field scale, pore sets and permeabilities also vary spatially.However, to use the multi-porous extension, the considered medium cannot be too large since the solid-matrix stiffness must be uniform, and quasi-static assumptions must be obeyed.To evaluate large spatial domains, like reservoirs, a discretization into smaller scales may be needed (Karimi-Fard et al., 2006).
Having discussed each porosity type and indicated studies, where the simplified dual-porosity extension has already been applied, we can now list the anisotropic multi-porous extension potential applications.We believe that the aforementioned extension might be successfully applied to structures having: • hierarchical porosity with more than two gradations (micropores, mesopores, macropores), • complex porosity with microcracks or fractures allowing preferential flow (background porosity, microcracks, fractures with slow flow, macro-fractures with fast flow), • porosity formed by more than two spatially-distributed pore sets, • mixtures of the porosity types listed above (micropores, macropores, fractures or distributed hierarchical porosity or distributed complex porosity).

Strain-stress relations
Herein, we study the constitutive equations that describe the properties of a deformable elastic medium containing pores.In these equations, the deformation is viewed at a given instant.Strains and fluid content changes are related to the changes in stresses and pore pressures through poroelastic constants.In other words, all discussed variables (like stress, pore pressure, etc.) are viewed as instantaneous changes, rather than as absolute values.The deformation is assumed to be small so that linear relations can be utilized.In this section, the process of material consolidation or swelling is not considered since time dependency is excluded.We should also note that the sign convention here follows that of elasticity; stresses and strains are positive in the tensile direction.A list of symbols can be found in Appendix A.
Before we describe the constitutive equations, we should define the strains of the porous medium.The bulk strain tensor, where ui are the displacements of a skeleton (solid with empty pores) in the i-th direction; xi denote the coordinate axes.Throughout the paper, i, j ∈ {1, 2, 3}.In the case of isotropic medium and hydrostatic confining pressure, the displacements in the principal directions are equal.Hence, the volumetric strain can be written as e := i ∂ui/∂xi = 3eii.The change in the fluid content is more difficult to describe due to the varying nature of material porosity.Let us use a superscript (p), where p ∈ {1, . . ., n}, to denote a specific set in a n-porosity medium.Following Biot (1962), we define where is the displacement of the fluid contained in the p-th pore set; the volume fraction occupied by such a set is denoted by φ (p) ≡ V (p) /V .Expression (2) describes a relative volumetric strain of a fluid with respect to the solid, loosely referred to as "fluid increment" or "fluid content change" of a p-th pore set.Importantly, the fluid content changes between sets are not taken into account yet.In other words, ζ (p) should be viewed as a quantity that depends on the external behaviour of the bulk medium considered.
First, let us consider the constitutive equations for isotropic single porosity (Biot, 1941).They relate volumetric strain e and fluid increment ζ ≡ ζ (1) , to the changes in confining pressure pc and changes in fluid pressure p f ≡ p (1) f , through drained bulk modulus K and poroelastic coefficients.Changes in pressure are positive in compression.In where S and B are the storage and Skempton coefficients, respectively.
Let us consider the constitutive equations for isotropic dual porosity (Berryman and Wang, 1995).Due to the various characteristic of pore sets, fluid content and pressure changes are not constant throughout the medium.The increments ζ (1) , ζ (2) are related to the pore pressure changes of the first p (1) f and second pore set p The aij coefficient matrix must be symmetric because the scalar of the two remaining matrices is the potential energy (Berryman and Wang, 2000).The existence of the potential energy function implies the invariance of the coefficient matrix under permutations of subscripts i and j.The above form was first given by Berryman and Wang (1995); the original sign convention is preserved.We can treat the strains and stresses as tensors-instead of scalarsand represent them as 6 × 1 vectors.They are related by the elastic compliances S ijk .Hence, without changing any assumptions, we can rewrite isotropic relations (3) as (2) =: −a13/3 .This form is analogous to the one shown by Berryman and Wang (2000).Note that pressure changes have the opposite sign to stress changes, σij.
The equations for isotropic dual porosity can be translated into a general anisotropic case as follows.
We see that coefficients b (p) transformed into second-rank tensors, and negative signs have cancelled.Factor 2 in front of certain compliances appeared due to the index symmetries of the stress tensor.The meaning of aij and b (p) ij will be explained in the next section.Straightforwardly, we generalise the above form to the multiple-porosity, n-set 1)  . . .
where e, b (p) , σ are 6 × 1 vectors, S is a 6 × 6 matrix and T denotes the transposition.In the next sections, we will strive to get more insight into the novel equations ( 4)-( 5) given above.
4 Determination of a ij and b

(p) ij
In their present form, aij and b (p) ij are difficult to measure in practice, and their physical meaning is unclear.Therefore, in this section, we describe them in terms of elastic compliances and poroelastic coefficients (Skempton-like and storages).We consider various poroelastic boundary conditions that lead to the aij and b (p) ij determination.This way, a few universal constraints for weakly-connected sets are also provided.
The standard boundary conditions of poroelasticity are referred to as the drained and undrained states, defined by no change in fluid pressure and no change in fluid content, respectively.Another possible limit is that of constant confining stress, σ = 0 (this is not an absolute value).Such a condition is strived to be achieved in the fluid injection tests.However, in the case of a multiset extension, more variables require more boundary conditions that need to be considered.Besides, if sets are weakly connected, the long-time limit leads to pressure equilibration.In turn, the constraints for weakly-connected sets can be formulated.Below, we invoke those scenarios (test types) that provide the aforementioned constraints and useful information on aij and b (p) ij .For simplicity, we assume dual porosity (two pore sets) that can be easily extended to n > 2 considerations.

Drained test, long-time limit
In this test type, pressure throughout the pores is in equilibrium and constant.In other words, p Hence, strains from equation ( 4) are simplified to e = Sσ.In such drained conditions, the elastic compliances (S ijk ) can be determined.In the case of isotropy, we get a11 = 1/K, where K is the drained bulk modulus.

Undrained test, long-time limit
This scenario provides us with the first constraint of the theory, assuming that weak connections exist between sets.
If a rock sample has a closed system-so that the total fluid content change, ζtot, is zero-we get, and after a sufficiently long time, p In this way, we obtain the following constraints and These are valid only if the sets are not isolated (pressure equilibration over time is possible).In the case of uniaxial stress, we can derive bulk Skempton coefficients from (7), and obtain the undrained elastic compliances S u ijk from (6), Note that due to the pressure equilibration, the aforementioned bulk Skempton coefficients are equivalent to the Skempton coefficients of a single-porosity system.

Fluid injection test, long-time limit
Another constraint may be derived if we perform a standard test of fluid injection, where the applied stress is constant (σ = 0).We can again assume that the pore pressure equilibrates throughout the medium.This way, we may get the third (and last) long-time constraint that describes the storage coefficient of the bulk medium equivalent to the storage of a single-porosity system.Note that upon combining ( 8)-( 10), the undrained compliances are where ∆ ijk is defined as the difference between the undrained and drained compliance tensors; hence, it accounts for the effect of fluids.

Fluid injection test, drained first set
Consider another scenario where again, the applied stress is constant σ = 0. Let the fluid be injected directly into the first set (e.g., background porosity) so that upon a longer time, it becomes drained, p (1) f = 0, since fluid migrated towards the second set.If sets are weakly connected, it may happen that the pressure in the second set is-on average-still changing due to the fluid outflow.On the other hand, the change of fluid content is almost null at the connection points between sets.Hence, we consider a period when pore pressures in a medium are not equilibrated yet, p (2) Another experiment can lead to σ = 0, p (1) f = 0, and p (2) f = 0 if the fluid was injected in the second set, and in a short period, it did not have enough time to migrate towards another set.This way, the first set could remain dry (or drained).
In the case of either experiment, we get f .
Thus, we can define eij p (2) f σ=0, p (1) where b (2) ij is viewed as a poroelastic expansion due to fluids in the second set, scalar S (1,2) is the storage coefficients of weak connections between both sets, and S (2) is the storage coefficients of the second set.If sets are isolated, a23 = 0.It makes sense that isolated sets lead to S (1,2) = 0 since there are no connections where the fluid could be stored.In Figure 3, we illustrate the physical meaning of aij by considering sets 1 and 2 as two independent fractures with storages S (1) = a22, and S (2) = a33, respectively.Their connection would correspond to their intersection line, which holds its own storage S (1,2) = a23.13)-( 14) and ( 16) relate to physical objects.The coefficient a 23 is the cross-term describing the storage coefficient of the intersections of set 1 and set 2.
A mapping of this term to the example of two intersecting fracture sets is shown in 2D (a) and 3D (b).

Fluid injection test, drained second set
If the set numbers are interchanged, then an analogous fluid injection experiments can be performed, where σ = 0, p f = 0, and p (1) where b (1) ij is the poroelastic expansion due to fluids in the first set, and S (1) is the storage coefficient of this set (see Figure 3).
Assuming the uniaxial stress, we can define where B (1) ij is the Skempton-like tensor of the first set and S u (1) ijk is the undrained compliance tensor of the first set.

Now the meaning of b
(1) ij becomes more tangible since knowing the definition of a22 we can write b (1) Having ( 19), we notice that definition of S u (1) ijk is analogous to the meaning of S u ijk from (11), namely, k , where ∆ (1) ijk accounts for the effect of fluid caused by the first set.

Undrained second set, drained first set
In the analogous test, the order of sets is switched.This way, p f = 0 and ζ (2) = 0.In the case of uniaxial stress, we get (2 k .(1) ij }, . . ., {a

Summary of conditions
In the fluid injection tests, all sets except one must be drained (n combinations in total).The alternative tests that determine Skempton-like tensors and undrained compliances also require the sets to be drained except for one that needs to be undrained (again n combinations in total).The exact amount of test types can be easily deduced and is listed in

Alternative formulations of strain-stress relations
It might be beneficial to reformulate strain-stress relations by describing aij and b (p) ij in terms of storages and Skempton-like coefficients, as shown in ( 13)-( 14), ( 16), (19), and (21).We can rewrite equation (4) in the tensorial notation as, where in the idealized case of perfectly isolated sets, S (1,2) = 0. Having storage and Skempton-like coefficients that describe each set, we rewrite multi-porous equation ( 5) as where q ∈ N, 1 ≤ q ≤ n, and q = p.If certain set p is isolated from q, then S (p,q) = 0.The equations for all isolated sets are given and discussed explicitly in our parallel paper (Adamus et al., 2023).
So far, we have presented two versions of strain-stress relations for a multiple-porosity system.Equations ( 5) used coefficients a ij , whereas equations ( 22)-( 23) utilized more tangible Skempton-like and storage coefficients.
Below, we introduce new parameters that facilitate the alternative description of equations relating strains and stresses useful for further, time-dependent analysis.After algebraic manipulations on equations ( 22)-( 23), we obtain mixed strain-stress formulations, where again q ∈ N, 1 ≤ q ≤ n, q = p, and C ijk denotes elasticity tensor.Further, is the Biot-like tensor for each set of pores and describe fluid storage under no frame deformation.To grasp the physical meaning of the above-mentioned coefficients, let us think of fluid injection tests analogous to the ones from Section 4.4-4.5,where e = 0 instead of σ = 0 is required.
Then, we obtain p respectively.Note that definitions ( 26) and ( 27) are analogous to the Cheng (1997) definitions for single porosity (his equations ( 20) and ( 22)).The mixed formulation will be explicitly used in the consolidation equations in the next section.Also, one can try to introduce stress-strain relations by switching ζ (p) with p (p) f and reformulating poroelastic coefficients (Mehrabian, 2018).However, in contrast to the original Biot theory, such a formulation may become complicated due to the existence of multiple pore pressures.In a simpler case of isolated sets, we obtain where the undrained elasticity parameters Stress-strain relations for a more complicated case of weakly-connected sets are given in Appendix B.

Governing equations of consolidation
Herein, we study the governing equations that describe the transient phenomenon of three-dimensional consolidations.
These are differential equations governing the distributions of stress and fluid content change and settlement as a function of time in a medium under given loads.If a fluid content is increased due to an increase in the volume of pores, the governing equations describe the opposite process of swelling instead of consolidation.For convenience, we derive the case of dual porosity.Nevertheless, at the end of the section, the multi-porous generalisation is additionally presented.
To obtain the governing equations, first, changes in externally applied stresses must satisfy the equilibrium conditions; namely, We assume that body forces can be neglected.Then, we insert ( 24) into (30) and rewrite strains in terms of displacements (1).We obtain three governing equations To obtain the rest of the governing equations, let us define the flux of the fluid through pore set p as Having defined the fluid content changes (2) and the fluid flux (32), we can formulate the equations of fluid continuity.
Assuming that the fluid is incompressible and that a certain amount diffuses internally, where the interset fluid flow ∂ζ (p,q)  ∂t := Γ (p,q) p (p) is assumed to be proportional to the difference in pore pressures.Equations ( 33)-( 34) state that the external change in the amount of fluid in a pore set at any instant is balanced by the fluid flowing into the set whilst a certain amount of fluid is diffused inside the medium.It is required that ∂ζ (1,2) /∂t = −∂ζ (2,1) /∂t since the internal diffusion between the sets must be equal; fluid volume is assumed to be conserved and not compressed.Note that in the case of n > 2 sets, there is more than a single leakage constant Γ.For instance, if n = 4, we get six possibilities, namely, Γ (1,2) , Γ (1,3) , Γ (1,4) , Γ (2,3) , Γ (2,4) , and Γ (3,4) .In total, there are n p=1 (n − p) leakage constants.Now, let us invoke a crucial constitutive relation.Darcy's law for dual porosity can be written as where µ denotes viscosity.Fluid is assumed to be of the same type throughout the medium.Note that ( 35)-( 36) are the anisotropic extensions of Aifantis (1980c) expressions.We can combine the aforementioned Darcy's law ( 35)-( 36) with the continuity equations ( 33)-( 34) and fluid content changes from strain-stress relations (25).This way, we obtain the remaining set of governing equations, namely, Three differential equations ( 31) for stress distribution and two diffusion equations ( 37)-( 38) for the fluid flow, determine five unknowns, ui, p f , and p (2) f .For n pore-sets, the governing equations are rewritten as There are n + 3 equations that determine n + 3 unknowns, ui, p f , . . ., p f .

Numerical examples
Let us consider numerical examples to demonstrate the usage of consolidation equations.In our simulations, we go beyond the typical setting of single or dual porosity assumed in the past.Herein, we consider triple porosity, where a medium contains micropores, macropores, and fractures.Also, we allow different mechanical and diffusion properties of a pore set.The solutions of the partial differential equations were obtained using finite element methods incorporated inside the Matlab Reservoir Simulation Toolbox, MRST (Lie and Møyner, 2021).Within MRST, we generalized a dual-porosity module provided by Ashworth and Doster (2019b) to fit our multi-porous extension.
Although our examples have illustration purposes only, they are selected in such a way to mimic a probable geological scenario.

Model conditions and parameters
We choose the following geometrical setting and stress-strain boundary conditions (Figure 4).The considered medium is 2 m×2 m with a regular 20 × 20 mesh.Changes in fluid pressures are induced by external compressional stress of 1 MPa imposed on the top side.Displacements are vertical only.The fluid is allowed to flow through the top boundary only.Medium is considered to be undrained at a time t = 0.At the beginning of the consolidation, the material is intact and has the properties discussed in the paragraph below and schematically presented in Figure 4.This triple-porosity case can be typical for, e.g., reservoir rocks.Although there is experimental evidence of multiporosity behaviour, the exact measurement of poroelastic coefficients in geomaterials remains challenging.To the best of our knowledge, the storages and Biot-like parameters extracted in the laboratory are available for the simplified isotropic dual porosity case only (e.g., Berryman and Wang, 1995).Therefore, similarly to other researchers, we must choose certain input values subjectively (e.g., Zhang et al., 2018).Our choices are based on the values proposed by Berryman and Pride (2002) for Weber sandstone.We assume three constituents, where the constituent (c) occupies a specific volume fraction v (c) and contains one pore set (p) with initial volume fraction φ (p) < v (c) .Set 1 corresponds to macroporosity, set 2 describes microporosity, whereas set 3 refers to fractures.Material and fluid characteristics can be found in Table 2.The majority of parameters corresponding to sets 1 and 3 are taken from Berryman and Pride (2002).The values of other coefficients-including all parameters of set 2-are chosen by us; they are bold in Table 2.
Let us discuss the volume fractions and stiffnesses chosen.We assumed that the original volume occupied by unfractured porous background given by Berryman and Pride (2002) is equally divided between constituents with macro and micropores, v (1) = v (2) , where the macroporosity itself (set 1) occupies almost twice the space of the microporosity (set 2), φ (1) ≈ 2φ (2) .In general, Young's modulus diminishes with a higher concentration of pores (e.g., Pabst and Gregorová, 2014).Due to a smaller void space in set 2, we choose Es > E (2) > E (1) , where Es is Young's
Table 2: Parameters for isotropic triple-porosity simulations.Specific values of unknowns x, y, z correspond to different cases discussed.Bolded values are chosen subjectively by the authors.Other values were provided by Berryman and Pride (2002) based on laboratory measurements.
modulus of the solid phase, E (2) is the drained Young's modulus of the microporosity constituent, and E (1) denotes drained Young's modulus of the macroporosity constituent.The fractured constituent consists of a volume given by Berryman and Pride (2002), where fractures occupy most of the phase that results in a very low Young's modulus, E 3) .According to Lutz and Zimmerman (2021), Poisson's ratio changes insignificantly when pores are spherical, and νs is close to 0.2.On the other hand, if cracks are considered, Poisson's ratio should diminish.Therefore, we assumed that νs = ν (1) = ν (2) > ν (3) .To obtain the effective elasticity needed in the consolidation equations, we calculate effective Young's modulus E and effective Poisson's ratio ν by employing a lower bound of Voigt's average, namely, Now, let us discuss the values of coefficients describing fluids and poroelastic properties.We chose the typical viscosity of water at µ = 1 cP, and we selected permeabilities corresponding to the well-known situation of highly permeable fractures and background porosity having very low permeability.We assume that the time-dependent external diffusion of the fluid directly from the background is negligible that is reflected by k (1) = k (2) = 0. Also, fluid does not exchange between micro and macro porosities, which is reflected by the intersection storage a23 = 0, interflow permeability k (1,2) = 0, and leakage Γ (1,2) = 0.However, fluid exchange is allowed between micropores (set 1) and fractures (set 3) as well as between macropores (set 2) and fractures.Therefore, inter-pore sets permeabilities and leakages k (1,3) , k (2,3) , Γ (1,3) and Γ (2,3) > 0. To check the influence that the relative differences in permeability and leakage between pore sets may inflict on internal and external flux, we introduced multiplier z equal to 10, 10 3 , or 10 5 , which relate k (p,q) and Γ (p,q) as indicated in Table 2. Leakage coefficients governing the internal flow are obtained as follows.
where δ is the shape factor (Warren and Root, 1963) assumed to be equal to π 2 , and L = 0.1m is the fracture spacing viewed as the parameter that corresponds to the fracture concentration.To check the influence of the relative difference between pore sets in poroelastic coefficients, aij, i.e., storages, and b (p) , i.e., poroelastic expansion coefficients, we introduced multipliers x and y that are respectively equal to 0.8 and 0.5, 0.5 and 0.5, or 0.5 and 0.8.
Based on the work of Selvadurai and Suvorov (2020), the Biot coefficient is usually lower for a lower concentration of pores that have similar shapes.Combining ( 19) and ( 26), we can express the isotropic expansion coefficient as, Hence, b (p) is strictly related to the Biot-like coefficient; its value is expected to be lower for lower pore concentrations.
This is why the chosen values of y are below 1 to relate poroelastic expansion coefficients between macro and micropores (sets 1 and 2, respectively).In the literature, the values of fluid storage are still weakly explored, but we suspect achieving lower storage for lower pore concentrations of similar shape, which is why x is also below 1.Since the shape of the inhomogeneity can greatly affect the poroelastic coefficient (e.g., Selvadurai and Suvorov, 2020), the interplay of x and y can mimic different pore configurations between background constituents, i.e., macro and micro pore sets.

Results
In this section, we discuss the results of the triple-porosity simulations.Specifically, we focus on the time-dependent fluid pressure changes caused by the fluid outflow and gradual drainage of the medium.Before we move to the results of the model configurations discussed previously, first, let us invoke an end member case of isolated pore sets.In other words, assume that the drainage is achieved for the fractured constituent only, whereas the background porosity clusters (sets 1 and 2) remain undrained.To obtain such scenario, we set a23 = a24 = a34 = 0, k (1,2) = k (1,3) = k (2,3) = 0, and Γ (1,2) = Γ (1,3) = Γ (2,3) = 0.The rest of the parameters from Table 2 remain unchanged.The results are provided for x = 0.8, y = 0.5, and various values of z; they are illustrated in Figure 5.We notice that the solid deformation and the drainage of fractures lead to an increase in fluid pressures of the background constituents and a decrease in fracture fluid pressure.The above phenomenon can be explained as follows.In the absence of fluid volume change within the porous background, ζ (1) = ζ (2) = 0, solid compression leads to the increase of a volumetric strain, e, that must be compensated by the increase of p (1) f and p (2) f ; see e.g., equation ( 25).Further, the increased aforementioned pressures equilibrize the decrease of p (3) f caused by the fluid outflow from the fractures; see e.g., equation ( 24).Due to the isolation of the background sets, their pressures must increase till the drainage process is finished.Naturally, the complete drainage appears earlier for larger k (3) .This is why p  2. In other words, we simulate the scenario of weakly-connected pore sets.Results for all combinations of x, y, and z multipliers are provided in Figure 6.
Particularly, Figures 6a-6c, have the same multipliers as in the previous case of isolated sets so that the isolated and connected scenarios can be compared directly.
Let us discuss the results.In Figures 6a-6c, background pressures increase gradually and, upon a certain time, they start to decrease abruptly.The decrease appears earlier for larger interflow permeabilities (stronger leakage).If the interflow is small, as illustrated in Figure 6c, pressures converge with the isolated sets scenario from Figure 5c up to t ≈ 10 s, after which macro and micropore fluid pressure abruptly decrease to zero.By contrast, we notice that a relatively small contrast in fracture and interflow permeabilities leads to pressure equilibration and single-porosity response over a short period of time, as shown in Figure 6a.An interesting bump of background pressures can be noticed in Figure 6b at t ≈ 1 s.This can be explained as follows.Initially, some fluid volume contained in fractures drains directly from them, and background sets behave as isolated since their pressures are lower than that in fractures, which impedes internal leakage (the background is already fully saturated).Subsequently, once pressures in the background become higher, the fluid starts to diffuse towards the fractures so that the pressures equilibrize throughout the medium.However, such an equilibration process does not happen instantaneously due to the much lower interflow permeability that causes the bump.In the case of low contrast in permeabilities, the bump is hardly visible.Figures 6d-6f, display the case in which multipliers x = y = 0.5 and therefore represent a higher contrast in poroelastic coefficients between micro and macro-pores sets.This situation exhibits similar pressure evolution to that shown in Figures 6a-6c, but with very small difference between background (macro and micro pore sets) pressures.
This small contrast in background pressures is representative of an almost dual-porosity response.This phenomenon of nearly dual-porosity behavior i.e., micro and macro pores acting as a single set, may appear counter-intuitive considering that both storage and Biot-like coefficients are greatly diminished for the microporosity (x = y = 0.5) compared to the macroporosity.However, it can be explained as follows.The absence of the external outflow from the background, ζ (1) = ζ (2) = 0, can be rewritten using (5) along with a33 = xa22, and b (2) = yb (1) as 0 = b (1) σ + a22p (1) + a23p (2) + a24p (3) , 0 = yb (1) σ + a23p (1) + xa22p (2) + a34p (3) .b 2) .In consequence, the ratio between x and y-not just their absolute valuesinfluences substantially the relationship between micro and macro-porosity pressure evolution.The importance of poroelastic coefficients in the consolidation process is again confirmed in Figures 6g-6i.These figures represent a case of relatively low storage and high Biot-like coefficients that cause high fluid pressure in microporosity.
Our simulations show the consolidation process in multi-porous media and provide certain insights on time-dependent changes in fluid pressures.We demonstrated that the pressure evolution and pore sets interplay is controlled by the relative differences in permeability and poroelastic coefficients between pore sets.By considering different combinations of pore sets permeabilities and poroelastic coefficients, we showed specific scenarios where the triple-porosity case converges to single or dual-porosity responses.In particular, a triple porosity medium will rapidly converge to a single porosity medium if the contrast between the permeability of the most permeable pore set (fractures in this example) and the inter-pore sets permeability is small (less than two orders of magnitude difference in this example) (Figures 6a,6d,6g).We also identified characteristic features, such as pressure "bumps" (Figures 6b, 6e ,6h) appearing in the consolidation process when the contrast between fracture and inter-pore sets permeability is sufficient to cause distinctive responses between pore sets but not large enough to cause pressure equilibration in the background for a significant amount of time before draining (Figures 6c,6f,6i).Our examples should not be treated as conclusive but rather encourage future investigations on other multi-porous scenarios.

Discussion
Our multi-porous extension provides a rigorous description of an instantaneous deformation along with time-dependent fluid flow and consolidation of a solid medium containing complex porous structures.The effects of pore sets and their anisotropic responses are not neglected.To obtain such an accurate description that accounts for mesoscopic inhomogeneities (pore-sets), many coefficients are needed (see e.g., Table 2).These may be difficult to obtain in both laboratory and field measurements.Nevertheless, recent studies (Brantut and Aben, 2021) revealed that an accurate determination of local variations in fluid pressure, at the scale of the laboratory sample, is possible.Given rapid developments in experimental geophysics, one should optimistically look to the future applicability of the multi-porous extension.Further, alternative tests to the ones from Sections 4-5 are possible.First, we can consider uniform-expansion thought experiments (Berryman, 2002).Second, the poroelastic coefficients can be obtained from micromechanical derivations, as presented in our parallel paper (Adamus et al., 2023).In the planned future developments, we will perform the experiments proposed in Sections 4-5.
Our analysis of time-dependent deformation for multi-porous anisotropic poroelasticity reveals potentially significant

Summary
We have proposed the multi-porous extension of anisotropic poroelasticity-quasi-static theory originally presented by Biot (1941).In the extended theory, pores or cracks form multiple sets (porosity clusters).Such sets are either weakly connected or isolated from each other, leading to non-uniform fluid content changes and fluid pressure changes within a bulk medium.Each set may induce anisotropy, meaning that pores in a set are not necessarily randomly oriented.Also, the solid matrix-in which the sets are embedded-is allowed to be anisotropic.
In Section 2, we have indicated practical scenarios pertinent to our extension: hierarchical porosity, complex porosity, and spatially-distributed pore sets.Also, we referred to the cases where a simpler theoretical extension of isotropic dual porosity was already utilized with satisfactory results.In Section 3, we formulated the strain-stress relations that employed novel poroelastic coefficients aij and b In Section 5, we discussed alternative formulations of strain-stress relations.Alternative coefficients 1/M (p) and α (p) ij were proposed and determined.1/M (p) can be viewed as storages under no frame deformation, whereas α (p) ij are the Biot-like coefficients.Relevant types of laboratory experiments were indicated.Finally, in Section 6, we introduced time dependency to obtain the governing equations of three-dimensional consolidation (39)-(40).To do so, we proposed the extended Darcy's law and fluid continuity equations.These equations were combined with mixed strain-stress formulations to obtain novel diffusion equations.The remaining governing equations were derived from the stress equilibrium conditions.In Section 7, the usage of consolidation equations was demonstrated on novel simulations of the triple-porosity case.These numerical examples demonstrated that the pore sets pressure evolution and interplay is controlled by the relative differences in permeability and poroelastic coefficients between pore sets.
The simulations of the time-dependent drained behaviour of a multi-porous anisotropic poroelastic material show that positive pore pressure transients are generated in the weakly connected pore sets, and these could potentially be of sufficient magnitude and duration to push the material towards brittle failure.

Figure 2 :
Figure 2: Literature examples (granted permissions by corresponding authors) with porosity types pertinent to the multiporous extension.a) Dolostone with irregular vuggy pores (arrows) and smaller intercrystalline pores from Wang et al. (2015).b) Porous sandstone with pores of different sizes from Farrell et al. (2014).c) Fractured sandstone with micropores from Rizzo et al. (2018).d) Coal with cleats and micropores from Panwar et al. (2017).e) Sandstone at a sample scale and f ) field scale taken from Heidsiek et al. (2020).

Figure 3 :
Figure 3: Schematic illustration showing different pore sets and how the terms in equations (13)-(14) and (16) relate to

4. 6
Undrained first set, drained second setConsider an abrupt change of stress imposed on a jacketed sample.Assume that fluid outflows from the second set due to the insertion of a tube (p (2) f = 0), whereas the first set remains approximately closed (ζ (1) = 0).Such a short-time experiment can work only if there is very little or no flow between sets-that is also required by the theory extension.We obtain, where a long-time limit is not required.In the case of dual porosity, to determine aij and b (p) ij , we need only two types of tests.A reader can choose to either perform the fluid injection tests presented in Sections 4.4-4.5 or perform the "drained-undrained" tests from Sections 4.6-4.7.In the first possibility, aij (defined as storages) and b (p) ij (defined as poroelastic expansions) are directly measured (12)-(17).In the second possibility, a combination of measured Skempton-like tensors and undrained compliances lead to indirect determination of aij and b (p) ij .Specifically, equation (18) gives a22 and b additional equation (10) from fluid-injection long-time test gives remaining a23.In the case of multiple pore sets, we require additional scenarios to determine all {a (1) ij , b be determined directly (fluid-injection tests) or indirectly ("drained-undrained" tests).

Figure 4 :
Figure4: On the left: sketch of a geometrical setting, stress, and flow conditions.On the right: sketch of a triple-porosity scenario with a mapping of selected poroelastic coefficients.Connected macropores (set 1) are in blue, micropores (set 2) in red, and fractures (set 3) in black.Macro-and micro-porosity clusters are treated as isolated from each other, whereas the other sets are weakly connected.Ouflow coming from the background sets (macro-and micro-porosity) is neglected.

Figure 5 :
Figure 5: Semi-log plots of the average fluid pressure changes with time.Initially undrained triple-porosity medium with isolated sets.Microporosity and macroporosity remain undrained, whereas the fractured constituent is drained with time.

f
approaches zero most rapidly in Fig 5c, where k (3) = 10 5 mD, and most slowly in Fig 5a, where k (3) = 10 mD.Now, let us consider examples with the input values given in Table

Figure 6 :
Figure 6: Semi-log plots of the average fluid pressure changes with time.Initially undrained triple-porosity medium with weakly-connected sets.Each set becomes drained with time.
Section 4, we determined the physical meaning of these coefficients.aij can be viewed as storages under constant stress, whereas b (p) ij are the poroelastic expansions or products of storages and Skempton-like coefficients.Several types of tests to obtain aij and b (p) ij were proposed.
Let us discuss each scenario (test type), where different boundary conditions were assumed.First, we point out test type only (drained, long-time).Note that the amount of required test types is not affected by the number of sets embedded in the solid matrix.Boundary conditions from Sections 4.1-4.3canbe straightforwardly translated into n > 2 considerations without the necessity of introducing new scenarios.Further, let us indicate those scenarios that determine aij and b three test types that require a long-time limit.Naturally, a drained, long-time test type is indispensable to determine compliances, S ijk .Further, two scenarios from Sections 4.2 and 4.3 are important to establish the constraints (6), (7), and (10).They are necessary only if sets are not perfectly isolated.In other words, in the case of isolated sets, we need one

Table 1 :
Number of boundary test types for n-set extension.The required number determines drained compli- ances (asterix), constraints coming from long-time pressure equilibration (dot), and coefficients aij and b