Multiphase ﬂow simulation through porous media with explicitly resolved fractures

Accurate simulation of multiphase ﬂow in fractured porous media remains a challenge. An important problem is the representation of the discontinuous or near discontinuous behaviour of saturation in real geological formations. In the classical continuum approach, a reﬁned mesh is required at the interface between fracture and porous media to capture the steep gradients in saturation and saturation-dependent transport properties. This dramatically increases the computational load when large numbers of fractures are present in the numerical model. A discontinuous ﬁnite element method is reported here to model ﬂow in fractured porous media. The governing multiphase porous media ﬂow equations are solved in the adaptive mesh computational ﬂuid dynamics code IC-FERST on unstructured meshes. The method is based on a mixed control volume – discontinuous ﬁnite element formulation. This is com-bined with the P N+1 DG-P N DG element pair, which has discontinuous (order N + 1) representation for velocity and discontinuous (order N) representation for pressure. A number of test cases are used to evaluate the method’s ability to model fracture ﬂow. The ﬁrst is used to verify the performance of the element pair on structured and unstructured meshes of different resolution. Multiphase ﬂow is then modelled in a range of idealised and simple fracture patterns. Solutions with sharp saturation fronts and computational economy in terms of mesh size are illustrated.


INTRODUCTION
The flow of multiphase fluids in fractured porous media affects a large number of industries that exploit geological resources, such as petroleum, CO 2 sequestration, geological disposal of nuclear waste, mining and tunnelling. These flows are also of scientific interest with respect to fracturing, faults, stress interaction and underground flows. For reservoir engineers and hydrogeologists concerned with modelling fluid transport, characterising fractured rocks is one of the main challenges, because a good description of the network connectivity and apertures in the subsurface is still elusive (Faybishenko et al. 2000;Bonnet et al. 2001;Berkowitz 2002). Another major challenge comes from the complexity of the coupled physical processes involved in fractured porous systems carrying multiphase fluids, and the need to predict the response when the system is perturbed by the action of engineers. An important coupled process that is often neglected in current modelling approaches takes place at the scale of the fracture aperture. The solid skeleton of rock is described by the fracture wall positions, which are held in equilibrium by the balancing combination of the local total stress field near the fracture walls and the local fluid pressure in the fractures, both of which can be subject to change. If they do not change, the stress-dependent fracture apertures are considered constant in what are sometimes described as stress-dependent permeability models, and the hydraulic apertures and effective flow properties are constant. However, the micro-seismicity associated with [The copyright line for this article was changed on 23 March after original online publication]. man-made perturbations of reservoir fluid pressures, or in situ stress relief in response to rock excavation, is clear evidence that the solid skeleton cannot always be regarded as static; furthermore, quasi-static changes in pore fluid pressures resulting from interventions by man can introduce imbalances leading to changes in flow properties. It is these complex coupled processes, in which dynamic unlocking of the fracture walls (e.g. in shear, or a fluid injection driven fracture propagation, i.e. hydro-fracturing) causes transient, potentially hazardous or beneficial responses that are proving especially challenging to model.
Fracture wall behaviour and fracture apertures are highly stress-dependent and must be considered explicitly with detailed geometric representation if the important physical phenomena are to be captured, especially those involving transient flow through high-permeability pathways on the timescale of seismically observed reactivation or fracture propagation events. There is an increasing recognition of the importance of geomechanics in reservoir behaviour and the fact that, for many important phenomena, hydromechanical coupling effects, which include not just how the flow is guided by the solid skeleton, but also how the fluid can alter the solid skeleton in a full two-way coupling, must also be modelled (Rutqvist & Stephansson 2003;Zoback 2010). This paper aims to contribute towards this latter challenge, by developing computational methods that better capture key aspects of the physics associated with multiphase flow in and around fractures for which fracture wall location change is permissible, including fluid transfer between fracture and matrix, and the response when the system is perturbed by the action of engineers. We outline a modelling methodology that has the capability to tackle the particular problems posed by the extreme aspect ratio of fractures and by the inherent solid-fluid coupling of fluid flow through fractured media, through application of an adaptively refining anisotropic mesh, as developed in the Imperial College Finite Element Reservoir SimulaTor (IC-FERST, see Jackson et al. 2013), a module of the general purpose CFD code fluidity developed at Imperial College (see AMCG 2014). This dynamically optimises mesh resolution to capture sharp saturation changes at fracture boundaries, as well as the moving saturation front within fractures and matrix. The method includes explicit representation of fractures as 3D domains with high permeability and extreme aspect ratio, bounded by fracture walls with potentially complex geometries embedded within a lower permeability matrix.
The governing equations for multiphase flow through porous media are solved with triangular/tetrahedral elements using a mixed control-volume finite-element (CVFE) method and the P N+1 DG-P N DG element pair, which has a discontinuous, order N+1, representation for velocity and a discontinuous, order N, representation for pressure (Jackson et al. 2013). The discontinuous representation for pressure allows us to use control volumes that are discontinuous across fracture walls, allowing better representation of sharp saturation changes between fracture and matrix. The underlying mass balance equations are solved in control volume (CV) space, and finite elements (FE) are used to obtain the high-order fluxes on CV boundaries (Jackson et al. 2013). A related family of elements, the P N DG-P N+1 element pair, originally introduced by Cotter et al. (2009) andHam (2011), is also used here for comparative purposes, as it has similar balance preserving properties for porous media flows with the P N+1 DG-P N DG element pair. These features are implemented for porous media flow in the open-source IC-FERST (Jackson et al. 2013) and, for more general CFD applications, in fluidity (AMCG 2014). Our new approach is compared against previous published methods for modelling fracture flow in the Discussion.

Governing equations
Assume a volume of an incompressible porous medium with porosity φ that does not vary in time, containing N k immiscible, incompressible fluid phases. The volume fraction of each phase is where the subscript, k 2 ½1; N k , denotes phase, and S k is the saturation of that phase such that Darcy's law states that for a phase k in a porous medium, where q k is the Darcy velocity, k rk is the relative permeability of the phase and is a function of S k , l k is the phase dynamic viscosity, K is a rank two tensor describing the permeability of the porous medium, p is the pressure, and ς uk is a source term (such as gravity forces and capillary pressure) associated with the force balance. The Darcy velocity can be rewritten in terms of the interstitial velocity of the fluid then we can rewrite Eq. 3 as: Conservation of mass for each fluid phase implies: where Q k is a volume source term.
The weak form of the discretised force balance equations for time level n and phase k are thus of the form: where the finite element pressure and velocity fields are given by: in which n is the normal to the element, and N p and N u are the degrees of freedom for the FE pressure and velocity representations, respectively (see Appendix 1 for further details). Here, Ω is the volume domain with Ω E the subspace associated with velocity element E, Γ E is the surface boundary of the element E, and Γ Ωp is the boundary of the domain on which pressure is set weakly to p bc (see Fig. 1). The interelement pressurep is the finite element pressure on the other side of the boundary of the element E. The continuity equation (Eq. 6) is discretised in space with CV basis functions M i and in time with an adaptive h-method which switches smoothly from Crank-Nicolson (h = 0.5) to backward-Euler (h = 1) (see Blunt & Rubin 1992; as implemented by Gomes et al. 2008). The Crank-Nicolson method has the simplicity of a two-level time-stepping method, which is second-order accurate and unconditionally stable. However, large grid Courant numbers can result in oscillations and unphysical solutions. For each time step, a value of h is calculated at each CV interface based on a Total Variation Diminishing (TVD) criterion. TVD schemes are computationally efficient and are often employed to solve the advection equation within CVFE methods and similar formulations. To determine where to apply high-order fluxes, extrema are detected and quantified according to a normalised variable diagram (NVD). We adopt this approach because it is computationally efficient and ensures a first-order, nonoscillatory method is applied at extrema with an upper bound that corresponds to a TVD condition in 1D. The resulting equations for time-step size Dt and phase k are given by: where the saturation: The system is closed by the summation constraint, Eq. 2. More details on the equation set, discretisation methods and overall solution method are given in Appendix 1.

Models and numerical solutions
In the models presented here, structured and unstructured meshes (shown in Table 1) for the fractures and the porous media were generated with the software packages AN-SYS (De Salvo et al. 1987) and Gmsh (Geuzaine & Remacle 2009), respectively. The fractures considered here are idealised and have fixed length, aperture and orientation. Structured and unstructured triangular elements have been employed to verify the degree of mesh dependence of the numerical solutions. All models and results are reported in dimensionless form, with model length normalised to 1, permeability contrasts reported as a ratio normalised to the matrix permeability, phase saturation (A) (B) Fig. 1. Discretisation scheme. (A) Discretisation scheme for P 1 DG-P 2 element pair. Shaded areas represent degrees of freedom for saturation fields, white nodes represent degrees of freedom for the velocity field and triangles represent degrees of freedom for the pressure field. Velocity is piecewise linear on the intersections between control volumes (shaded areas) and the triangular elements (as indicated by hatchings), but discontinuous across all boundaries. Pressure is piecewise quadratic and continuous between elements. (B) Discretisation scheme for P 2 DG-P 1 DG element pair. Control volumes of both elements are outlined by dashed lines. Ω Cv is the volume of a control volume, and Γ CV is its boundary.
normalised to vary between 0 and 1 (i.e. removing irreducible and residual saturations) and time-scaled in terms of pore-volumes injected (PVI). Porosity is assumed constant, which is a simplification but does not affect the significant aspects of the results presented nor is constant porosity a required limitation of the methodology. Gravity and capillary forces are neglected; it is assumed viscous forces dominate. A Corey-type relative permeability model is chosen for both wetting (water) and nonwetting (oil) phases given by The viscosity ratio of the two fluid phases is constant with the resident nonwetting ('oil') phase 10 times more viscous than the injected wetting ('water') phase. Twophase incompressible flow is simulated assuming the wetting phase ('water') is injected at a constant interstitial velocity (of dimensionless value one) over one model boundary, into a domain initially saturated with the nonwetting phase ('oil'), and fluids are produced at constant (zero) pressure over the opposite boundary; all other boundaries have no flow across them. The simulations are designed to be analogous to an oil production scenario in which water is injected to displace oil towards the production well(s). Model 1 contains a low-permeability barrier; models 2-5 include high-permeability fractures with varying geometries. Fluid transport is tracked by monitoring the volume fraction of the injected wetting phase in each element. The time step is constant and corresponds to 4 9 10 À4 PVI for Model 1 and 2 9 10 À5 PVI in all other simulations.

Model 1
Model 1 contains a low-permeability barrier (Fig. 2B). Space is discretised with several different unstructured mesh resolutions ( Fig. 2A). The results from several simulations are compared to illustrate the differences between the continuous and the novel discontinuous method formulated here.
For the continuous method model, as fluid flows through the rock domain, there is no sharp saturation interface between the lower permeability barrier and the more permeable matrix (Fig. 3). The wetting (water) saturation varies smoothly towards, and within, the barrier. With decreasing element size around the barrier (from mesh A to mesh D), the interface between the inclusion and the background becomes sharper. These results confirm that a considerably refined mesh is required to accurately model flow around a barrier when a continuous numerical method is employed. In this case, the wetting phase spreads into the low-permeability barrier because it is defined CV-wise and the CVs span the boundary between the barrier and the matrix. As the mesh resolution increases and the CVs become smaller, the effect becomes less significant.
In contrast, as shown in Fig. 4, even the coarsest mesh (mesh A) preserves the expected sharp saturation discontinuity between the low-and high-permeability regions when the novel discontinuous method presented here is employed; recall capillary forces are neglected, so only viscous forces drive the wetting phase into the low-permeability domain. Flow around the barrier is better captured with the discontinuous method than with the continuous one when the same mesh is adopted. This has an effect on produced volumes of the nonwetting (oil) phase and also on the time at which the injected (water) phase traverses the domain and breaks through at the outlet face: water breakthrough occurs after 0.24 PVI in the discontinuous model (mesh A), but later in the continuous model because of the nonphysical migration of water into the low-permeability barrier. The accuracy of the results of the continuous method with mesh D, whose elements' size is about one twelfth of mesh A, is approaching the accuracy of results of the discontinuous pressure method with mesh A. As far as computational speed is concerned, the discontinuous method is slower than the continuous one for the same mesh, by a factor of approximately 1.5 (see Table 2). Nevertheless, the discontinuous method with mesh A is significantly faster (by a factor of 3) than the continuous method with the much finer mesh D needed to achieve comparable overall accuracy. We return to the issue of computational efficiency in the discussion section.  Model 2 contains a single high-permeability fracture that traverses the model domain from injection boundary to production boundary, embedded in a lower permeability matrix (Fig. 5). The discontinuous method presented here very accurately preserves flow in the fracture despite its very high aspect ratio and minimises nonphysical leakage (in the absence of capillary forces) of the injected wetting phase into the surrounding matrix; the correct physical solution here is a sharp saturation discontinuity at the boundary of the fracture, with a much slower movement of the injected phase into the lower permeability matrix adjacent to the inlet boundary (Fig. 6). Note that flow in the fracture does not match the 1-D Buckley-Leverett solution for twophase, incompressible, viscous-dominated flow (Buckley and Leverett, 1941 as reported in Dake 1978) because there is cross-flow between the fracture and matrix, as shown by the pressure contours (Fig. 7). Close to the inlet, there is flow from the matrix into the fracture, driven by the constant injection rate across the inlet boundary and the higher permeability in the fracture. Conversely, close to the saturation front in the fracture, there is flow from the fracture into the matrix, driven by the higher pressure in the fracture arising from the differing mobility of the injected and displaced phases. This cross-flow arises only because of viscous forces and is often observed in systems of contrasting permeability (see, for example, Zapata & Lake 1981); it is not a numerical artefact in the discontinuous pressure method and does not arise from capillary forces, which are neglected here. Note that the flow of the injected phase into the matrix adjacent to the fracture is retarded relative to the location of the displacement front away from the fracture; this reflects the lateral movement of the injected phase into, and preferential flow through, the high-permeability fracture. Channelling of the injected phase into the fracture leads to very early breakthrough at the outlet face after just 0.025 PVI.

Model 3
Model 3 has two perpendicularly intersecting fractures with the same permeability ratio and model set-up as in section 3.2 (Fig. 8). As in Model 2, the discontinuous method is able to accurately capture flow through the fractures, including a small amount of viscous cross-flow into the fracture parallel to the inlet and outlet boundaries (Fig. 9). This causes breakthrough to be slightly delayed relative to the single fracture case, at 0.03 PVI.

Model 4
Model 4 has a regular five by five network of perpendicular intersecting fractures. The mesh and permeability map are shown in Fig. 10. As expected, each fracture has the same velocity and phase volume fraction at each time level (Fig. 11). At each intersection, there is a small amount of viscous cross-flow into the intersecting fracture set, driven by the (small) pressure difference between the flowing and stagnant fractures that results from the changing total mobility. Flow into the fractures parallel to the inlet and outlet boundaries causes breakthrough to be further delayed relative to the single and double fracture cases, at 0.035 PVI.

Model 5
In this model, 20 pseudo-randomly distributed fractures of variable length and orientation are inserted in the rock matrix (Fig. 12). The model domain is discretised using a fully unstructured mesh (in contrast to models 2-4) because of the varying orientation of the fractures. Simulation results for two different mesh resolutions are reported; the coarser mesh E contains 1308 elements and the finer mesh F contains 5316 elements (see Table 1). The discontinuous method is able to capture flow through the fractures with similar accuracy regardless of mesh resolution, but solutions obtained using the continuous method are strongly meshdependent (Fig. 13). In particular, the injected (water)  phase contacts more of the resident (oil) phase in the matrix as the mesh resolution decreases and penetrates a smaller distance across the domain. In an oilfield context, the continuous method predicts higher oil recovery (because more oil is displaced from the matrix) and later water breakthrough at production wells than the discontinuous method. Moreover, the discontinuous results show the fractures aligned with the water flow direction have contributed most to the saturation transport (Fig. 14). These fractures are dominant with respect to fluid transport for this particular test case. However, note that fluid transfer between fractures and matrix does occur despite the absence of capillary or gravity forces in these simulations, driven by the viscous pressure drop across the model domain and between upstream fractures and downstream matrix. Simulations on fracture networks with realistic connectivity and fracture apertures more typical of in situ stress conditions will further highlight the benefit of such discontinuous methods. The discontinuous method is presented here in the context of a static solid fractured porous skeleton through which fluids are transported. However, one future goal and motivation for developing this approach, as discussed below, is to capture dynamic effects which are known to be important in both the fluid and solid in the near wellbore region.

DISCUSSION
Research on fluid flow in fractures and in fractured porous media has a history that spans nearly four decades (Bear & Berkowitz 1987). Of the various conceptual models for flow and transport in fractured porous rock, there are three that dominate. One depicts the rock as a network of discrete fractures (a discrete fracture network or DFN), which can be applicable in cases where the matrix rock has no significant contribution to fluid storage or flow (e.g. Min et al. 2004). A second ('dual porosity' approach) assumes flow occurs primarily through a connected fracture network, which is modelled as a continuum with equivalent porosity and permeability, and accounts for fluid transport between stagnant matrix and flowing fractures using transfer functions that appear as source terms in the continuum description of fracture flow (e.g. Barenblatt et al., 1960;Warren and Root, 1963). An alternative continuum approach within this conceptual model is the continuous-time random-walk method (Berkovitz and Sher, 1998), which has proven to be very successful when calibrated using a discrete fracture and matrix approach (Geiger et al. 2010). A third approach (the discrete fracture and matrix or DFM method) models flow explicitly in both fractures and matrix, using a continuum representation with sharply varying material properties such as porosity and permeability to describe discrete fractures within the rock matrix (e.g. Kim & Deo 2000;Juanes et al. 2002;Matthai et al. 2007).
In most continuum models, the sharp material interfaces that are the characteristic of fractured porous media are difficult to represent unless they are assumed to have very simple geometries, such as simple planar surfaces, which are either parallel or orthogonal. Yet accurate simulation of flow through fractured rock requires models that can capture both the physics of the flow process and the flow geometry. Errors and inaccuracies often result from discretising the discontinua with continuum mechanics equations for the purpose of numerical simulation (Nick & Matthai 2011b). Numerical methods based on the finite-difference-finite-volume approach and the two-point flux approximation require a spatial discretisation that is socalled k-orthogonal (typically grids that are Cartesian or provide the perpendicular bisector (PEBI) property are used; see, for example, Wu & Parashkevov 2009) which greatly limits the complexity of fracture geometries that  can be represented. Methods based on the CVFE approach (and variants thereof) allow the use of flexible, non-korthogonal meshes to discretise space, which capture complex fracture geometries more accurately using fewer degrees of freedom (e.g. Geiger et al. 2004). However, even these approaches usually require very fine meshes to adequately represent sharp gradients in flow-related fields such as velocity or saturation that occur at material property boundaries such as fracture walls; moreover, CVFE approaches typically use control volumes that span material property boundaries which can lead to 'smearing' of CVbased fields across the boundary. Yet capturing flow in and around such discontinuities is of first-order significance if fluid flow in fractured porous media is to be represented accurately. Flow within the fracture, and in the matrix adjacent to the fracture, must be properly captured.
Owing to its importance, several recent numerical studies have addressed alternative methods to solve this 'discontinuous flow' problem for fractured porous media; that is, the problem of how to mathematically represent and solve the multiphase porous media flow equations in domains with sharp boundaries for saturation and velocity. A vertex-cen-tred finite-volume method for two-phase flow in fractured porous media was presented by Reichenberger et al. (2006), Epshteyn & Riviere (2007) and Jaffre et al. (2011), in which fractures are modelled as lower dimensional finite element/ finite volumes (lines in 2D models and surfaces in 3D models) with their own permeability and transmissibility that correspond to a prescribed fracture aperture. These methods for solving DFM models were illustrated for 2D problems by Latham et al. (2013) and Lei et al. (2014) who used a novel geomechanical simulation method to obtain the local variability of fracture apertures. They modelled the fracture network response to stress and produced a heterogeneously deformed solid skeleton; thus, the explicit local fracture aperture geometries were obtained from geomechanical simulations. However, to solve for the average flow properties for their one-way coupled, that is static fracture network, they too found it convenient to collapse the solid geometry and the fracture flow representation to lower dimensional elements, that is they modelled fracture flow by extracting the variable mechanical apertures and mapping the corresponding variable hydraulic apertures along the lower dimension median line of the fractures. Lower dimensional approaches may also be able to address solid deformation with full two-way coupling to capture highly transient phenomena. However, when fluid pressures are highly sensitive to the geometry of the solid skeleton, there could be major disadvantages for the lower dimensional approach and the unrealistic collapsing of the solid skeleton to form a solid that assumes parallel-sided zero-thickness fractures. To avoid errors in such an approach would require local compensatory changes in the stress field to account for unrealistic solid geometry and for such coupled transient behaviour, this is likely to be especially problematic.
For any formulation wishing to retain actual 3D topologies of voids or high aspect ratio fractures with sharp boundaries representing discontinuous saturations, these essentially continuous finite element approaches are computationally very challenging. Monteagudo & Firoozabadi (2007) extended the control volume discrete fracture method (CVDF) to incorporate heterogeneity in rock properties, and the new extended algorithms were employed to model flow with zero capillary pressure in the fracture. Hoteit & Firoozabadi (2008a,b) proposed an implicit, mixed finite element-finite-volume formulation with discontinuous Galerkin (DG). This method aims to capture discontinuities in saturation arising from capillary effects. Chen et al. (2008) illustrated a DG method in which saturation gradients were captured better than the standard Galerkin method. Nick & Matthai (2011a) developed a hybrid element discretisation with embedded discontinuities to simulate single-phase flow and transport through two-dimensional models of heterogeneous media, in which the effect of sharp material interfaces were taken into account. Eymard et al. (2012) presented the Vertex Approximate Gradient (VAG) scheme for discretisation of multiphase compositional Darcy flow.
The method we report here differs from these previous approaches by employing a 3D fracture representation. This allows the method to model two-way coupled transient behaviour, as it retains the explicit 3D representation of the solid and fracture wall, which is key to modelling shear reactivation or propagation of fractures. Clearly, there are likely to be computational speed advantages with lower dimensional models. However, valid comparisons of accuracy and CPU time for the two approaches (full and lower dimensional geometry) for both static and transient solid-fluid interaction fracture flow problems have yet to be performed to explore such differences.
The flow model proposed here offers the future possibility to couple with geomechanical models that track the solid velocities and stresses, including fracture wall deformation and fracture growth. Such coupling might occur as the solid medium responds to changes in fluid pressures and in boundary stress changes. This is the wider scope into which we introduce these methods; however, here we have reported the method in the context of fractured reservoirs where the solid skeleton is held static. The focus in this initial paper is on the novel discontinuous method developed to accurately and efficiently capture discontinuities in velocity, saturation and other flow-related fields across fracture walls and other such boundaries in material properties, as compared with the continuous methods presented previously and discussed above. The formulation includes gravity and capillary forces but these have yet to be fully tested, so here we have focused on viscous forces, which often dominate flow, especially at high rate and close to production or injection wells. Results including gravity and capillarity will be presented in future papers. The key, interrelated research objectives we have addressed are (1) the development of numerical methods that preserve sharp discontinuities in flow-based fields such as velocity and saturation at fracture walls and other boundaries in material properties in heterogeneous porous media; (2) the development of numerical methods that allow 3D representation of extreme aspect ratio features such as fractures, to facilitate coupled models of geomechanics and flow; and (3) the development of methods that achieve (1) and (2) in a computationally efficient manner. Moreover, we have presented results for a number of test cases that demonstrate the ability of the discontinuous CVFE method to accurately simulate multiphase flow through fractured porous media without artificially spreading saturation across the fracture walls. The first example illustrated the performance of different element pairs of structured and unstructured meshes of different resolution. Two-phase flow was then modelled in a range of idealised fracture networks. An issue raised by Nick & Matthai (2011b) is that discontinuous methods typically use more computational degrees of freedom than continuous methods and thus require longer run times when the same mesh is used. This is confirmed here (see Table 2); however, results suggest that much coarser meshes can be used with the discontinuous method to achieve similar or better levels of accuracy. This in turn impacts on the computation time. Overall, results here suggest that the discontinuous method is more efficient and computationally cheaper than the continuous one (CPU time ratio, 1:3.37). Further discussion of the performance and merits of continuous anisotropic adaptive mesh methods have been reported by Mostaghimi et al. (2015), who applied a similar mesh optimisation algorithm to that presented here but used a different approach to discretise and solve the governing equations.
In the examples shown here, the novel discontinuous formulation has been used to simulate multiphase flow through fractured porous media. However, the fracture patterns and apertures were idealised. They were set to be constant and were given constant permeability values, and this is not generally the case in hydro-mechanically coupled rock masses. In future work, it is envisaged that the solid skeleton and fracture network will be modelled with geomechanical methods (Latham et al. 2013;Lei et al. 2014). The solid skeleton, fracture walls and apertures, especially in the region of a wellbore, will be free to respond dynamically to the imposed effective stresses that may in turn be varying dynamically through solving coupled solid-fluid equations. In addition to modelling multiphase flows in fractured porous media in a solid skeleton assumed to be fixed, the methods presented here are likely to be particularly advantageous to the community addressing flows with dynamic fracture wall behaviour, such as in hydraulic fracturing problems.

CONCLUSIONS
(1) A novel implicit discontinuous CVFE method to model multiphase flow in fractured porous media based on the new P 2 DG-P 1 DG element pair has been described. This takes into account the jump in saturation at the interface of different materials. (2) Simulations with different mesh resolutions for a number of idealised fracture networks are undertaken to compare the new fully discontinuous formulation against conventional CVFE methods in which CVs span material property boundaries defined element-wise. (3) Results presented here show that, for a fixed mesh, the discontinuous method simulation run times are slightly longer than the continuous method for the same mesh and that the numerical solutions are stable and provide consistent results. (4) The small computational overhead is trivial given that much coarser meshes can be used with the discontinuous method to achieve similar accuracy with the continuous method and run times are significantly reduced.
mesh are contained within a single element. However, by