Numerical manifold method modeling of coupled processes in fractured geological media at multiple scales

Mengsu Hu, Jonny Rutqvist

Lawrence Berkeley National Laboratory, Berkeley, CA, 94720, USA

ABSTRACT The greatest challenges of rigorously modeling coupled hydro-mechanical (HM) processes in fractured geological media at different scales are associated with computational geometry. These challenges include dynamic shearing and opening of intersecting fractures at discrete fracture scales as a result of coupled processes, and contact alteration along rough fracture surfaces that triggers structural and physical changes of fractures at micro-asperity scale. In this paper, these challenges are tackled by developing a comprehensive modeling approach for coupled processes in fractured geological media based on numerical manifold method (NMM) at multiple scales. Based on their distinct geometric features, fractures are categorized into three different scales: dominant fracture, discrete fracture, and discontinuum asperity scales.Here the scale is relative,that of the fracture relative to that of the research interest or domain. Different geometric representations of fractures at different scales are used, and different governing equations and constitutive relationships are applied. For dominant fractures, a finite thickness zone model is developed to treat a fracture as a porous nonlinear domain. Nonlinear fracture mechanical behavior is accurately modeled with an implicit approach based on strain energy. For discrete fractures, a zero-dimensional model was developed for analyzing fluid flow and mechanics in fractures that are geometrically treated as boundaries of the rock matrix. With the zero-dimensional model, these fractures can be modeled with arbitrary orientations and intersections. They can be fluid conduits or seals,and can be open,bonded or sliding.For the discontinuum asperity scale,the geometry of rough fracture surfaces is explicitly represented and contacts involving dynamic alteration of contacts among asperities are rigorously calculated. Using this approach, fracture alteration caused by deformation,re-arrangement and sliding of rough surfaces can be captured.Our comprehensive model is able to handle the computational challenges with accurate representation of intersections and shearing of fractures at the discrete fracture scale and rigorously treats contacts along rough fracture surfaces at the discontinuum asperity scale. With future development of three-dimensional (3D) geometric representation of discrete fracture networks in porous rock and contacts among multi-body systems,this model is promising as a basis of 3D fully coupled analysis of fractures at multiple scales, for advancing understanding and optimizing energy recovery and storage in fractured geological media.

Keywords:Dominant fractures Discrete fractures Discontinuum asperity scale Coupled processes Numerical manifold method (NMM)

1. Introduction

Fractures play key roles in subsurface energy recovery and storage,including hydrocarbon and geothermal energy production,and nuclear waste disposal. Fractures, with sizes ranging from microns to kilometers,may act as conduits or seals for fluid flow in these various subsurface energy activities (Rutqvist and Stephansson, 2003). In subsurface energy recovery (e.g. hydrocarbon recovery and geothermal exploitation), effectively creating a fracture network or making use of natural fractures is the key for efficient production. In contrast, in subsurface energy storage systems (such as carbon sequestration and nuclear waste disposal),fractures or faults may act as unfavorable flow conduits that may compromise the seal integrity of the storage facility. In geological systems involving fluid-filled fractures and porous media,complex coupled hydro-mechanical (HM) and thermo-hydro-mechanical(THM) processes can occur, including alteration of natural fractures and creation of new fractures that may intersect with natural fractures. Because fractures play key roles in energy recovery and storage and because they can be created and altered as a result of coupled processes, understanding coupled processes in fractured geological media is essential.

At reservoir scale, fractures are often not isolated entities, but form networks of interacting fractures. These fractures are usually very thin (e.g. microns to millimeters) relative to their length(meters). These fractures may dynamically alter in terms of their dimensions and physical properties, and new fractures may be created as a result of coupled processes. Since 1980s, a number of numerical models have been developed for modeling coupled HM processes in fractured rock,including equivalent continuum,dualcontinuum, and discontinuous models. As the arbitrarily oriented fractures may influence the deformation and fluid distribution in a complex way that cannot be simplified as a continuum, numerical modeling considering discrete fractures with full coupling capability is of great importance. Noorishad et al. (1982, 1992) developed a finite element model for fully coupled HM and THM processes in deformable fractured rock masses. After that,increasing demands for engineering solutions have inspired development of many computer codes capable of modeling HM and THM behaviors of fractured rock at various levels of sophistication,including those applied for the analysis of nuclear waste disposal and geothermal energy recovery(Rutqvist et al.,2001).Most of the aforementioned models were developed based on finite element method (FEM). With the development of discontinuous methods,fractures could be explicitly represented as a displacement discontinuity,as they are modeled as interfaces between contacting individual blocks. This includes models based on the distinct element method (DEM), including the commercially available universal distinct element code (UDEC) (Itasca, 2011) and 3DEC(Itasca, 2013), and models based on discontinuous deformation analysis (DDA) that treat coupled fluid flow and deformation in discrete fractures, but with rock blocks being assumed impermeable (Kim et al., 1999; Jing et al., 2001). Later, models based on enriched FEM were developed (e.g. Silvestre et al., 2015) that included simplified jump terms to capture the mechanical displacement discontinuity and hydraulic pressure continuity associated with fractures.In this method,however,the mechanical effects on hydraulic properties were not considered.

Models where fractures are explicitly simulated can be categorized depending on the geometric representation of the fractures for fluid flow and mechanics. For fluid flow in fractures, there are three types of models:n-dimensional (i.e. for two-dimensional(2D) models, fractures are represented by 2D elements),n-1 dimensional (i.e. for 2D models, fractures are represented by onedimensional (1D) elements), and zero-dimensional models recently developed by Hu et al. (2016, 2017b). As dimensions are associated primarily with degrees of freedom(DOFs)to be solved in the numerical modeling, the zero-dimensional model for discrete fractures means that no additional DOFs are required for the fractures because they are treated as discontinuous boundaries of rock matrix. For mechanics of fractures, there are two types of models:n-dimensional solid element models and discontinuous models which treat fractures as discontinuities between contacting rock blocks.The major disadvantage ofn-dimensional fluid flow models is the inaccuracy of representation of fracture apertures and the changes they are subjected to, in spite of the high computational efforts associated with the number of elements required to approximate the real geometry of intersecting fractures. The disadvantage ofn-1 dimensional fluid flow models is that they neglect fracture thickness and this results in additional required DOFs that make them incompatible with mechanical models for fractures(e.g.difficult to model initiation,shearing or re-opening of fractures). Though discrete fracture network (DFN) models have been widely developed, accurate and efficient modeling of interaction between fractures and rock matrix is lacking, especially when geomechanics plays an important role. Mechanically,ndimensional solid element models (Rutqvist et al., 2009) are excellent for representing fractures with certain apertures, but when fracture apertures are far less than their lengths,the use ofndimensional solid elements becomes too computationally expensive. Discontinuous mechanical models, on the other hand, are promising for representing thin fractures where fracture surfaces could be open, bonded or slip/shearing, but very few numerical models rigorously consider all these mechanisms.Because of these limitations, most models are not suitable for analyzing coupled processes at discrete fracture scales without rigorous consideration of intersections and shearing of the discrete fractures, particularly those that may be evolving dynamically as a result of coupled processes.

For single fractures, several models have been developed recently for analyzing micro-scale coupled processes in fractures as a part of the international DECOVALEX project (Bond et al., 2016).The models were categorized into 2D simplified models,statistical models, and homogenized (in which a rough fracture surface is treated as a single entity) models. For example, McDermott et al.(2015) presented a model that combined numerical (for flow) and analytical(for chemo-mechanical)methods to simulate small-scale coupled chemo-hydro-mechanical processes, but the model is limited to parallel edge to edge surface contacts. Pan et al. (2016)developed an elasto-plastic cellular automaton model for a single fracture using grids with different apertures to represent contacts and voids. In their model, flow and transport are calculated using conventional numerical schemes. However, the mechanics and contact alteration of the fracture are not considered in their analysis.In these models,geometric features(such as rough boundaries along fractures, pores and grains) are either not represented explicitly,or they are approximated by spheres or rectangular grids.Thus,contacts along rough surfaces cannot be accurately captured.Because of these limitations, numerical modeling of coupled processes has to the author’s knowledge never been attempted at the microscopic scale where discontinuities are important.

In this paper, numerical approaches are presented using numerical manifold method (NMM) that is able to overcome the limitations as discussed above, in order to simulate coupled processes at multiple scales.In Section 2,fractures are categorized into three different scales based on their geometric features:dominant fracture, discrete fracture, and discontinuum asperity scales. The analyses of the fundamentally different physical constitutive behaviors and coupled processes of fractures at each of these three scales are presented. In Section 3, NMM is introduced that allows simulating fractures as continuities as well as discontinuities. In Section 4,three different models for these three scales of fractures are presented. These are (1) a finite thickness zone model for dominant porous fractures with nonlinear mechanical behavior;(2)a zero-dimensional fracture model for fluid flow and mechanics in fractures networks, where fractures are boundaries of the rock matrix, and can be open, bonded, or sliding while being fluid conduits or seals; and (3) explicit model for micro-scale fractures with rigorous treatment of contacts between fracture asperities and their dynamic alteration. In Section 5, differences among these different models are summarized, and the challenges associated with numerical modeling of coupled processes in fracture at multiple scales are concluded.

2. Fracture across scales: Geometric and physical features

2.1. Three scales of fractures based on geometric features

Depending on multiple scales,fractures appear to have different geometric and physical features. Fig.1 shows discrete fractures at reservoir(100 m-10 km),outcrop(1 m—100 m),and core(1 mm—1 cm, Ajo-Franklin et al., 2018) scales. Herein, it can be seen that at the reservoir and outcrop scales, fractures appear in groups, arbitrarily oriented and intersecting with each other.At the core scale,a dominant fracture appears,possibly filled with mineral fillings and connected with smaller fractures at the surrounding rock.If looking at the microscale (1 μm—1 cm) or nanoscale, a single fracture appears to be quite rough and its tensile failure could be controlled by failure of carbon nanotubes(Taloni et al.,2018).Therefore,the three types of fractures can be geometrically summarized: dominant fractures with a certain width, discrete thin fractures, and discontinuum fractures with asperities forming rough surfaces.

Note that the scale is relative, meaning that the scales of fractures are relative to the scale of interest of the problem or domain.For example, at 1 m-1 km scale, fractures/faults could also appear with similar geometric features to that at core scale, and those different sized fractures could be considered as dominant features of their respective scale of interest.Another example is that fracture networks could exist at micro-scale,such as micro-fractures in core samples,micro-crack populations associated with subcritical crack growth. Regardless of different formation mechanisms(Voigtl?nder et al., 2018), micro-fractures in rock cores appear to have similar geometric features to meter-scale discrete fractures at reservoir.Based on this ‘relative scale’concept,it is also possible to label micro-crack populations as discrete fracture scale, which makes it reasonable to use similar types of numerical models to simulate discrete fracture networks at reservoir.

Corresponding to geometric features, physical features of fractures at different scales can also differ. For a dominant fracture, it can be porous,connected with small cracks,or filled with minerals;across its width, compression of the fracture can be more difficult for a given increment of stress, resulting potentially in a fracture with a residual aperture. For discrete fractures, they may serve as fluid conduits or seals depending on their connectivity and flux exchange with the rock matrix. They may be mechanically open,bonded, or shearing due to fluid pressure or stress. If looking at a fracture at micro-scale, on the two sides of the rough surfaces, a significant difference of the open fracture channel and rock matrix can be identified.Fluid flow mostly occurs within the open channel,while the mechanical deformation of the rock matrix impacts the geometry of the open channel.Furthermore,the deformation of the fracture surfaces can be significant, due to the lack of constraints even though the surface-to-surface contacts exist.

2.2. General governing equations for coupled processes in porous media

Regardless of the scales considered, in a fluid-saturated porous system (e.g. rock matrix or a dominant fracture), coupled HM processes satisfy conservation of momentum and mass, described by Biot’s general theory of three-dimensional (3D) consolidation(Biot,1941):

where σ is the total stress tensor,fis the body force vector,vis the fluid velocity vector,α is the Biot—Willis coefficient(usually ranges between 0 and 1),εvis the volumetric strain of the porous media,Mis the Biot’s modulus,pis the fluid pressure,andtis time.The Biot—Willis coefficient as a factor multiplied by fluid pressure in Eq. (1)signifies a modification and generalization of Terzaghi’s effective stress law to

where σ′is the effective stress tensor,mT=[1,1,1,0,0,0]for three dimensions ormT= [1,1, 0] for two dimensions.

For mechanical analysis of linearly elastic porous media with small-deformation, we have

whereEis the elastic constitutive tensor,ε is the strain tensor,Ais the strain-displacement matrix, anduis the displacement vector.

For fluid flow in porous media,it is assumed that the fluid flow satisfies Darcy’s law:

Fig.1. Fractures at reservoir, outcrop, and core (Ajo-Franklin et al., 2018) scales.

whereKis the tensor of permeability coefficient, andhis the hydraulic head or the piezometric head (i.e. the sum of the pressure head and the elevation head).

In the case of fractures, in contrast, different governing equations and constitutive relationships are applied for fluid flow and mechanics.These will be presented in detail in Section 4.In general,two types of couplings between fluid flow and mechanics exist:direct and indirect couplings(see Fig.2,as a concept presented by Rutqvist and Stephansson(2003)).Direct couplings refer to as solid deformation perturbing conservation of mass that impacts pore fluid pressure (i.e. pore-volume coupling), while fluid pressure impacts the effective stress.Direct couplings are included in Biot’s general theory(Eqs.(1)—(3)).When volumetric strain is very small,however, direct coupling will be reduced to one way. Indirect coupling means that hydraulic or mechanical properties change with deformation and/or fluid pressure, where the fracture permeability is very sensitive to deformation, and this requires accurate calculation of this indirect coupling.

3. Fundamentals of NMM

NMM (Shi, 1991) is based on the concept of “manifold” in topology. In NMM, independent meshes for interpolation and integration are defined separately. Based on this approach, an initially one-time generated, non-conforming mesh (not necessarily conforming to the physical boundaries) can be used and flexible local approximations can be constructed and averaged to establish global approximations for both continuous and discontinuous analyses.

In NMM, independent mathematical and physical covers are defined. A mathematical cover is a set of connected patches covering the entire material domain. For example, a quadrilateral/circular/rectangular patch can be used as a mathematical cover(see A,B,C in Fig.3,respectively).Features such as density and shape of these mathematical patches define the precision of the interpolation. The physical patches are mathematical patches divided by boundaries and discontinuities,determining the integration fields.The union of all the physical patches forms a physical cover. For example, physical patch C is the entire model domain, while physical patch B is divided from mathematical patch B by boundaries. Physical patch A (divided from mathematical patch A by boundaries) is further divided into physical patches A1and A2by the inner discontinuity.The overlapping areas by multiple physical patches are defined as elements.As a result,the model domain Ω is discretized into five elements: A1BC (the overlap of physical patches A1,B and C),A1C,A2C,BC,and C.From Fig.3,it can be seen that the shape of the mathematical patches can be arbitrary; the relative location of the mathematical patches to the model domain can also be arbitrary (only if satisfying Ω?A∪B∪C), and the number of physical patches on each element can be arbitrary.

Fig.2. Hydro-mechanical direct(I)and indirect(II)couplings(Modified after Rutqvist and Stephansson, 2003).

Fig. 3. NMM mathematical and physical meshes.

On each physical patch, a local function is assigned, such as constant one, linear one, or anyone that is able to capture the solution behavior on the patch. The weighted average of the local patch functions forms the global approximation. For example, if using linear local functions, a global second-order approximation could be constructed (Fig. 4a, Wang et al., 2016); if using a local function with a jump of the first derivative, a material interface crossing patches and elements could be simulated(Fig.4b,Hu et al.,2015b). Or most commonly, if using discontinuous local functions,fractures can be simulated(Fig.4c,Hu et al.,2016,2017b).With this dual-mesh concept,NMM is capable of simulating both continuum and discontinuum problems with accurate geometric representation and flexible numerical approximation.With the advantages of dual-mesh concept, NMM also has been successfully applied to analyzing moving interface problems, such as free surface flow(Wang et al., 2014, 2016; Zheng et al., 2015; Yang et al., 2019).

NMM with the concept of global approximation can be related to other numerical methods, as shown in Fig. 5. If using a bilinear weight function on rectangles with a constant local function,NMM can be simplified to the FEM. If using a piecewise linear weight function in each direction with a constant local function, NMM is then simplified to the finite volume method(FVM,Hu and Rutqvist,2020). If using a constant weight function with a constant local function,NMM is simplified to the DEM.If using a constant weight function with a linear patch function (resulting in a notable difference from DEM with constant patch function), NMM is simplified to DDA. Nevertheless, a comparison of various numerical methods on all aspects is not attempted to be made herein,including (1) interpolation/approximation, (2) construction of global equilibrium(transforming differential to integral equations),(3)approaches of integration,and(4)solving of linear or nonlinear global equations. Only the first aspect is compared, i.e. interpolation/approximation, as it defines fundamentals of a numerical method.NMM provides a flexible and general approach to include continuous and discontinuous methods in a unified form.

In this study, constant patch functions and linear weight functions composed of shape functions of triangular mathematical meshes are used to approximate temperature, fluid pressure and displacements, which are generally expressed as follows:

where φ,w and φpcare the field variables(such as hydraulic head or fluid pressure and displacements),weight function,and physical patch functions,respectively.

Fig. 4. Flexible choice of local approximation functions: (a) linear function, (b) a jump junction for a weak discontinuity, (c) a discontinuous function for a fracture.

Fig. 5. Relating NMM to other numerical methods.

4. Modeling coupled processes in fractures at multiple scales:numerical models and applications

4.1. A finite-thickness nonlinear poroelastic model for dominant fractures

In fractured rock masses,it is common to see individual dominant fractures or faults with asperities and mineral fillings that cannot be easily simplified to parallel planes. The main flow feature may be a complex geologic feature,consisting of multiple branching fractures intermingled with mineral-filled sections and damaged host rocks adjacent to fracture surfaces(Fig.6a).Another related key property is the nonlinear relationship between stress and fracture aperture.Moreover, the flow feature is also associated with a mechanical weakness that may allow for inelastic shear slip along its plane.

Fig. 6. Schematic of the simplified porous fractured rock model.

In a study focusingonthe largerscale(larger thanasinglefracture/fault),a porous finite thickness zone was developed to deal with this type of system(Hu et al.,2017a).This finite thickness zone is porous and has strong nonlinear properties reflecting fracture flow and fracture opening and/or shear behavior,in consideration of fracture filling effects. The thickness of this equivalent porous media flow feature in the model may far exceed the real fracture width including open fracture parts and filling.It can include part of the host rocks on each side of the flow feature, still retaining the key features of potential fracture flow and nonlinear deformation behavior.The model for such a feature is depicted in Fig. 6b. It is a porous medium of thicknessIdthat includes a dominant fracture flow path as well as other materials such as fracture filling and part of the host rocks.

With the concept of the porous finite thickness zone that has nonlinear flow and mechanical features, this section describes its mathematical formulation, numerical implementation, and shows an application with such a model.

Fig.7. Nonlinear relationship between σ′n and bm (Modified from Rutqvist et al.,2000).

4.1.1. Mathematical statement

In order to account for the nonlinear feature of the porous fracture,a reformulation of Bandis et al.(1983)’s equation(Rutqvist et al.,1998,2000) is used to describe the nonlinear relationship of the effective fracture normal stresswith the mechanical aperturebm(as shown in Fig. 7):

The following relationship describes the behavior of fracture shear displacement under shear stress:

where ζ and ψ are the constants;is the shear stress; and Δusis the shear displacement. Examining Eq. (11), it is concluded that when ψ = 0, the linear relation between shear stress and shear displacement is retained.

For fluid flow in fractures, the hydraulic conductivitykfof a fracture is related to a hydraulic fracture aperturebh,which can be defined according to Witherspoon et al. (1980):

where ρfand μfare the fluid density and dynamic viscosity,respectively;gis the gravitational acceleration; andbhis the hydraulic aperture assumed to be

wherebhris the residual hydraulic aperture when the fracture is mechanically closed, andfis a factor that compensates for the deviation of flow in a natural rough fracture from the ideal parallel smooth fracture surfaces.

With the above concepts and equations, the aperture of the dominant fracture flow path is used to calculate the hydraulic conductivity as given in Eq. (12). The deformation behavior of the finite thickness zone is contributed to two parts: the nonlinear behavior of the fracture described in Eq. (9), and linear elastic deformation of the solid fracture fillings and adjacent host rocks.The HM coupling includes direct pore-volume coupling, as well as indirect coupling with changes of mechanical and hydraulic properties induced by flow and deformation, respectively.

4.1.2. Numerical approach for the nonlinear constitutive behavior

As shown in Fig.6,dominant fractures are represented as a finite thickness zone with a given width and constitutive behavior differing from rock matrix.Geometrically,the fracture zone has the same dimension as the rock matrix (i.e.n-dimensional representation). Therefore, the geometric discretization for dominant fractures is the same as that for the rock matrix. Here an approach is described to model the nonlinear physical constitutive behavior of this zone.

As the deformation of the finite thickness zones consists of linear deformation of the mineral fillings and adjacent host rocks and nonlinear deformation of the fracture, the normal strain is expressed as the summation of the two:

where the superscriptiirepresents theii-th time step, and η represents the compliance of fillings and adjacent host rock within the fracture zone. Combined with Eq. (9), Eq. (14) becomes

In order to accurately account for the nonlinear behavior of the finite thickness fracture zone,an implicit approach was developed by Hu et al. (2017a). In this approach, instead of using an approximation with linear segments to the nonlinear relationship, the nonlinear mechanical behavior is fully incorporated using strain energy for the material under deformation. Therefore, those nonlinear relationships are directly introduced to strain energy as described in the following subsections for normal and shear deformations.

The normal constitutive model expressed in Eq. (15) can be rewritten as

where

The strain energy in the porous medium representing a fracture zone Πefnis expressed as

Combined with Eq. (16), Eq. (18) becomes

where

Transforming coordinates from globalx-yto locals-ncoordinate system yields

The shear constitutive model expressed by Eq. (11) can be further expressed as

where

Similar to the approach for fracture normal mechanical behavior,the associated strain energy using Eq.(11)can be written as

With Taylor expansion and coordinate system transformation,we have

where

Fluid flow in the fracture zone is governed by Eq.(2),where the volumetric strain perturbs conservation of mass (direct coupling).In addition, the hydraulic conductivity of the fracture zone is expressed by Eqs.(12)and(13),where hydraulic aperture should be updated dynamically as a result of the deformation of the fracture zone (indirect coupling).

As the fracture zones are modeled in this approach as porous media with nonlinear properties that differ from the surrounding rock,the boundaries of the fracture zones are regarded as material interfaces. The displacement continuity across these material interfaces is realized by the penalty method (Shi, 1991) combined with the continuity of hydraulic head, and the normal flux is fulfilled by the Lagrange multiplier method developed by Hu et al.(2015a).

With the energy-work based theorem(Wang et al.,2014),work and energy associated with fluid flow and mechanical deformation are implicitly constructed and updated at each time step. Minimization of the total energy at each time step leads to the global equilibrium equations.

4.1.3. Example: Coupled processes at a single dominant fracture

In order to demonstrate the formulation in consideration of both direct and indirect coupled HM processes in rock with fractures, a rectangular rock domain is simulated, which contains a fracture zone subjected to instantaneous vertical load and a constant pressure fluid injection.The simulated domain is 1 m wide and 2 m high, with a 1 m × 0.1 m horizontal fracture zone at the vertical center (Fig. 8). The mechanical fracture aperture for the assumed dominant fracture is 1 ×10-4m (0.1 mm) with an equivalent hydraulic aperture of 5 × 10-5m (50 μm). For the rock matrix, the Young’s modulus is 4 GPa, the Poisson’s ratio is 0.2, and the permeability coefficient is 5 × 10-9m/s. For the fracture zone, Bandis’parameteris 5 MPa,the shear constants ζ is 10-11Pa-1and ψ is 0, and the factorfis 0.5. Initially, the total vertical stress is -8 MPa(compressive stress)and fluid pressure is 0.

First, an instantaneous 10 MPa vertical downward load is applied on the top of the domain. A mechanical analysis without fluid injection was conducted. This results in an instantaneous closure of the fracture as a result of its nonlinear normal closure behavior with changing normal stiffness.The simulation indicates a mechanical fracture aperture of 6 × 10-5m (60 μm) at the final steady state, which is accurate according to Eq. (9). Then a simulation was conducted considering only indirect coupling, i.e. the fluid—solid interaction terms for direct coupling is deactivated. In this case, the coupling occurs only in one way, i.e. mechanical deformation affects permeability but there is no influence of fluid pressure on the mechanical field. The changes of mechanical and hydraulic properties of the fracture under load and injection with constant pressure of 8 MPa at the left end of the fracture zone and the pressure at the right end of the fracture zone are fixed at zero.Finally,the full package is run,considering both direct and indirect couplings with results shown in Fig. 8.

The distribution of fluid pressure in the following three cases is compared: (1) without considering coupling, (2) only considering indirect coupling, and (3) considering both direct and indirect couplings in Fig. 8. The difference of fluid pressure distribution in Fig.8a and b is not obvious,indicating that a steady state is reached for the case of indirect coupling only after 30-d injection.However,in Fig.8c,a steady state has not been reached and fluid continues to dissipate from the left to right.This difference can be explained by the fact that in Fig.8b(with only indirect coupling),a steady state is reached when mechanical deformation no longer occurs, whereas in the case for Fig. 8c, the final steady state will be reached only after the interaction between the mechanical and fluid flow fields is balanced.The slower convergence to steady state shown in Fig.8c is caused primarily by pore-volume direct coupling in the relatively low permeability rock matrix. As observed in the simulation, fluid flow reaches steady state much earlier in the fracture zone (see Fig. 8c). Overall, the effect of pressure on solid deformation is not obvious,which is also reflected by the rapid convergence of fracture aperture to its final values for both cases (2) and (3). Further, the aperture change is compared with time at the injection point under these two conditions. The aperture at the final stage reduces to 6 × 10-5m(60 μm)when only considering indirect coupling.This value is the same as the one in the case of pure mechanical analysis,proving its verification.However,when considering both direct and indirect couplings, the aperture remains steady at 6.5 × 10-5m(65 μm) under the effect of fluid pressure on the solid skeleton.Based on this simple example,it is known that pore-volume direct coupling may play a significant role for dominant fractures, thereforen-dimensional models which consider the direct coupling and nonlinear behavior of the fracture zone are necessary for analyzing coupled processes for this type of fractures.

Fig.8. Distribution of fluid pressure(MPa)for(a)flow analysis without considering coupled effects,30 d after injection;(b)only considering indirect coupling;and(c)considering both direct and indirect couplings. The dimensions of the simulated domain are in meter.

4.2. Zero-dimensional model for discrete fractures

At discrete fracture scales where fractures may be oriented and intersecting arbitrarily,n-dimensional andn-1 dimensional models have limitations.A zero-dimensional model was developed to simulate fluid flow, geomechanics, and their couplings. In this section, the general mathematical formulation, geometric representation, and numerical implementation are summarized and an example is given to demonstrate its efficiency and capability.

Fig. 9. Fluid flow (a) along fracture and (b) normal to fracture directions in zerodimensional model.

4.2.1. General mathematical formulation

The basic idea of a zero-dimensional fracture model is to use the surrounding elements of a fracture to represent the fracture,avoiding the need for additional dimensions or DOFs for the fractures themselves. Since the fractures are treated as discontinuous boundaries of the rock matrix, no additional DOFs are required to represent them.For fluid flow and mechanics in porous media,Eqs.(1)—(7) are still applicable. As for fractures, flux occurs both along fracture and normal to fracture directions.Fig.9 shows how this is accomplished numerically for the two flow directions (where φ represents variables, i.e. hydraulic head for flow and temperature for heat transfer).

Fig.9a shows fluid flow along fractures.Because the fractures in this case are very thin relative to their length (e.g. microns to millimeters aperture for meter-sized fractures), it can be assumed that hydraulic head within the fracture is uniform across its thickness:

where φ and φ′denote the field variables on different elements that are divided by the fracture, φf(shuō)(s) is field variable within a fracture,andsrepresents local coordinate along the fracture. Thereby,fluid along a very thin fracture is represented by flow along its two surfaces:

where αfis the permeability coefficient for Darcy flow,andis the flux.Here parallel plate flow in fractures is assumed as given in Eq.(12).

Fig. 9b shows the fracture-matrix interaction, i.e. fluid or heat exchange normal to a fracture.A fracture that is bounded by its two surfaces from the surrounding rock is explicitly represented by these two surfaces, which belongs to different elements. As the fractures are assumed to be thin and unfilled,the distribution of the field variable on each surface of a fracture and within a fracture is uniform in the direction normal to the fracture surfaces.Therefore,these two fracture surfaces are considered to consist of two Dirichlet boundaries:

Fig.10. Mechanical states of each segment of a fracture in networks: open, bonded and sliding.

The Eqs.(30)—(32)include all the possibilities of fluid flow in a fracture, which may act as a fluid conduit or seal. In contrast, the mechanical state of a fracture is more complicated.A fracture may have several segments and each of the segments could be open,bonded,or sliding,as shown in Fig.10.The contact state is impacted by fluid flow and mechanical deformation,and may be dynamically changing.In view of the mechanical states,these three states have different boundary constraints.

When the segment of a fracture is open, it is considered that it has a linear constitutive behavior:

When a segment of a fracture is bonded,the distances between the two sides of the segment should be zero:

wheredfis the vector of distance, including the normal distancedf，nbetween the two surfaces of the fracture segment,and relative distancedf，salong the two surfaces.dfis calculated by accounting for original distances, and the relative displacements of the two surfaces while they are moving and deforming.

When a segment of a fracture is sliding, the Coulomb’s law of friction is satisfied in the tangential direction, while the normal distance between the two surfaces of the fracture segment should be zero:

4.2.2. Implementation of zero-dimensional model in NMM

The key issues for modeling coupled processes in discrete fractures are how to handle (1) the geometrical representation of intersecting fractures, (2) fracture flow and fracture-matrix flow interaction, (3) deformation and dynamic contacts involving slip and opening of the fractures, and (4) couplings between flow and mechanics.

(1) Geometric representation

In order to explicitly simulate fluid flow and mechanics of discrete fractures as well as their interactions (flux exchange with rock matrix), both fractures and rock matrix need to be geometrically represented. In NMM, a triangular mathematical mesh is used to overlap the entire simulation domain. Once the mathematical mesh is generated, the same tree-cutting algorithm could be used(Shi,1991;Hu et al.,2017b)to generate the physical covers and elements, divided by fractures and boundaries,considering locations of triangle edges. Fig. 11 shows different meshes for fractured porous rock domains with 20, 60 and 150 fractures. The triangular mathematical meshes are independent from the fracture geometry. Indeed, fractures and boundaries divide mathematical meshes into arbitrarily shaped physical covers and elements.By using two coincided lines to represent a fracture,all fractures are accounted for in the calculation,including isolated fractures that do not intersect with other fractures where their interaction with the rock matrix cannot be ignored. It is also interesting that with the tree-cutting algorithm,cases of fluid flow in porous media with sealed fractures are readily modeled.

Corresponding to the mechanical states of fractures described by Eqs. (33)—(35), each fracture is discretized into several line segments and these segments may have different contact states.Thus,it is possible to consider complicated behaviors such as shear dilation or uneven opening of a fracture. The density of fracture discretization is consistent with the global meshing density (see Fig.11), which can be selected flexibly.

An important issue for calculating discrete fractures is how to simulate intersections of fractures. Fig. 12 demonstrates a geometric representation of two fractures that intersect with each other as well as with one triangular mathematical mesh.As can be seen, the two intersecting fractures divide the triangle into four different parts (A, B, C, D). Then fluid flow and contact states(satisfying constraints described by Eqs.(33)—(35))will be applied on the four pairs of parallel interfaces(interfaces between A and B,C and D,B and C,and A and D)to account for the opening,bonded and sliding states of the surfaces of each fracture.

Fractures in open and bonded states or alteration between these two states are easier to simulate because this does not require changes of contact pairs.Zero-dimensional fracture model assumes that at the initial stage,a fracture is approximated by two surfaces parallel with each other at the beginning, but these two surfaces can be non-parallel after deformation and motion,or opening of the fracture caused by fluid pressure. This capability is included in the algorithm.

A significant challenge for modeling coupled processes in fracture network is to simulate shearing along fractures,as this leads to dynamic changes of contacts between different elements.As shown in Fig.13,when the four blocks A, B,C and D are in contact (when the fractures are completely bonded),the contact pairs are A and B,B and C,C and D,and A and D.But when sliding(slip)occurs at one of the fractures,the contact pairs become A and B,B and C,B and D,C and D, and D and A. By using a rigorous contact algorithm that updates contact pairs at each time step,sliding along fractures can be rigorously and explicitly represented.

Fig.11. NMM mesh generated for fractured porous rock with (a) 20, (b) 60, and (c) 150 fractures.

Fig. 12. Geometric representation of an element with one intersection by two fractures.

The concept of a zero-dimensional fracture model was proposed by Hu et al.(2016) for the first time and it is a practical method to simulate discrete fracture networks.The zero-dimensional fracture model is distinct from zero-thickness models, as the zerodimensional model considers the thickness of a fracture and its change under deformation(i.e.the fractures have a width).With no need to introduce additional DOFs for fractures,such an approach is very flexible to model existing and initiated fractures as well as shearing and re-opening of fractures.

(2) Zero-dimensional model for coupled fluid flow and mechanics in discrete fractures

The explicit geometric representation of discrete fractures with no need to introduce additional DOFs lays the basis for the zerodimensional model for fluid flow and mechanics of discrete fractures embedded in a porous rock matrix. By using a zerodimensional fluid model, both along-fracture flow and fracturematrix interaction can be handled. In this model, effect of a residual hydraulic aperture on fluid conduction is considered even when the fractures are mechanically closed. When the mechanical aperture of a fracture is changing due to continuous opening or the change of contact state, its hydraulic aperture and hydraulic conductivity are updated. Moreover, as a result of deformation, fluid flow and their coupling,the contact state may change dynamically,and the corresponding constitutive behavior should be described differently.Such a stringent scheme that guarantees the accuracy of the results when considering complex fractured rock and HM coupling is important. The scheme also accounts for indirect couplings. Thus, coupled HM responses are able to be modeled in a complex intersecting fracture network.

An energy-work based formulation for flow is used to establish the global equilibrium equations by minimizing the total potential energy associated with each component of work. The work associated with along-fracture flowis represented as

Combining with Eqs. (7) and (36) yields

The work associated with fracture-matrix interaction flowis

where dsdenotes integration on each discretized segment of a fracture; andkf，nis the permeability coefficient in the direction normal to the fracture.

For calculation of the contact between two fracture surfaces,there are different terms of potential energy associated with different contact states described in Eqs. (33)—(35).

When a fracture segment is open,the linear fracture constitutive behavior (Eq. (33)) is represented by its corresponding strain energy:

When a segment of a fracture is bonded, the continuity of relative distance described as Eq. (36) is enforced by using the penalty method:

Fig.13. Geometric representation of open, bonded and sliding contact states for elements divided by intersecting fractures.

wheregfrepresents the stiffness of penalty springs in normal and tangential directions.

When the two surfaces of a fracture segment are sliding, the continuity of displacement in the normal direction should be satisfied and this is realized with normal component of Eq.(40).In the tangential direction, potential energy associated with sliding state for the two surfaces corresponding to Eq. (35) is

where the vector signs of Πf，slidinganddf，srepresent the two sides of one fracture segment that may have different shearing displacements. By accounting for this together with geometric representation as shown in Fig.13,sliding(slip)along fracture surfaces can be accurately simulated.

4.2.3. Example:Coupled processes in discrete fractured porous rock

The zero-dimensional fracture model is applied to a rather soft 100 m × 100 m porous domain containing 25 arbitrary fractures.The domain is subjected to a 10 kPa traction on the top surface with a hydraulic head of 0 m. The other three boundaries are impermeable.The initial hydraulic head is assumed to beh=100 m over the entire domain. For the rock matrix, the Young’s modulus is 4 MPa, the Poisson’s ratio is 0, and permeability coefficient is 2.5 × 10-8m/s. For the discrete fractures, the initial mechanical aperture is 0,the shear and normal stiffnesses are both 1×106Pa/m, the factorfis 0.5, and the residual aperturebhris 10 μm.

Inthisexample,onlythemajorresponsestotheloadare presented as a demonstration of the zero-dimensional fracture model. The simulation ran for 1000 d with the results shown in Fig. 14. The tractionloadfromthetophascausedasignificantverticalcompaction of the fracture porous media (Fig.14b). Complex hydromechanical responses in the fractured porous media result in heterogeneous vertical compaction,with subsidence of the top surface varying from about 0.3 m to 0.8 m. This vertical compaction does not occur instantaneously,butgraduallyaswaterneedstobesqueezedfromthe fully saturated porous fractured media and out through the top surface. The heterogeneously and dynamically changing hydraulic properties along with strong hydromechanical pore-volume coupling are reflected in the heterogeneous pressure field(see Fig.14c).Finally,Fig.14a indicates the locations of fracture opening marked with red circles. Thus, this application demonstrates a case of load and drainage of a saturated fractured media that may lead to long-term large and heterogeneous ground surface subsidence.

4.3. Explicit modeling of contacting rough fracture surfaces at discontinuum asperity scales

In order to understand how fractures respond to coupled processes and thereby derive reasonable mechanical and hydraulic constitutive laws,it is necessary to model coupled processes with a detailed representation of surface geometry. In this context, two challenges are settled. First, fractures with asperities (for example at a microscale) have more complicated geometric features that cannot be simplified easily;therefore,these geometric features lead to discontinuities in each physical field (flow, transport, mechanics). Second,pore-scale processes are described by a different set of governing equations for flow, transport, and mechanics. In this section, a discontinuum asperity model with explicit representation of geometric features is developed for accurate modeling of coupled processes in fractures. This approach can be used for micro-scale analysis.

4.3.1. Mathematical statement

At micro-scale,Darcy’s law is insufficient to describe flow in the open channels of fractures, and the Navier—Stokes equation in combination with mass conservation equation are required:

where v is the fluid velocity.

In addition to the force balance, small-strain assumption, and continuum constitutive equations, translational and rotational displacements need to be considered. These processes are described by the following equations:

In addition, the force is a result of not only the internal or external load,but also the interactive forces between elements(i.e.the contact force):

Fig.14. Simulated (a) horizontal and (b) vertical displacements (m), and (c) pore pressure (Pa) at t = 1000 d.

Correspondingly, continuum constitutive laws (Eq. (4), (9) and(11)) are not sufficient to describe the stress—strain/stress—displacement relationships. Instead, contact forces are the functions of displacements/location discontinuity:

Depending on the states of contact, different constitutive relationships may apply,such as Eqs.(33)—(35).The most challenging additional equations involve calculation of contact forces between multiple blocks. Different from discrete fracture analysis where fractures are initially assumed as combinations of parallel surfaces’segments, rough fractures often contain multiple asperities with arbitrary shapes.The challenge of computing contacts is to identify when and where contacts occur among many blocks, which is made complicated by the fact that these blocks are moving,deforming, and in some cases breaking apart. In turn, motion,deformation and breakage of blocks are impacted by contact forces,which constitute a serial process. Thus, inaccurate contact calculation may lead to a completely different overall system behavior.

In order to carry out contact calculation, simplifications are generally made,by either assuming contact pairs are fixed,or using simple shapes such as spheres or rectangles to approximate them.The mechanisms involved as well as the errors caused by these geometric approximations are shown in Fig.15.

The coupling between fluid flow and geomechanics within open channels (pores, channels of fractures) is mostly in one way: mechanical deformation leads to boundary changes of these fluid channels. The coupling between fluid flow and mechanics in rock matrix and open channels is carried out by ensuring the continuity of pore pressure at the fracture surfaces between the channel Navier—Stokes flow and the porous Darcy’s flow through flux terms. Therefore, fluid impact on mechanics in the rock matrix is considered using the Biot’s equation.

Fig.15. Simplified contact calculation using (a) spheres and (b) predefined rectangles for contacts.

Fig.16. Schematic of contact calculation for micro-scale mechanical analysis.

4.3.2. Numerical implementation

In NMM(Shi,1991),an algorithm was developed that rigorously incorporates contact detection, contact enforcement and openclose iteration. As shown in Fig. 16, among a number of moving and deforming blocks, NMM first computes the possible contact blocks. This involves certain estimation because the motion and deformation of blocks may occur at any time. Setting a range of possibilities enables precise and complete detection of all possible contact blocks. Between each two potentially contacting blocks,there are many possibilities where exact contact occurs.All of these possibilities are accounted in the code. They are categorized into three possibilities: vertex-to-vertex, vertex-to-edge, and edge-toedge contacts. Edge-to-edge contact is a special case where the surfaces of two contacting blocks are parallel(e.g.discrete fractures where the surfaces are parallel). Vertex-to-vertex contacts can be transformed to vertex-to-edge contact with the possible contact pairs in Fig. 16. Criteria for identifying the contact pairs from all these possible contact pairs are included in the code.After contact pairs are identified, the contact pairs may be open, bonded or sliding, as described by Eqs. (33)—(35). For bonded and sliding states,Eqs.(34)and(35)can be used to calculate the contact forces.For the open state, however, it is assumed that there is no interaction between contact pairs (i.e. complete open space between boundaries). Therefore,there are no constraints applied.

At each time step, open-close iteration may be carried out several times until the enforced contacts reach convergence. As a result of the coupled processes, contact pairs may change; for the same contact pairs, these three states may transfer dynamically.Thus, the iterations involving detection, enforcement of contact constraints,and open-close iteration need to be executed to ensure convergence.

4.3.3. Example: Calculation of rough fractures with explicit geometric representation

By using the contact model with explicit geometric representation of the fracture asperities, an example of contact alteration along rough surfaces under the impact of load is given. In this example, a single fracture partially contacting along their rough surfaces is laterally confined with two plates and fixed on the bottom. It is assumed that there is no flux exchange between the fracture channel and the rock matrix on the top and bottom, and fluid dissipates from both sides. Therefore, only one-way coupling exists: the deformation and contacts of fracture asperities impact the geometry of the fluid channel.Such an assumption allows us to decouple the problem.Fig.17 shows a mechanical simulation of the fracture alteration impacted by load. When load is applied(Fig.17a), the upper rough surface deforms and moves toward the bottom rough surface until they fully contact (Fig. 17b). This example involves large deformation (the upper surface), dynamic change of contact pairs (such as contact pairs alteration between the left plate and the upper surface),and contact states transferring(such as the contact state transferring from open to bonded between the upper and bottom surfaces).

Fig.18a—c show results of vertical displacement,vertical stress,and shear stress at the steady state. Stress concentration occurs at both contacting asperities, which can lead to different responses depending on the materials, such as plasticity, growth of new fractures, or coupled thermal and chemical responses. Such an extreme case demonstrates the capability of NMM for tackling problems of contacts between multiple bodies and/or asperities.From this example,it can be seen that contact alteration could lead to dramatic structural changes of fractures, which further impacts upscale physical(mechanical and flow)features of the fractures.On the other hand, localized stress concentration can lead to further geometric and geophysical changes of the fracture. Therefore, this modeling capability is essential for understanding fundamental behavior of fractures at micro-scale.

Fig.17. Modeling of a rough fracture: (a) before load and (b) after load reaching steady state.

Fig.18. Results of (a) vertical displacement, (b) vertical stress, and (c) shear stress at the steady state.

Table 1Overview of NMM modeling of coupled processes in fractures at multiple scales.

5. Discussion

It is quite challenging to rigorously simulate coupled HM processes in fractured geological media because of computational geometry associated with fractures at different scales. Because of limitations in the current models,there is a gap insofar,as coupled processes at discrete fracture scales could not be analyzed without consideration of intersections and shearing of the discrete fractures that might be dynamically evolving as a result of coupled processes.In addition, numerical modeling of coupled processes has never been attempted at micro-scale when the contact evolution of discontinuities is important. In this study, these computational challenges are tackled by developing a comprehensive numerical model with NMM for analyzing coupled processes in fractures embedded in porous geological media at multiple scales. Depending on the ‘relative scale’of fractures(the scale of fractures relative to the scale of interest of the problem or research),the fractures are categorized into three different scales: dominant fracture,discrete fracture, and discontinuum asperity scales. For these different scales,different governing equations as well as fracture constitutive behaviors are applied in terms of HM coupling. Correspondingly,numerical implementation varies between different scales,including implementation of those different governing equations,constitutive relationships, and HM couplings. With all these components, the NMM model presented herein is able to simulate fractures ranging from micro-scale to reservoir scale (summarized in Table 1).

For dominant fractures, Biot’s equation in combination with nonlinear constitutive relationships is used.Such a system involves both direct coupling in rock matrix and fractures, and indirect coupling in fractures.As fractures are geometrically represented by solid elements, intersections of fractures are intersections of solid elements, which is straightforward to treat. Shearing of fractures can be implicitly realized with different nonlinear laws (Eq. (11)).The most important requirement of a numerical model for dominant fractures is to accurately represent their nonlinear behavior.To accomplish this, an implicit approach that accounts for their strain energy was developed and verified.

For discrete fractures,Biot’s equation can be used for describing fluid mass conservation and balance of force in the rock matrix,and the associated direct coupling. However, the fractures, which may be oriented or intersected arbitrarily, are very thin and may be open, bonded or sliding dynamically. These fractures have flux exchange with rock matrix. Fluid flow in these fractures, which is described by a reformulated cubic law, is highly sensitive to mechanical changes (open, bonded or sliding). Similarly, the mechanical changes are highly sensitive to the fluid flow, thus a twoway indirect coupling should be considered. In order to account for such complex behavior,a zero-dimensional fracture model was developed by considering fractures as boundaries of solid rock matrix. Fluid flow in fractures and real-time flux exchange with the rock matrix are implicitly considered. Permeability is updated each time as a function of mechanical aperture, depending on the mechanical states. The mechanical states of each fracture segments are rigorously considered in three types: open, bonded and sliding with different constraints and constitutive behaviors. Such a complex system is complicated further with fracture intersections and shearing. By using a tree-cutting algorithm with discontinuous surfaces approximating each fracture, each fracture intersection is able to be calculated considering contacts between each two sides around the intersection. As each contact pair (two surfaces of each fracture segment) is updated at each time step,shearing along a fracture segment can be explicitly calculated. As each fracture is discretized into several line segments and these segments may have different contact states, complex behaviors such as shear dilation or uneven opening of a fracture can be calculated.

For fractures where the geometry of asperities cannot be simplified (e.g. at a microscale), a discontinuum asperity model is developed that explicitly represents the boundaries of asperities.For matrix of fractures, Biot’s equation is used. But for the open channel bounded by the rough surfaces of a fracture,Navier—Stokes equation in combination with conservation of mass is recommended to calculate fluid flow. Direct coupling occurs at the rock matrix.However,in the fluid channel,the coupling is majorly in one way:mechanical deformation and contacts impact the geometry of the fluid channel and thus the flow. In this situation, it is more reasonable to treat the hydrological and mechanical process as decoupled. Intersections of fractures are explicitly represented as for a single fracture.As for shearing in fractures,it occurs mostly at the asperity scale rather than at the scale of a single fracture. The challenge of modeling this micro-scale behavior is to capture when and where contacts occur. By using a rigorous contact detection algorithm, fracture alteration and shearing due to severe contact alteration are able to be simulated.

6. Summary and perspectives

In summary,modeling at different scales requires a different set of governing equations, constitutive behaviors, and geometric representations.By defining relative scales,it is possible to use the same type of model at different scales as long as the physical behaviors are described well.Between those differences for different relative scales, computational geometry plays an important role and provides the basis for numerical modeling, determining discretization to be used, accuracy for modeling of fracture intersections, and capability to treat complex behaviors such as shearing,or multi-body contacts.To date,due to the limitations of computational geometry for describing 3D fracture networks in porous rock and the limitations of geometric representation of contacts among multi-body systems, 3D computation of fracture modeling at discrete fracture or micro-scale is still rare. Development of 3D computational geometry targeting at these two problems will be essential and promising for 3D fully coupled analysis of fractures at multiple scales for advancing fundamental understanding and optimized control of energy recovery from fractured geological systems.

Declaration of Competing Interest

The authors wish to confirm that there are no known conflicts of interests associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgements

This work was supported by Laboratory Directed Research and Development(LDRD)funding from Berkeley Lab,and by the Spent Fuel and Waste Disposition Campaign, Office of Nuclear Energy of the U.S.Department of Energy,and the Office of Science,of the U.S.Department of Energy under Contract No. DE-AC02-05CH11231 with Berkeley Lab.This work was partially supported by Open Fund of the State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences (Grant No. Z017004). Editorial review by Dr. Carl Steefel at the Berkeley Lab is greatly appreciated.