Modelling stress-dependent single and multi-phase ﬂ ows in fractured porous media based on an immersed-body method with mesh adaptivity

This paper presents a novel approach for hydromechanical modelling of fractured rocks by linking a ﬁ nite- discrete element solid model with a control volume- ﬁ nite element ﬂ uid model based on an immersed-body approach. The adaptive meshing capability permits ﬂ ow within/near fractures to be accurately captured by locally-re ﬁ ned mesh. The model is validated against analytical solutions for single-phase ﬂ ow through a smooth/ rough fracture and reported numerical solutions for multi-phase ﬂ ow through intersecting fractures. Examples of modelling single- and multi-phase ﬂ ows through fracture networks under in situ stresses are further presented, illustrating the important geomechanical e ﬀ ects on the hydrological behaviour of fractured porous media.


Introduction
Discontinuities are ubiquitous in crustal rocks as a result of geological processes such as faulting and jointing. Fractures often dominate the strength [1] and deformation properties [2] of geological formations. Interconnected fractures can also serve as conduits or barriers for fluid and chemical migration in the subsurface [3,4]. Understanding the important role of fractures in the hydromechanical processes of geological media is a challenging issue which is relevant to a variety of engineering applications, including oil and gas exploitation, mining operation, geothermal production, CO 2 sequestration, and geological disposal of radioactive waste. Thus, more efforts are required to investigate the fluid flow behaviour in subsurface fractured rocks subjected to complex geomechanical effects (e.g. in situ stresses and tectonic deformation). It is considered to be necessary to understand the properties and mechanisms that govern this problem at different scales, e.g. the scale of individual fractures and that of the fracture network.
The hydromechanical behaviour of natural fractures is significantly affected by the surface roughness properties. The opposing rough surfaces are in contact through some of the asperities which control the mechanical deformation of the fractures. The separation of the fracture walls (i.e. aperture) determines the fluid flow and transport properties. However, the apertures that form channels for fluid flow can be highly variable due to the normal and/or shear deformation of fractures under in situ stress conditions. Attempts have been made to predict the behaviour of fluid flow in single fractures. The classic parallel plate solution of the Navier-Stokes equations led to the 'cubic law' to describe incompressible flow through two idealised, smooth walls [5]. The validity of the cubic law for describing fluid flow through rough fractures under normal confinement has been examined by Witherspoon et al. [6]. Tsang and Witherspoon [7] further studied the effects of fracture surface roughness and normal stress on fluid flow through single fractures and suggested that the fracture aperture for cubic law calculation can be replaced by an equivalent hydraulic aperture. Barton et al. [8] developed an empirical formulation to describe the strength, deformation, and conductivity of rough fractures based on extensive experimental data. More recently, Pryak-Nolte and Nolte [9] presented a universal scaling relationship between fluid flow and fracture stiffness to capture the aperture distribution and flow channelling within single rough fractures.
The role of the topology and connectivity of fracture networks in fluid flow has also been well recognised. de Dreuzy et al. [10][11][12] presented a series of numerical studies of flow in a fracture network with different distribution of fracture locations, lengths and apertures. Davy et al. [13] investigated and showed a strong dependency of the flow behaviour in multi-scale natural fracture systems to the fracture length distribution and transmissivity distribution. Sanderson and Zhang [14], using 2D natural fracture patterns, investigated the role of changes in the stress state of the network in varying the overall conductivity of a network as a result of localisation of flow and deformation. Min et al. [15], and Baghbanan and Jing [16] studied the effect of stress on the permeability of a fractured rock with the matrix rock as- finite-discrete element method (FEMDEM), which can capture the heterogeneity of stress distributions, interaction of matrix blocks, deformation of rough fractures, and propagation of new fractures [17][18][19][20][21].
It is noted that attempts at modelling the stress-dependent fluid flow in fracture networks have been largely focused on single-phase fluid flow problem. However, consideration of multi-phase flow in fractured rocks is important for many practical applications, such as enhanced oil-recovery (EOR) and CO 2 storage, where interaction between water, gas and/or oil is prevalent. Thus, developing a robust numerical approach to study the complex geomechanical and multi-phase flow processes in fractured rocks is demanded. To achieve this, we need to link the geomechanical model with multi-phase flow solver to be able to capture the key hydromechanical processes (see Fig. 1).
Coupling approaches for hydromechanical coupling can be generalised into two main categories: implicit and explicit coupling approaches. In the framework of implicitly coupled methods define solutions where the fluid governing equations and solid governing equations are solved together in a tightly-coupled (or 'monolithic') formulation for the flow and mechanical properties [22,23]. Implicit approaches with a monolithic formulation can be more computationally expensive per time iteration [24] but have the advantage of being able to model large time-steps. Explicit coupling methods refer to frameworks where the fluid and solid governing equations are solved in a time-marching, iterative manner, communicating coupled fields between the reservoir model and geomechanics solver. Iterative approaches (such as [25][26][27]) are often simpler mathematically and easier to implement with existing computational models. However, they are often associated with longer computational simulation time to find a solution due to the necessity of having small time-steps to achieve convergence of two or more models in the coupling scheme [24,22].
Using the immersed body method alongside a conservative mesh-tomesh projection scheme, an iterative coupled modelling scheme can be easily implemented with existing numerical frameworks and tools (e.g. commercial or research codes embodying decades of complex physics code development in the separate disciplines for solids and for fluids) in much less time than the implementation of a full, tightly-coupled formulation. Key aspects of the physics that govern single and multi-phase fluid flow in fractured media can be calculated on separate non-conforming grids and represented conservatively in an extended fluid-domain [28,29] with the use of a novel ring-mesh surrounding the solid domain. Additionally, other characteristic processes that are known to contribute to flow dynamics (such as fracture wall roughness-induced dilation and stiffness) can be easily integrated into the FEMDEM solver without the need to modify a monolithic coupled formulation.
Additionally, use of an iterative scheme based on immersed body approach and projection between non-conforming grids can enable the modelling of large-scale domains with geomechanics. Modelling of localised regions of interest in the full domain without the need to solve for coupled interactions everywhere, as in a full tightly-coupled formulation, can be accommodated with this coupling approach.
This new alternative approach combines the capability of the FEMDEM method in capturing geomechanical behaviour of rock matrix and natural fractures [19], and that of the control-volume finite-element method (CV-FEM) in simulating multi-phase flows with sharp interfaces well preserved by mesh adaptation [30]. The remainder of the paper is organised as follows. The governing equations and models used in the hydro-mechanical modelling are outlined in Section 2, followed by a series of validation tests for flow through single smoothor rough-walled fractures in Section 3. In Section 4, examples of simulating stress-dependent single-and multi-phase flow through fracture networks are further presented. Finally, a brief discussion is given and conclusions are drawn.

Solid model
The mechanical behaviour of fractured rocks is solved based on the combined finite-discrete element method (FEMDEM) [31], which can capture large strain deformation, multi-body interaction, fracturing and fragmentation [32,17,33,34]. The 2D FEMDEM model uses an unstructured, fully-discontinuous mesh of three-noded triangular finite elements with their interfaces connected by four-noded joint elements. There are two types of joint elements: unbroken joint elements embedded inside the matrix and broken joint elements placed along existing fractures [19]. The joint elements (both broken and unbroken) are inserted between the finite element pairs before the numerical simulation and no remeshing is performed in the solid computation. The equation of motion of each finite element is given by where M s is the lumped nodal mass matrix, u s is the velocity vector, t is the time, F c is the contact force, F int is internal nodal forces, F l is the boundary load, F p is the pore fluid pressure, and F b is the body force. The deformation of the bulk material is captured by linear-elastic constant-strain finite elements with the continuity constrained by the bonding forces of unbroken joint elements [35]. The interaction of discrete matrix bodies is calculated based on the penetration of finite elements via broken joint elements [36]. The elasto-plastic fracturing of geological rock media is modelled by a smeared crack (i.e. cohesive zone) model [35] that can consider the non-linear stress-strain behaviour in the plastic zone ahead of a fracture tip. The propagation of new fractures is further represented as the transition of unbroken joint elements to broken ones in the unstructured grid.
The empirical joint constitutive model proposed by [37,8,38] has been implemented into the FEMDEM model [19]. The non-linear closure of a fracture under compression is characterised by a hyperbolic equation [37]: where ν n is the current normal closure (mm), σ n is the local normal stress (MPa) that is derived from the Cauchy stress tensor of adjacent finite elements, k n0 is the initial normal stiffness (MPa/mm), and ν m is the maximum allowable closure (mm). Values of k n0 and ν m are given by [37]: where a 0 is the initial aperture (mm), JRC is the joint roughness coefficient, and JCS is the joint compressive strength (MPa). Both JRC and JCS are scale-dependent parameters [39] and their field-scale values, i.e. JRC n and JCS n , can be estimated using [8]: where L n is the field-scale effective fracture length (i.e. the length between fracture intersections) calculated through a binary-tree search of connected broken joint element chains [19], JRC 0 and JCS 0 are measured based on the laboratory sample with length L 0 . For unweathered fractures, the initial aperture a 0 may be estimated as [37]: During the shearing process under a normal compression, fractures first contract due to the compressibility of asperities and then dilate with roughness damaged and destroyed [8]. Dilational displacement can be related to the shear displacement based on an incremental formulation given by [38]: where v d s is the increment of normal displacement caused by shear dilation, u d is the increment of shear displacement, and d mob is the mobilised tangential dilation angle given by [38]: where M d is a damage coefficient given by [40]: The mobilised joint roughness coefficient JRC mob can be calculated using a dimensionless model [8] based on the ratio of the current shear displacement u c to the peak shear displacement u p , which is given by [8]: The coupled normal and shear deformation can thus be modelled by an incremental formulation: The mechanical aperture a m is derived by combing the effects of mesoscopic opening (induced by fracture network deformation and explicitly resolved in the FEMDEM) and microscopic closure (controlled by fracture roughness and implicitly captured by the JCM) as given by [19]: where w is the mesoscopic normal separation of the opposite walls of broken joint elements in the deformed FEMDEM mesh, and v is the microscopic accumulative closure derived from the incremental formulation. The first part of the piecewise function corresponds to the scenario that the broken joint element is mesoscopically opened, while the second part models the condition that the two opposite walls of the fracture are in contact at the FEMDEM discretisation scale. The hydraulic aperture a h defined as an equivalent aperture for laminar flow is further calculated based on an empirical relation with the mechanical aperture a m [38]: A linear interpolation is used to determine the value of hydraulic aperture in the transition phase, i.e.
Given the inertia effect of fluid flow is negligible and the ratio of joint segment length to aperture is high [41], the permeability of the local broken joint element k f can be calculated based on the cubic law [6]: Its permeability tensor k f is thus given by: where θ is the angle between the mid-line of the broken joint and the xaxis (Fig. 2).

Fluid model
The incompressible and isothermal single and multi-phase flow through fractured porous media is solved based on the control-volume finite-element method (CV-FEM) [30,[42][43][44]. The equations below are presented in the form for multiple phases, but can be easily simplified for single-phase flow. The equation of motion of fluid flowing through porous media is governed by Darcy's law: denotes different fluid phases (n is the total number of phases), u α is the phase saturation-weighted Darcy velocity of phase α p , α is the pressure, and s uα is a source term (e.g. gravity and/or capillary pressure). The term on the left-hand side σ α is treated as an 'absorption term': where K is the permeability tensor, μ , r α α K and S α are the relative permeability, viscosity and saturation of phase α, respectively. The mass conservation equation is given by: where ϕ is the porosity and s cty is the volume source term. The phase saturations are further constrained by: Time discretisation of the problem domain is based on the 'θ-scheme' which permits a smooth switch from θ = 0.5 (Crank-Nicholson, semi-implicit) to 1 (backward-Euler). The backward-Euler (θ = 1) discretisation is used here for simplicity. Further details about the fluid model can be found in [30].

Immersed-body approach
The solid domain is linked to the fluid field using an immersed-body approach [29]. The solid mesh discretised independently is fully  immersed in a 'super' fluid mesh with the information of the two fields exchanged via a 'ring' mesh surrounding the solid (at both the model boundary and fracture walls) ( Fig. 3a and b). The ring mesh is generated or re-generated (if deformation occurs) automatically along the fracture walls or the solid boundaries (Fig. 3c). The thicknesses of the ring mesh elements along the broken joints are constrained to the scale of the local hydraulic aperture to accurately characterise the flow complexity near/in the fracture, as shown in Fig. 3c.
A conservative Galerkin interpolation [45] is used to project local fields between the solid domain mesh (Fig. 4a) and the fluid mesh (Fig. 4b) using the generated ring mesh surrounding the solid domain. The interpolation can project the geometry of the solid and ring meshes, as well as the local permeability of fractures to the fluid solver. The overlap between the projected meshes and the fluid mesh (Fig. 4c) is computed to scale the field values for conservation purposes: wher Ω is the domain and q is a particular field for projection (such as nodal coordinates, permeability, etc.). In this approach, the projection of permeability from the ring mesh to the fluid mesh is thus governed by: where k f is the fracture local permeability tensor projected through the aperture-scaled ring mesh, N i f and N i r are the basis functions of the fluid  and ring mesh respectively, and n f and n r are the total number of elements in the fluid mesh and ring mesh, respectively. The objective of the projection is to minimise the L 2 norm of the error [46]. The technique of adaptive mesh refinement (and coarsening) is employed in the fluid flow calculation, so that the stress-induced aperture and permeability evolution of high aspect ratio fractures can be well respected (Fig. 5). The metric for identification of fractures in the fluid solver is achieved using the projected permeability field from the solid domain. The Hessian of the field [47] is calculated to get local curvature; a high curvature means that the mesh needs to be refined (e.g. the region near a fracture where a sharp gradient in the permeability field is expected), while low curvature suggests that the mesh can be coarsened (e.g. the region inside the matrix): where M a is the metric used for adaptivity computation based on the Hessian ∊ H x x ( ), ( ) is a user-defined weight, n is the dimension of the space. The metric provides a value for the bound of the L 2 -norm interpolation error with a target value of ∊ x ( ) [48,47]. The adaptive mesh is also used to focus the resolution at the sharp interface between different fluid phases (Fig. 5c-d). More details about the mesh adaptivity method can be found in [47].

Coupling overview
The iterative coupling framework is generalised as a schematic and presented in Fig. 6.

Validation
We validate the immersed-body approach based on benchmarking tests of single and multi-phase flow through single or multiple fractures with prescribed apertures. Here, we only validate the model for simulating flow through a fracture with a prescribed aperture, while the capability for simulating stress-dependent aperture variation in single fractures has been validated in the previous work [19].

Single-phase flow through fractures
The validation test is based on a single pre-existing fracture spanning a square rock domain (0.1 m × 0.1 m), as shown in Fig. 7. The rock has a matrix porosity of 0.02 and permeability of × − 1 10 17 m 2 (i.e. 0.01mD). The solid is immersed within a larger fluid domain (0.2 m × 0.1 m). The solid and fluid properties are listed in Table 1. The simulation results of flux through the fracture will be compared with the analytical solution based on the cubic law [5]: where Q is the volumetric flux through the fracture, P in and P out are the pressures at the inlet and outlet, respectively, μ is the viscosity, a h is the prescribed hydraulic aperture, and L 0 is the sample fracture length.
Two scenarios are considered. In the first case, the single-phase fluid flows through a smooth-walled fracture with a constant aperture of different prescribed values. Significant consistency is obtained between the numerical and analytical results of the flux through the fracture (Fig. 8). In the second case, the fluid flows through single rough-walled fractures associated with different JRC 0 values, i.e. 5-20. As shown in Fig. 9, the numerical results are also closely matched to the analytical solutions.

Multi-phase flow through fractures
We also validate the proposed approach for simulating multi-phase flow through multiple fractures in a porous medium. Refined mesh is not only placed in/near fractures, but also used to capture the sharp interface between different fluid phases. The modelling example corresponds to a scenario of water-flooding into an initially oil-saturated rock embedded with an idealised system of multiple fractures. We denote water as phase 1 and oil as phase 2. The rock matrix has a porosity of 0.2 and a permeability of × − m 1 10 15 2 . The simulation results of phase  volume fraction are compared with those from Eclipse and GPRS reported by [49]. In this simulation, water is injected at the left bottom corner of a 1 m 2 domain of fractured porous rocks initially saturated with oil (Fig. 10). The water viscosity is × − 1 10 Pa·s 4 and the oil viscosity is × − 0. 45 10 Pa·s 4 . The effects of gravity and pore compressibility are neglected. Simulations are conducted with two different meshes: a fixed fine mesh with an element size of 8 mm (Fig. 11a), and an adaptive mesh with the minimum and maximum element size set as 1 mm and 100 mm, respectively (Fig. 11b). A 37 % reduction in the number of total elements in the mesh is achieved when enabling adaptive remeshing. A good consistency is found between our modelling results and the simulation results in [49] using the commercial simulator Eclipse or the discrete-fracture based GPRS model (Fig. 12) for the oil production driven by the injected water. However, as shown in Fig. 13c-d, our method can better capture the saturation evolution of the water and oil phases in the fractured porous medium. Especially, the mesh adaptivity technique allows a highly-reduced number of elements, while still well capturing the sharp interface characteristics between different phases (Fig. 13d). A reduction of × 4.5 is seen for the simulation time with use of adaptive mesh refinement, compared with the fixed mesh example of Fig. 11(a).

Fluid flow through stressed fractured porous rocks
The numerical modelling approach developed in this research has been applied to simulate single-and multi-phase flow through a natural fracture network. The fracture pattern is based on an outcrop mapped at Kilve in the southern basin of the Bristol Channel, UK (Fig. 14) [50]. This fracture network includes two distinct sets of fractures, orientated approximately 100°and 140°to the x-axis, respectively. The rock has a density of 2700   (a) (b) Fig. 11. Two different meshes used to illustrate the advantages of adaptive mesh simulations: (a) a fine non-adaptive mesh, and (b) an adaptive mesh. Note the significant reduction in the total element number by using the adaptivity technique. Fig. 12. Comparison of the cumulative oil production results calculated using the approach presented in this paper (solid line), and those using GPRS (dashed lines) and Eclipse (triangles) as reported in [49]. is used for water and oil, respectively. Capillary pressure is not included.

Stress effects on single-phase flow
The distributions of shear displacement and hydraulic aperture of the fracture network in four different in situ stress ratio cases ( = σ σ / 1 x y to 4) are shown in Fig. 15a and b, respectively. An increase in both shear displacement and hydraulic apertures is observed with the increased stress ratio. In particular, the largest increase in shear displacement is found in the joint set orientated at 140°(see Fig. 15a). A significant increase in hydraulic aperture also occurs in the same joint set due to the shear dilation effect. Due to the high contrast between the fracture and matrix permeabilities, the deformation of fractures also leads to variation in the fluid pressure field (Fig. 15c-d). More localised pressure distribution is observed in the high stress ratio case. Fig. 16 Fig . 13. Simulation of water flooding through an initially oil-saturated fractured porous medium. Two previously published results using (a) the commercial software Eclipse and (b) the discrete-fracture based GPRS model [49]. Our simulation results for two different meshes using (c) a very fine, fixed mesh, and d) an adaptive mesh. Note the significant reduction in the total number of elements by using the adaptivity technique. further gives the probability density distribution of shear dilation and hydraulic aperture for all broken joint elements (representing all preexisting and new fractures) under the four different in situ stress ratio conditions. Shear dilation is observed to be noticeably increased in higher stress ratios, especially for the highest stress ratio of 4 (see the maximum dilation of 0.027 mm (A in Fig. 16a) and 0.29 mm (B in Fig. 16a) for the stress ratio of 1 and 4, respectively), which therefore results in increased hydraulic apertures (Fig. 16c). The equivalent permeability of the fractured rock in the x and y directions, i.e. − k eq xx and − k eq yy , are derived from the simulation of single-phase flow through the system subjected to different in situ stress conditions (Fig. 17). With the increase of the stress ratio, − k eq xx decreases slightly first due to the contraction of the fractures under increased mean stress level, and then increases as a result of shearing and dilation that mainly occurred in the fracture set oriented at 140°. Note that the variation of − k eq xx is not considerable. However, − k eq yy exhibits significant increase with the increase stress ratio, because of the dominant role of the 140°oriented fracture set in flow in the y direction. In all the stress conditions, − k eq yy is higher than − k eq xx .

Stress effects on multi-phase flow
To further demonstrate the capability of the immersed-body method for modelling more complicated multi-phase flow through geomechanically-constrained fracture networks, simulations of water flooding into an initially oil-saturated fractured rock are conducted. The fractured rock is loaded by a series of different in situ stress conditions, i.e.  injection rate is set at the left boundary of the domain and a fixed pressure is defined for the right boundary, while the top and bottom boundaries have a no-flow condition. Fig. 18 shows the evolution of the phase volume fraction field in the fractured rock under the different in situ stress ratio conditions. In addition to using adaptive mesh refinement around the fractures in the fluid model, in order to reduce aritificial saturation, adaptive mesh refinement is also employed to capture the sharp interface between water (phase 1) and oil (phase 2), as shown in Fig. 19. It can be seen that the injected water quickly penetrates through the domain via major spanning fractures, before gradually migrating into small fractures and the porous matrix, during which the oil is progressively pushed out from the right boundary (see Fig. 18). It seems that the stress conditions with > σ σ x y ( Fig. 18b and d) tend to exhibit more significant channelized flow than those with < σ σ x y ( Fig. 18c and e). For example, see the greater penetration of water into the fracture one third way down the left hand boundary of the domain in Fig. 18d. This is because the spanning fractures are less closed when > σ σ x y . Fig. 20 gives the cumulative production of oil from the system driven by the injected water. The fractured rock under the isotropic stress condition produced much less oil, because of the smaller amount of oil stored in the isotropically compressed fractures. However, the fractured rock under the anisotropic stress conditions accommodated more oil in the fractures associated with larger apertures as a result of shear dilation, which thus leads to more oil production. The diversion of the production curve of the isotropic stress case from others of the anisotropic cases at the point where the pore volume of injected water is 0.12 is considered to be due to the effect of there being more highly compressed apertures under the isotropic stress loading. It can also be noted that for the cases with > σ σ x y , more oil is produced than for the cases with < σ σ x y , because fractures of Set 1 are less compressed and fractures of Set 2 are highly sheared and dilated when > σ σ x y . However, it can be expected that Set 2 is also significantly dilated for the case = σ σ 3 y x , and hence, production of this case is also high.

Discussion
The paper set out to introduce a novel method for modelling hydrological behaviour of fractured porous media subjected to geomechanical effects. The model can capture the fracture behaviour at both the roughness scale, e.g. fracture closure and shear dilation, and the network scale with complicated connectivity properties. The methods and validation results shown for single fractures and idealised networks together with examples of complex natural fracture networks suggest that this generic framework for combining fracture and matrix flow calculations has the potential to improve flow prediction accuracy over more traditional methods used to account for in situ stress effects. The linkage of the solid and fluid domains through a conservative field projection via the ring mesh combined with the adaptive mesh refinement permits very accurate simulation of fluid flow through the high aspect ratio natural fractures associated with stress-dependent variable apertures. The multi-phase fluids mesh occupies the entire computational domain and by reference to the solid domain with its highly resolved opening, closing and shearing fractures, carries the solid volume fraction of each element increasingly accurately as the multi-phase mesh refines in the vicinity of fractures. Thus, the realistic behaviour of the stressed solid skeleton containing the fractures can be well respected in the fluid flow simulation. Another key feature of this research is the study of stress-dependent multi-phase flow processes in fractured rocks. Most of previous works on stress-dependent subsurface flow [14][15][16][17][18]21,32,33] were mainly focused on single-phase flow problems. However, the stress effects on multi-phase flows can also be significant, as illustrated in the water flooding simulation in this paper. The benefits of the mesh adaptivity technique are also manifest in the capturing of the sharp interface between different fluid phases. Use of Discontinuous-Galerkin (DG) discretisation methods in the CV-FEM formulation, alongside adaptive mesh refinement, artificial diffusion of saturation across the fracture and solid matrix can be significantly reduced with small computational overhead [51]. This approach is an improvement of Su et al. [51], taking into account stress-dependent behaviour of fractures at the local and network scale, and is an alternative technique to conventional lower-dimensional representation of fractures as used in [52,53]. Modelling of fractured porous rock systems with our approach can have very useful implications for many engineering applications, such as oil production in fractured carbonate reservoirs and CO 2 storage in geological formations [54].
In the current work, the effects of pore structure compressibility and variation in matrix permeability under stress loading were not taken into account, which may be considered minor given the high contrast between fracture and matrix permeabilities [55]. Furthermore, in order not to over-complicate the explanation of the immersed-body coupling method, this paper emphasises only the one-way coupling scenarios. That is to say, the solid skeleton deformation and aperture development computed for different in situ stress scenarios are first derived from the solid solver as an approximate equilibrium state and then, for simplicity, the solid is constrained thereafter to be fixed.
In the examples, there is therefore no opportunity for the fluid pressures associated with the flows to dynamically alter the solid skeleton geometry. Of course, such changes would be expected if the pore fluid pressure has a significant transient variation that causes hydraulic fracturing and/or pore structure dilation/compaction. Thus, more effort is needed to incorporate the influence of fluid pressure on fracture deformation and propagation and the effect of poroelasticity [54], which will be the key focus of our future research. One of the purposes of writing this paper, whereby a general iterative coupling framework is introduced, in spite of reserving twoway coupled applications for subsequent papers, is to emphasise the range of one-way coupled reservoir flow problems that also benefit from this adaptive mesh optimisation approach. The approach presented enables the fluid mesh to automatically identify the location of fractures and focus extra mesh resolution there. Furthermore, the key role of incorporating geomechanics in coupled frameworks for fracture flow is to account for local in situ stress modification and the associated fracture aperture dependencies. In this context, accounting for the roughness of the fracture walls in the hydraulic aperture calculation is important. Hence, the inclusion of a validation of the implementation of the empirical fracture roughness model within this coupled framework. Actually, capability of the immersed-body coupling method developed by our research group for handling the two-way coupling problems has already been demonstrated in the context of solids interacting with inertia-dominated flows [28,56]. In the next paper, we will 'switch on' the two-way coupling function and show how it can be applied to simulate the complex coupled solid-fluid interactions involved in hydraulic fracturing.

Conclusions
A novel approach for hydromechanical modelling of fractured porous media based on an immersed-body method has been presented in this paper. The approach has been used to investigate single-and multi-phase flow in complex fracture networks under geomechanical effects. The key achievements are summarised as follows: (1) Use of adaptive unstructured mesh refinement, and coarsening to handle flows along high aspect ratio geological structures, ensuring that artificial diffusion is low or negligible when large variations and discontinuous variations in hydrological rock properties and fluid velocities are present. (2) Illustration of an immersed-body technique that allows leak-off to be captured automatically for all cases including whether matrix permeabilities and porosities and fluid matrix storage are low or high. (3) As an alternative to the modelling approach to fracture matrix flow that is based on a lower-dimensional representation of fractures, the new approach generates a local anisotropic permeability tensor (informed by fracture orientation) to represent the fracture flow while retaining the geometry of linking and dead-end fractures, and their variable apertures. (4) Application to realistic locally-varying fracture apertures based on a natural fracture network i.e. a mapped outcrop of a fracture network (the method can equally well be applied to a statistically generated DFN) that is subsequently brought to equilibrium with an applied in situ stress state. (5) A preliminary exploration of multi-phase behaviour in water injection scenarios for a fractured oil-saturated rock mass with anisotropic permeability arising as a consequence of local in situ stress modification of the fracture system is also presented. (6) The approach presented can be readily extended to dynamic twoway coupling simply by allowing an iterative loop for the updating and passing changes in fluid pressure as a result of hydro-geomechanical interactions.