Segura 2008
Segura 2008
Segura 2008
SUMMARY
In a companion Part I of this paper (Int. J. Numer. Anal. Meth. Geomech. 2008; DOI: 10.1002/nag.735),
a coupled hydro-mechanical (HM) formulation for geomaterials with discontinuities based on the finite
element method (FEM) with double-node, zero-thickness interface elements was developed and presented.
This Part II paper includes the numerical solution of basic practical problems using both the staggered and
the fully coupled approaches. A first group of simulations, based on the classical consolidation problem
with an added vertical discontinuity, is used to compare both the approaches in terms of accuracy and
convergence. The monolithic or fully coupled scheme is also used in an application example studying
the influence of a horizontal joint in the performance of a reservoir subject to fluid extraction. Results
include a comparison with other numerical solutions from the literature and a sensitivity analysis of the
mechanical parameters of the discontinuity. Some simulations are also run using both a full non-symmetric
and a simplified symmetric Jacobian matrix. On top of verifying the model developed and its capability
to reflect the conductivity changes of the interface with aperture changes, the results presented also lead
to interesting observations of the numerical performance of the methods implemented. Copyright q 2008
John Wiley & Sons, Ltd.
KEY WORDS: geomechanics; reservoir engineering; discontinuities; discrete crack; zero-thickness inter-
face; HM coupling; consolidation
∗ Correspondence to: I. Carol, ETSECCPB (School of Civil Engineering), UPC (Technical University of Catalonia),
Jordi Girona 1-3, E-08034 Barcelona, Spain.
†
E-mail: ignacio.carol@upc.edu
‡
Formerly at ETSECCPB (School of Civil Engineering), UPC (Technical University of Catalonia), Jordi Girona 1-3,
E-08034 Barcelona, Spain.
1. INTRODUCTION
The effect of discontinuities may be important in many geomechanical problems, and it is essential
in most reservoir engineering problems, especially in the analysis of naturally fractured (NFRs)
reservoirs and hydraulic fracture processes.
NFRs are known since the 19th century, although it is in the second half of the 20th century
when the effect of discontinuities on oil production has been studied more in depth due to the
discovery of the giant fields in the Middle East and Texas [1].
Areas of intense fracturing will be an objective of reservoir engineers since in most cases the
discontinuities will increase the permeability of the system. Flow behaviour of NFR is usually
described through the ‘double’ or ‘dual’ porosity idealization, in which two different porosities
are differentiated: the porosity of the rock matrix itself and the porosity of the fracture network.
The latter provides a preferential pathway for the flow of the fluid stored in the former. In other
words, most of the hydrocarbon storage is provided by the matrix porosity, whereas most of the
permeability is provided by the fracture porosity. In reservoirs where the permeability of the matrix
is also important compared with the permeability of the discontinuities, the ‘dual-permeability’
conceptualization is used, in which flow takes place along the fracture network and also through
the porous medium.
NFRs are usually stress-sensitive, and their fluid flow characteristics strongly depend on changes
in effective stresses. This is because when fluid pressure in the reservoir decreases due to fluid
production, the effective overburden load decreases (if tensile stresses are considered positive
as in Part I of the paper [2]) and causes fractures to close, with a subsequent reduction in the
permeability capabilities of the system and a significant impact on the reservoir productivity. This
major problem is more exacerbated in initially over-pressured reservoirs, where fracture closure
can be very significant during hydrocarbon recovery but less important in initially under-pressured
reservoirs because most of the closure at reservoir depth has already occurred [3]. Besides the
major impact of compaction of both discontinuities and hosting rock on fluid production and
management, it can also generate important subsidence scenarios on the ground surface that can
cause important damage to surface facilities. A well-known example of subsidence due to oil
recovery in an NFR occurs at the fractured chalk formations of the Ekofisk field in the North
Sea [4].
The existence or lack of discontinuities in a reservoir influences the techniques used for oil
production. For instance, oriented boreholes can be drilled to connect sub-parallel systems of
discontinuities that do not intersect in the site [5]. On the other hand, sites with very low perme-
abilities can improve their production capabilities if new discontinuities are created extending from
the wellbore into the reservoir using hydraulic fracture techniques [6].
Uncoupled models usually lead to inaccurate predictions, and practical problems can be
only analysed using numerical techniques that capture the strong interactions between flow
and stress/strain behaviour in both the porous medium and the discontinuity (pre-existing or
developing discontinuity) levels. The companion Part I of this paper [2] presents a finite element
method (FEM)-based formulation for the hydro-mechanical (HM) coupled problem in jointed or
susceptible-of-fracturing geomaterials using ‘zero-thickness’ interface elements and double nodes
to discretize and describe the HM behaviour of discontinuities. The resulting system of coupled
equations is solved following two different numerical approaches: staggered and fully coupled.
In Part II of this paper, a series of example cases is used to verify the model and to illustrate
the results that may be obtained in practical simulations, as well as to extract some conclusions
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
COUPLED HM ANALYSIS 2105
on the various numerical techniques employed. A first simple academic case is used to verify
that the staggered procedure proposed in Part I of the paper effectively reaches the real solution
of the problem, given by the fully coupled approach. In such a task, included in Section 2, the
consolidation of a porous block is used as a reference case for which a theoretical solution is known
[7]. An analogous procedure is followed to obtain a theoretical solution for the consolidation of a
discontinuity that deforms linear-elastically. Finally, the consolidation of a porous block crossed
by a vertical joint [8] is analysed.
In Section 3, the fully coupled numerical model is used to analyse the effect of a major joint on
the HM performance of a reservoir as proposed by Guiducci et al. [9, 10]. One of the advantages
of the coupled model proposed in Part I is that it may be used to describe both pre-existing and
developing discontinuities depending on the interface constitutive model. In this case, since the
example involves a pre-existing discontinuity, the constitutive model is the one proposed by Gens
et al. [11] for rock joint. In Segura [12], the model is also applied to fluid-pressurized developing
fractures by using an interface constitutive model based on fracture mechanics concepts [13].
Finally, Section 4 includes some concluding remarks on the overall study.
Several numerical tests based on the well-known problem of consolidation are performed to verify
the staggered and the fully coupled schemes described in Part I of this paper [2]. The traditional
1D consolidation case, for which a theoretical solution exists [7], is first solved for a porous block
and for a discontinuity separately using both the fully coupled and the staggered procedures, and
these solutions are compared with the theoretical one. Later, a vertical joint is considered in the
middle of the block domain and its effect is analysed considering its permeability constant or
aperture-dependence. In this case, the problem is not 1D anymore and a theoretical solution is not
available, so that verification of both the staggered and the fully coupled formulations is made
through comparison of one with the other. The problem configuration for the block consolidation
is taken from Ng and Small [8] and it is shown in Figures 1 and 6 (without and with the vertical
joint, respectively).
10 KPa
p=0
1m
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
2106 J. M. SEGURA AND I. CAROL
∞ 2 pop
pop = ini
sin[MY] exp[−M 2 T ] (1)
m=0 M
∞ 2
UY (%) = 1− sin[MY] exp[−M 2 T ] (2)
m=0 M
where pop is the pressure exceeding the hydrostatic values (i.e. over-pressure) and Cvcm is the
consolidation coefficient:
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
COUPLED HM ANALYSIS 2107
1 Analyt t=0.0007
Cont
0.9 Cont STAGGERED
Analyt t=0.0021
0.8 Cont
Cont STAGGERED
0.7 Analyt t=0.0035
Cont
0.6 Cont STAGGERED
Analyt t=0.0070
Cont
z 0.5
Cont STAGGERED
0.4
0.3
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
excess pore water pressure/initial excess pore water pressure
0.0
0.1
0.2
settlement/ult. settlement
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
-2.2 -1.8 -1.4 -1 -0.6 -0.2 0.2 0.6
log T
For a linear–elastic discontinuity, the aperture and the effective normal stress relate through the
normal elastic stiffness:
1
An = (6)
Kn n
The Biot modulus for the discontinuity (M j ) is considered to be sufficiently high, and if the
definition of effective stress is introduced into Equation (5) together with Equation (6), the following
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
2108 J. M. SEGURA AND I. CAROL
1 Analyt t=0.0007
disc FULL
0.9 disc STAG
Analyt t=0.0021
0.8 disc FULL
disc STAG
0.7 Analyt t=0.0035
disc FULL
0.6 disc STAG
Analyt t=0.0070
z 0.5 disc FULL
disc STAG
0.4
0.3
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
excess pore water pressure/initial excess pore water pressure
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
COUPLED HM ANALYSIS 2109
10 KPa
p=0
1m
both the porous medium and the discontinuity. From a hydraulic point of view, a high longitudinal
permeability (in comparison with the permeability of the porous medium) is considered in order
to highlight the effect of the joint. This permeability is kept constant in one analysis, while in
another one it is considered as the initial value for the aperture-dependent permeability.
Comparison with the previous results of 1D porous block consolidation allows us to identify
the effect of the discontinuity in the evolution of the system and whether consideration of the
coupling between the transmissivity and the discontinuity aperture is really important or not. All
the numerical tests are run using both the staggered and the fully coupled approaches. Since in
this case no theoretical solution exists, the fully coupled (monolithic) solution is considered as the
reference solution and the staggered procedure is verified through comparison with it. The mesh
considered is a combination of two previous meshes in which the interface elements are inserted
in the continuum finite element mesh along the central line of the block.
The initial conditions do not change with respect to the classical consolidation case and an
increase in the pore fluid pressure equal to the upper distributed load occurs over the whole domain
at the beginning of the analysis. Immediately after the load has been applied, the system starts
evacuating fluid from the upper boundary where the fluid pressure is zero. The excess fluid pressure
diminishes downward from the upper boundary and therefore the effective stress decreases (positive
tension considered) promoting the porous medium compaction as well as the discontinuity closure,
both propagating downward. If the discontinuity permeability does not vary with closure, the joint
acts as a preferential path for the depletion of the system: the fluid pressure rapidly decreases
along the discontinuity, creating a pressure gradient between the joint and the surrounding porous
medium that facilitates the depletion of the porous block. The consolidation process is not 1D and
the porous block deforms more rapidly in the central (jointed) zone.
A coupled permeability (i.e. aperture-dependent) strongly influences the evolution of the system.
In this case, the fluid pressure initially diminishes faster along the joint than in the surrounding
porous medium too, but as the fluid pressure decreases also the permeability does due to the
closure of the joint. This closure appears first in the upper zone and propagates downward, so
that at early stages it occurs that the joint has consolidated (i.e. closed) on the upper part and it
cannot evacuate fluid anymore. However, it is still a preferential path for flow in the lower part
of the domain. Compared with the 1D porous medium consolidation, this situation translates into
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
2110 J. M. SEGURA AND I. CAROL
an increase of the fluid pressure at the upper zone and in a decrease of the fluid pressure in the
lower part during the first stages. The initially fast consolidation slows down after the upper part
of the discontinuity closes, consolidating at a pace similar to the porous-medium-only case. All
these considerations are observed in Figures 7 and 8, which, respectively, depict the fluid pressure
profiles along the central line of the domain (along the joint) and the settlements at the central
1 1
0.9 0.9
0.8 0.8
0.7 0.7
0.6 0.6
z 0.5 z 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
(a) excess pore water pressure/initial excess pore water pressure (b) excess pore water pressure/initial excess pore water pressure
1 1
0.9 0.9
0.8 0.8
0.7 0.7
0.6 0.6
z 0.5 z 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
(c) excess pore water pressure/initial excess pore water pressure (d) excess pore water pressure/initial excess pore water pressure
Figure 7. Fluid pressure distributions at: (a) 0.0007 d; (b) 0.0021 d; (c) 0.0035 d; and (d) 0.007 d.
0.0
without joints FULL
0.1 without joints S TAGGERED
with joints Kctt FULL
0.2 with joints Kctt S TAGGERED
with joints K (u) FULL
settlement/ult. settlement
0.4
0.5
0.6
0.7
0.8
0.9
1.0
-2.2 -1.8 -1.4 -1 -0.6 -0.2 0.2 0.6
log T
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
COUPLED HM ANALYSIS 2111
zone for the three cases analysed (no joint, joint with constant permeability, and joint with variable
permeability), and according to the two solution procedures (staggered and fully coupled).
The staggered and the fully coupled approaches give very similar solutions as shown in
Figures 7 and 8. The number of iterations in the staggered strategy is higher when the joint perme-
ability depends on the joint aperture and lower if no joints are considered. This is consistent with
the idea that the number of iterations depends on the intensity of the coupling between the pres-
sure and the displacements [14] since in the parallel plate model the permeability depends on the
cube of the aperture [15]. As the intensity of coupling increases, lower values of the stabilization
parameter for the staggered procedure, which was introduced in Part I of the paper, are needed
to avoid oscillations in the iterative solution procedure. Nevertheless, the accuracy of the method
and its convergence to the fully coupled solution are demonstrated.
Figures 9 and 10, respectively, show, for different values of time, the deformed mesh (discon-
tinuity closure shown as overlap of the continuum elements on both sides) and the fluid pressure
Figure 9. Mesh deformation (magnification factor ×50): (a) porous block alone; (b) uncoupled
transmissivity; and (c) coupled transmissivity.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
2112 J. M. SEGURA AND I. CAROL
Figure 10. Fluid pressure distributions: (a) porous block alone; (b) uncoupled
transmissivity; and (c) coupled transmissivity.
distribution for the consolidation of the porous medium block without discontinuity (a), with a
constant transmissivity discontinuity (b), and with a discontinuity that has an initial permeability
as (b) but dependent on the aperture (c). It is observed that the joint acts as an attractive path for
the depletion and that it consolidates much faster in (b). Case (c) situates somewhere in between
cases (a) and (b), and once the joint has closed in its upper part, it does not act as a preferential
path for the depletion of the system anymore.
Discontinuities usually act as preferential paths for the flow of fluid and the exploitation of reser-
voirs. However, as fluid pressure is decreased due to fluid withdrawal, the effective stresses decrease
(i.e. the medium is compressed), and therefore a progressive closure of discontinuities takes place,
especially in the pumping region. In stress-sensitive reservoirs, this closure dramatically decreases
the discontinuities’ transmissivity, and the production capabilities of the reservoir diminish.
This section analyses the effect of a pre-existing horizontal discontinuity on the depletion of
a single-phase reservoir as proposed by Guiducci et al. [9, 10]. The same problem is solved
considering different assumptions for the HM behaviour of the discontinuity, which makes it more
or less stress-sensitive to changes in fluid pressure. In particular, and as it will be shown, the
results are dramatically different if the longitudinal transmissivity is considered constant from the
beginning of the analysis, predicting much higher production rates and also higher settlements.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
COUPLED HM ANALYSIS 2113
Comparison is made between the results presented in [9, 10] and those obtained using the fully
coupled formulation presented in Part I of the paper.
U
n = −K ni Umc (13)
Umc −U
where K ni is the initial normal stiffness and Umc is the maximum discontinuity closure (Figure 12).
The discontinuity mechanical constitutive relationship used in [9, 10] is an empirical hyperbolic
Figure 11. Geometry, initial and boundary conditions for the fractured reservoir simulation.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
2114 J. M. SEGURA AND I. CAROL
Figure 12. Hyperbolic relationship between closure and effective normal stress.
Figure 13. Mesh and detail of the refinement at the pumping and the discontinuity zones.
function [16] with an analogous expression to [11]. The proposed values for the parameters are
Umc = 0.000217 m and K ni = 100 000 MPa/m.
The mesh shown in Figure 13 is used in all the analyses. It is composed of 2180 nodes, 2154
finite elements for the continuum, and 162 interface elements for the discontinuity. Strong fluid
pressure gradients are expected in the pumping area as well as around the discontinuity, and
therefore a more refined mesh is considered in these areas. With regard to the time discretization,
a backward scheme is used with time steps 0.1 years long during the first 9 years and 1 year long
until the analysis is finished.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
COUPLED HM ANALYSIS 2115
50
48
46
44
42
p (MPa)
40
38
coupled 7.5 years
36
uncoupled 7.5 years
34
coupled 20 years
32 uncoupled 20 years
30
0 500 1000 1500 2000 2500
x (m)
Figure 14. Fluid pressure profiles along the discontinuity at 7.5 and 20 years, and for coupled and
uncoupled transmissivity. Comparison to results from [9, 10] (dashed lines).
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
2116 J. M. SEGURA AND I. CAROL
q (m^2/s)
1.2E-05
1.0E-05
8.0E-06
6.0E-06
4.0E-06
2.0E-06
0.0E+00
0 500 1000 1500 2000 2500
x (m)
Figure 15. Fluid flux profiles along the discontinuity at 7.5 and 20 years, and for coupled and uncoupled
transmissivity. Comparison to results from [9, 10] (dashed lines).
uncoupled discontinuity
1.6 coupled discontinuity
uncoupled lower layer
1.4 uncoupled upper layer
coupled lower layer
1.2 coupled upper layer
Q (m^3/day)
0.8
0.6
0.4
0.2
0
0 2000 4000 6000 8000
t (days)
Figure 16. Evolution with time of the flow rate leaving the discontinuity and the porous blocks at the
pumping zone for coupled and uncoupled transmissivity. Comparison to results from [9, 10] (dashed lines).
shown in Figure 16, which could point out that the differences stated are due to differences in the
discontinuity zone refinement.
The distribution of fluid pressure within the reservoir is shown in detail for the coupled case at
7.5 and 20 years in Figure 17. It is clearly observed how the fluid pressure first decreases along
the discontinuity, creating a pressure gradient with the surrounding porous blocks that facilitates
its depletion through the discontinuity.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
COUPLED HM ANALYSIS 2117
Figure 17. Fluid pressure distributions at: (a) 7.5 years and (b) 20 years.
45
K = 1·10 Mpa/m
K = 1·10 Mpa/m
40 K = 1·10 Mpa/m
35
K = 1·10 Mpa/m
30
- σ’n (Mpa)
25
20
15
10
0
0.00E+00 5.00E-05 1.00E-04 1.50E-04 2.00E-04
U (m)
Figure 18. Initial and final pairs (n ,U ) for different values of K ni .
the previous analysis, that under the initial isostatic stress conditions (x = y = −13.3 MPa), the
initial opening of the discontinuity is coincident with the asymptotic closure. The cubic law [15]
for the discontinuity flow performance is considered in all the cases.
3.3.1. Influence of K ni . Figure 18 shows the expected path that the hyperbolic function will follow
at the pumping zone for different values of K ni and a fixed value of the asymptotic closure
Umc = 0.217 mm.
Initially, the discontinuity is highly permeable in all the cases. As the fluid pressure is decreased,
the stiffer the discontinuity initially is, the less it closes. In fact, as shown in Figures 19–21,
the charts obtained with K ni = 1×106 MPa/m are very similar to those obtained in the previous
section using a constant transmissivity, which makes sense, since with a high value of K ni the
discontinuity does not almost close during reservoir depletion, and its transmissivity alters very
little. On the other hand, a value of K ni = 1×103 MPa/m almost entirely closes the discontinuity
at the well zone during pore pressure drawdown and the discontinuity does not influence the flow
behaviour of the system anymore. The cases with K ni = 1×104 , 1×105 , and 1×106 MPa/m are
analysed and compared in Figures 19 and 20.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
2118 J. M. SEGURA AND I. CAROL
50 50
48 48
46 46
44 44
42 42
p (MPa)
p (MPa)
40 K_n = 1E+04 MPa/m 40
K_n = 1E+04 MPa/m
38 K_n = 1E+05 MPa/m 38
K_n = 1E+05 MPa/m
36 36
K_n = 1E+06 MPa/m
34 34 K_n = 1E+06 MPa/m
32 32
30 30
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
(a) x (m) (b) x (m)
Figure 19. Fluid pressure profiles along the discontinuity at: (a) 7.5 years and (b) 20 years.
1.4 discontinuity
lower layer
1.2
upper layer
1
Kn = 1e+06 MPa/m
Q(m^3/day)
0.8
0.4
Kn = 1e+04 MPa/m
0.2
0
0 1000 2000 3000 4000 5000 6000 7000 8000
t (days)
Figure 20. Evolution with time of the flow rate leaving the discontinuity and the
porous blocks at the pumping zone.
U = 0.300 mm
U = 0.217 mm
53 U = 0.134 mm
48
43
38
-σ’n
33
28
23
18
13
0.00E+00 5.00E-05 1.00E-04 1.50E-04 2.00E-04 2.50E-04 3.00E-04
U (m)
Figure 21. Initial and final pairs (n ,U ) for different values of Umc .
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
COUPLED HM ANALYSIS 2119
Figure 19 shows how the fluid pressure profiles along the discontinuity are higher for a lower
initial stiffness because the discontinuity rapidly closes and it depletes the system at a slower pace,
giving similar fluid pressure profiles at 7.5 and 20 years. On the contrary, a higher initial stiffness
implies no important variation of the discontinuity opening and, therefore, a faster depletion of
the system that gives lower pressure profiles.
Figure 20 shows that the flow rate discharged by the discontinuity at the beginning of the
analysis is the same in the three cases analysed, but it diminishes as the reservoir is depleted and
the discontinuity closes. This decrease is more dramatic for lower values of K ni because of the
faster discontinuity closure in these cases. In fact, the maximum flow rate is attained at t = 7.5 years
in all the curves but for K ni = 1×104 MPa/m. This is due to the rapid closure of the discontinuity,
which in this case makes the transmissivity to decrease at a rate greater than the imposed increase
in pressure gradient.
The case in which K ni = 1×104 MPa/m would correspond to a ‘stress-sensitive’ reservoir, since
the hydraulic characteristics of the system extremely depend on the stress level, whereas for
K ni = 1×106 MPa/m and at the stress levels studied in this case, the flow characteristics of the
discontinuity, and hence of the reservoir, suffer almost no alteration.
3.3.2. Influence of Umc . In this case, a fixed initial discontinuity stiffness of K ni = 1×105 MPa/m
is used, whereas the asymptotic closure is varied. In particular, three different cases are studied:
a higher (Umc = 0.3 mm), a mid (Umc = 0.217 mm), and a lower (Umc = 0.134 mm) asymptotic
closure. In all the cases, the asymptotic closure is considered to be also the initial discontinuity
aperture, so that contrary to the previous case, the initial aperture is different for the three cases.
The expected path that the normal stress acting on the discontinuity and the discontinuity closure
will follow is shown in Figure 21 for the different values of Umc .
Note that if the same initial aperture were considered for the three cases (in this case, the largest
value is 0.3 mm), then the results would be completely different, since the attained closure would
be lower for lower values of Umc . In this case, if the asymptotic closure did not match with the
initial aperture, a residual aperture would be reached for high compressive stresses, which would
be the difference between the initial aperture and the asymptotic closure.
The three cases are highly conductive at the beginning of the analysis in comparison to the
continuum medium, although of course the higher the initial aperture is, the more permeable the
discontinuity will be. Moving to n = −28.3 MPa, the mouth of the discontinuity (i.e. at the well)
will reach low aperture values for Umc = 0.134 mm, still high aperture values for Umc = 0.300 mm,
and the discontinuity with Umc = 0.217 mm will locate in between these two behaviours as deducted
from the numerical results plotted in Figures 22 and 23.
An important difference from the analysis in Section 3.3.1 exists, since now the aperture is
significantly reduced in the three situations, even in the more permeable case. This makes the three
situations to exhibit a similar degree of non-linearity, which will affect the number of iterations till
convergence as discussed in the next section. This is not the case with the analyses in Section 3.3.1,
in which the higher the K ni is, the more linear the normal stress–closure relationship is, and the
lower the aperture variation is.
Figure 22 shows the fluid pressure profiles along the discontinuity mid-plane. As already
described, the lower the Umc is, the more impervious the discontinuity is during the whole analysis,
and therefore the higher the fluid pressure profiles are at a given value of time and the lower the
fluid rate discharged by the discontinuity as shown in Figure 23. It is important to note that in
this case the flow rate discharged by the discontinuity in each case is different from the beginning
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
2120 J. M. SEGURA AND I. CAROL
50 50
48 48
46 46
44 44
42 42
p (MPa)
p (MPa)
40 u_mc = 0.134 mm 40
u_mc = 0.134 mm
38 u_mc = 0.217 mm 38
u_mc = 0.217 mm
36 u_mc = 0.300 mm 36
34 34 u_mc = 0.300 mm
32 32
30 30
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
(a) x (m) (b) x (m)
Figure 22. Fluid pressure profiles along the discontinuity at: (a) 7.5 years and (b) 20 years.
1.6 discontinuity
lower layer
1.4
upper layer
1.2
1
Q(m^3/day)
umc = 0.3 mm
0.8
0.6
umc = 0.217 mm
0.4
umc = 0.134 mm
0.2
0
0 1000 2000 3000 4000 5000 6000 7000 8000
t (days)
Figure 23. Evolution with time of the flow rate leaving the discontinuity and
the porous blocks at the pumping zone.
of the analysis, since a different initial aperture is considered for each study. Therefore, the three
curves present a different slope already from the beginning. This is not the case in Figure 21,
where the initial aperture is the same for the three values of K ni , and therefore the flow rates are
very similar during the first time-steps.
3.3.3. Symmetric vs non-symmetric Jacobian matrix. The same problems solved for different
values of K ni in Section 3.3.1 and different values of Umc in Section 3.3.2 have been solved
using two different Jacobian matrices, both described in Part I of the paper: a full non-symmetric
Jacobian matrix and a symmetric Jacobian matrix. The number of iterations till convergence using
the full Jacobian matrix should be lower than using the symmetric matrix. In contrast, computation
of the full Jacobian matrix is more time consuming due to the number of operations involved in
the calculations as well as due to the higher memory requirements. However, in these analyses,
since the number of degrees of freedom is not very large, the difference in computational time per
iteration between both methods is not important.
Figure 24 shows the number of iterations till convergence for the simulations with a fixed asymp-
totic aperture Umc = 0.217 mm and different K ni , and Figure 25 shows the number of iterations for
the analyses using a fixed K ni = 2×105 MPa/m and different Umc .
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
COUPLED HM ANALYSIS 2121
25 25
20 20
sym Jacobian
iterations
iterations
10 sym Jacobian 10
full non-sym Jacobian
5 5
0 0
0 5 10 15 20 0 5 10 15 20
(a) t (years) (b) t (years)
25
20
iterations
15 sym Jacobian
full non-sym Jacobian
10
0
0 5 10 15 20
(c) t (years)
Figure 24. Number of iterations until convergence at each time-step for: (a) K ni = 1×104 MPa/m;
(b) K ni = 1×105 MPa/m; and (c) K ni = 1×106 MPa/m.
12 12
10 10
8 8
iterations
iterations
6 6
4 4
sym Jacobian sym Jacobian
2 2
full non-sym Jacobian full non-sym Jacobian
0 0
0 5 10 15 20 0 5 10 15 20
(a) t (years) (b) t (years)
12
10
8
iterations
4
sym Jacobian
2
full non-sym Jacobian
0
0 5 10 15 20
(c) t (years)
Figure 25. Number of iterations until convergence at each time-step for: (a) Umc = 0.134 mm;
(b) Umc = 0.217 mm; and (c) Umc = 0.3 mm.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
2122 J. M. SEGURA AND I. CAROL
If the full non-symmetric Jacobian matrix is used in the calculations, a similar number of
iterations is needed independent of K ni and Umc . However, if the symmetric matrix is considered, the
number of iterations depends on the problem under consideration. Comparison of Figures 24 and 25
shows how for a given time-step, the number of iterations drastically increases as the initial
stiffness is reduced according to the values considered, whereas this increase is not so pronounced
in the case of decreasing Umc . Looking at Figure 21, the three curves show a similar non-
linear shape, which is not the case in Figure 18. Therefore, the mechanical non-linear behaviour
substantially affects the number of iterations till convergence when using a symmetric Jacobian.
A given variation in normal stress gives a more important change in aperture (and, therefore, in
longitudinal transmissivity) for the case with a lowest K ni , and this translates into a higher number
of iterations for this case. As the depletion takes place, the number of iterations also increases.
The number of iterations at later time-steps (from t = 9 years on) increases because the time-steps
considered are longer. In Figure 24(a) the use of the full Jacobian matrix would be recommended
due to the great difference in the number of iterations, whereas in Figure 24(c) the symmetric
Jacobian matrix can reach the solution faster due to a lower computing time spent on each
iteration. With respect to Figure 25, all the cases could be solved with a non-symmetric Jacobian
matrix, although the difference in the number of iterations is not very important, specially in
case (c).
4. CONCLUDING REMARKS
The results reported in this Part II paper validate the formulation described in Part I [2] and illustrate
some of its capabilities. Additionally, they clearly show the great influence of discontinuities in
some coupled HM problems. The consolidation of a block crossed by a vertical joint in the centre
shows how the presence of a discontinuity accelerates the consolidation process. However, the rate
at which the consolidation of the system takes place depends on the closure of the discontinuity
due to fluid pressure decrease, which reduces the flow capabilities of the discontinuity and slows
down the consolidation process compared with an uncoupled transmissivity calculation. The second
example representing the behaviour of a producing reservoir crossed by a major horizontal joint
also shows a major influence of the mechanical and coupling properties of the interface (cubic law)
on the overall reservoir performance. From the numerical viewpoint, some interesting conclusions
may also be drawn from this study, among them (i) if converged, the staggered approach leads
to the correct solution, although strong coupling such as caused by cubic law in discontinuities
usually requires the monolithic approach for convergence, and (ii) in the monolithic approach,
a full non-symmetric Jacobian matrix is especially effective when the discontinuity behaviour is
markedly non-linear.
ACKNOWLEDGEMENTS
The first author wishes to acknowledge the FI doctoral fellowship from AGAUR-DURSI (Generalitat de
Catalunya, Barcelona) and support from ETSECCPB and UPC for completing his doctoral thesis. Partial
support was also obtained from grant MAT2005-24522 and research project BIA2006-12717, funded by
MEC (Madrid), as well as grant 80015/A04 funded by MFOM (Madrid). The early stages of this research
were also supported by Schlumberger.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag
COUPLED HM ANALYSIS 2123
REFERENCES
1. Kazemi H, Gilman JR. Multiphase flow in fractured petroleum reservoirs. In Flow and Contaminant Transport
in Fractured Rock, Bear J, Tsang CF, de Marsily G (eds). Academic Press Inc.: New York, 1992; 267–323.
2. Segura JM, Carol I. Coupled HM analysis using zero-thickness interface elements with double nodes. Part I:
theoretical model. International Journal for Numerical and Analytical Methods in Geomechanics 2008; DOI:
10.1002/nag.735.
3. Aguilera R. Geologic and engineering aspects of naturally fractured reservoirs. CSEG Recorder 2003; 28:44–49.
4. Visser J. Stabilized finite element methods for coupled geomechanics and multiphase flow. Ph.D. Thesis,
Department of Petroleum Engineering, Stanford University, Palo Alto, CA, U.S.A., 2002.
5. National Research Council. Rock Fractures and Fluid Flow: Contemporary Understanding and Applications.
National Academy Press: Washington, DC, U.S.A., 1996.
6. Valkó P, Economides MJ. Hydraulic Fracture Mechanics. Wiley: New York, 1995.
7. Taylor DW. Fundamentals of Soil Mechanics. Wiley: New York, 1948.
8. Ng KLA, Small JC. Behavior of joints and interfaces subjected to water pressure. Computers and Geotechnics
1997; 20(1):71–93.
9. Guiducci C, Pellegrino A, Radu JP, Collin F, Charlier R. Numerical modeling of hydro-mechanical fracture
behavior. In Numerical Models in Geomechanics-NUMOG VIII, Pande GN, Pietruszczak S (eds). Swets &
Zeitlinger: Lisse, 2002; 293–299.
10. Guiducci C. Numerical modeling of hydro-mechanical fracture behavior. Ph.D. Thesis, Scuola Normale Superiore
de Pisa, Pisa, Italy, 2003.
11. Gens A, Carol I, Alonso EE. A constitutive model for rock joints: formulation and numerical implementation.
Computers and Geotechnics 1990; 9:3–20.
12. Segura JM. Coupled HM analysis using zero-thickness interface elements with double nodes. Ph.D. Thesis,
Technical University of Catalonia (UPC), Barcelona, Spain, 2007.
13. Carol I, Prat P, López CM. A normal/shear cracking model. Application to discrete crack analysis. Journal of
Engineering Mechanics (ASCE) 1997; 123(8):765–773.
14. Lewis RW, Schrefler BA. The Finite Element Method in the Static and Dynamic Deformation and Consolidation
of Porous Media. Wiley: Chichester, 1998.
15. Snow D. A parallel plate model of fractured permeable media. Ph.D. Dissertation, University of California,
Berkeley, 1965.
16. Bart M. Contribution a la modlisation du comportement hydromcanique des massifs rocheux avec fracture. Ph.D.
Thesis, Univ. des Sciences et Technologies de Lille, Lille, France, 2000.
Copyright q 2008 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2008; 32:2103–2123
DOI: 10.1002/nag