Flow in fractured porous media: A review of conceptual models and discretization approaches

The last decade has seen a strong increase of research into flows in fractured porous media, mainly related to subsurface processes, but also in materials science and biological applications. Connected fractures totally dominate flow-patterns, and their representation is therefore a critical part in model design. Due to the fracture's characteristics as approximately planar discontinuities with an extreme size to width ratio, they challenge standard macroscale mathematical and numerical modeling of flow based on averaging. Thus, over the last decades, various, and also fundamentally different, approaches have been developed. This paper reviews common conceptual models and discretization approaches for flow in fractured porous media, with an emphasis on the dominating effects the fractures have on flow processes. In this context, the paper discuss the tight connection between physical and mathematical modeling and simulation approaches. Extensions and research challenges related to transport, multi-phase flow and fluid-solid interaction are also commented on.


Introduction
Fractures, both natural and engineered, provide major conduits or barriers for fluid flows in many types of media. The existence of fractures that affect flow and transport are characteristic of different types of porous media for domains ranging from millimeters to hundreds of kilometers. Subsurface rocks are, for example, always fractured to some degree due to tectonic movement. Fractures are present in porous media like soil (Council 2001) and glaciers (Foun-tain and Walder 1998), and in porous materials such as wood (Smith et al. 2003;Watanabe et al. 1998). Several anthropogenic materials, like concrete (Carmeliet et al. 2004;Rouchier et al. 2012), have similar characteristics, intendedly constructed or not. Thin membranes or other thin layers in porous materials, for example in fuel cells Hassanizadeh 2014, 2015) or biological materials, can conceptually also be modeled in the same manner.
For a porous medium, secondary permeability provided by conductive fractures can result in an effective permeability several magnitudes larger than the primary permeability provided by connected pores. Moreover, there is generally no separation of length scales in fracture sizes; thus, upscaling into average properties often results in poor model quality. As a supplement, a wide range of methodologies tailored for fractured media have been proposed in recent years, aiming to accurately capture the effect of complex fracture networks on flow patterns. Aided by increases in computational power that make higher resolution models feasible, developments typically focus on one or more of the following aspects: (1) accurate representation of flow in the fractures, (2) accurate representation of flow in the interaction between fractures and the surrounding porous media, and (3) ability to handle complex fracture network structures.
This paper reviews the most important modeling approaches for flow in fractured porous media, from physical, conceptual and mathematical models, to discretization approaches. It provides an integrated view on modeling and discretization and therefore does not give details on specific numerical schemes but rather present broader classes of methods, covering strengths and weaknesses in terms of modeling capabilities, computational accuracy and cost. Approaches to be covered include methods that are designed to work with traditional simulators for flow and transport as well as models that are based on explicit representation of fractures. Discretization approaches have different characteristics and features, depending on how fractures are represented in the computational domain, and how flow between fractures and between fractures and surrounding porous media is approximated.
Researchers new to the topic of modeling and simulation of flow in fractured porous media are provided with a broad introduction, with emphasis on dealing with challenges arising from the strong effects networks of fractures have on aggregated flow patterns. The general overview of physical, mathematical and conceptual models and discretization approaches aim to contribute in bridging research developments across different fields and application areas. The paper is not intended as a comprehensive review of the topic, and for the sake of brevity, the references are not exhaustive. More detailed discussion of some of the aspects addressed herein, and related topics can be found in previous books (Adler et al. 2012;Dietrich 2005;Sahimi 2011;Faybishenko et al. 2015;Bear et al. 2012) and review articles (Berkowitz 2002;Neuman 2005).
The paper starts by discussing the representation of fractures and fracture networks in porous media in Sect. 2, while models for flow are presented in Sect. 3. In Sect. 4, conceptual models for representing porous media with complex networks of fractures are reviewed. Discretization approaches are discussed in Sect. 5. In Sect. 6, the choices between different modeling approaches are discussed as well as research challenges that remain, before the paper concludes in Sect. 7.

Geometric Representation of Fractured Porous Media
Fractures in a porous medium are discontinuities in the medium in the form of narrow zones that are locally approximately planar and have distinctly different characteristics than the medium itself.
The detailed geometry of a fracture is represented by its boundaries to the medium on each of its sides. While fractures are generally thin compared to their extension, they are much wider than a typical pore diameter, and their length can span orders of magnitude, with a largest extent that can be comparable to the size of the domain of interest. At the same time, the total volume corresponding to the fractures is small compared to the volume of the surrounding medium. Due to these characteristics, fractures in a porous medium are typically approximated by two-dimensional inclusions in macroscale modeling of flow, as shown in Fig. 1a.
To the utmost, a fracture can be empty or it can contain an impermeable filling material. As a consequence, fractures can be both highly conducting or highly sealing compared to the host material. Networks of multiple, connected fractures structurally dominate flow patterns based on network geometries that can be highly complex. A simulation using the software PorePy , shown in Fig. 1c, illustrates this.
Modeling flow in porous media with complex networks of fractures calls for adequate representations. A major challenge is the difficulty of defining and measuring model parameters that allow for continuum-scale modeling, both in situations where fractures are treated explicitly and when they are integrated as part of the porous medium itself. In modeling of standard porous media, parameters used in macroscale modeling are related to pore-scale quantities by averaging over domains of increasing size. If, with increasing size of the averaging domain, the averaged quantity becomes approximately constant, the corresponding domain is denoted a representative elementary volume (REV) (Bear 1972). Parameters defined by the REV concept allow for building of flow models with quantities that vary continuously in space. Compared to the standard descriptions of a porous medium, fractures introduce intermediate length scales, but there may be no clear scale separation between the pore scale, fracture widths, fracture lengths and the macroscale of interest. Hence, the presence of fractures challenges the existence of an REV in porous media modeling. The resulting difficulty in defining averaged quantities is a recurring topic in this paper. In the following, a brief discussion of the representation of individual fractures and fracture networks in a fractured porous medium is given. Fig. 1 a A 2D conceptual drawing of a fracture network, including different intersections and non-planar fractures; b notation for fractures, intersections and interfaces; c simulation of heat transfer by convection, dominated by fracture flow 123

Individual Fractures
A fracture is bounded by two surfaces, which may be rough and locally oscillating. A thorough description of the geometrical and statistical modeling of a single fracture is given by Adler et al. (2012). Herein, the focus is not on the representation of the detailed geometry of a single fracture but rather a representation of the fracture at a scale much larger than the distance between the fracture's bounding surfaces, which is suitable for modeling of flow through the fractured medium at the network scale.
If a two-dimensional coordinate system is defined parallel to the plane defined by the fracture, the pointwise distance between the two bounding surfaces of the fracture defines its local aperture. If the fracture's bounding surfaces are locally fluctuating, the local aperture varies significantly over short distances. By representing the fracture's bounding surfaces by average planes, the average aperture, e V , of the fracture can be defined as the distance between these two planes. This gives an idealized measure of fracture width that does not account for the possible fine-scale fluctuations of the fracture's bounding surfaces.
In modeling of different physical processes, such as fluid flow and mechanical deformation of a fracture, different concepts for defining an idealized parameter that represents the width of the fracture in the context of the given process are commonly employed. This results in definitions of apertures that are conceptually and numerically different from the average aperture defined above. An example is given in Sect. 3, which introduces the hydraulic aperture of a fracture.

Fracture Networks
In networks of fractures, the size, orientation, location and aperture that characterize individual fractures can be complemented by measures of fracture density. For fractured rocks, both properties of individual fractures and networks can be characterized by various statistical distributions (lognormal, gamma, exponential, power law) and also by fractal descriptions (Bonnet et al. 2001;Liu et al. 2016). An argument for power law and fractal scaling is the lack of characteristic length scales in fracture growth. In the case fractures exist on multiple scales, even if the most dominating fractures may be represented explicitly, as discussed in Sect. 4, there may not exist an REV for the remaining part of the domain.
The connectivity of a fracture network is a central defining property considering flow and transport processes and can be studied by percolation theory (Robinson 1983;Berkowitz 1995;Thovert et al. 2017). When a network is statistically homogeneous with well-defined parameters, percolation thresholds can be studied and predicted. Close to percolation thresholds, connectivity is highly sensitive to fracture density. For general networks of fractures that are not well described statistically, individual members may be decisive for the connectivity of the network.

Modeling of Flow in Individual Fractures and the Porous Matrix
Fractured media provide particular challenges to modeling and simulation that are not present in standard porous media. These can be explained by the interaction between structural properties of fractures and fracture networks with dynamic processes that takes place in the domain. The fractures do, for example, strongly affect the nature of fluid velocities. Leading processes may not be well described in a continuum manner, as the fractures by definition introduce strong discontinuities that are not well represented by averaged descriptions (Long et al. 1982). Moreover, flow processes can also affect the nature of the fracture networks, for example due to mechanical or chemical fluid-solid interactions.
In this section, the basic equations for modeling of flow in fractured porous media are presented, considering a fractured porous medium domain . The modeling is based on a domain-decomposition approach, considering the different structural features of the fractured medium. Following standard conventions, the continuum host medium is referred to as the matrix, noting the ambiguity of the use of the term matrix, as this term is also used for the solid skeleton when describing the microscale structure of a porous medium. The matrix, M , can be porous, and it can also possibly integrate the effect of fine-scale fractures. The part of the domain consisting of explicitly represented fractures is divided into fracture intersections, I , and the remaining part of the fractured domain denoted F . The interface between F and M is denoted M F , while the interface between I and F is denoted F I . A conceptual drawing of a fracture network is provided in Fig. 1a, while the division of the fractured part of the domain and defined interfaces are illustrated in Fig. 1b. We assume that flow only takes place between objects one dimension apart, that is, between M and F and between F and I but not directly between M and I .

Matrix Flow
Assuming that there exists an REV for the matrix part of the domain, single-phase flow at the associated continuum scale is governed by mass conservation and Darcy's law. The mass conservation equation is given by where ρ is fluid density, v is the flux, and q represents a source term. In M , Darcy's law relates the flux v to the gradient of the pressure p by where K is permeability, μ is viscosity, and g is the gravitational acceleration.

Flow in an Individual Fracture
A fracture in a porous medium can be empty, or contain filling material, which makes the fracture itself a porous medium. If a fracture is treated as a thin but fully three-dimensional inclusion with a 3 × 3 permeability tensor, flow is governed by Eq.
(2) with a discontinuity in permeability associated with the fracture. However, models that exploit the representation of a fracture as a lower-dimensional surface are more feasible when considering multiple fractures that span a complex network. An idealized situation is the case of laminar flow in a fracture residing in an impermeable matrix. By averaging of the Stokes equations in the case of an empty fracture (Adler et al. 2012), or of Darcy's law in the case the fracture is filled with a porous material (Martin et al. 2005;Angot et al. 2009), the average specific discharge (flux) is where K is a 2 × 2 tensor and g is the component of g in the fracture plane. The pressure, p, should be interpreted as a macroscale pressure, and ∇ is the in-plane gradient operator. Due to the resemblance of Eq. (3) with Darcy's law (2), K is denoted the fracture permeability. The form of K depends on the configuration of the fracture's bounding surfaces and on possible filling material. An idealized representation of a conductive fracture is an open channel between two parallel planes separated by a constant aperture e V = e h . With further assumptions of incompressible flow and no-slip conditions on the confining planes, the fracture permeability is isotropic, K = κ f I, and is given by ; (4) see, e.g., Adler et al. (2012) for details. Equation (4) illustrates the strongly nonlinear relationship between fracture aperture and flow. Considering flow across a cross-section of a fracture, the flow rate through a fracture per unit length is given by e V v F ∝ e V K . That is, in the case of an empty fracture between parallel plates, Eq. (4) results in a cubic relation between aperture and flow, and, hence, Eq. (4) is often denoted the cubic law.
To accommodate realistic fractures, the fracture permeability must account for the roughness of fracture walls and the presence of gauge or other filling material inside the fractures. A result is that the fracture's average aperture, e V , may be much larger than its hydraulic aperture, e h , defined by Eq. (4). Still, there is in general a nonlinear dependence of the permeability on the aperture, and an important consideration in the construction of a simulation model is thus how accurately the aperture, including its variation along the fracture, and the associated fracture permeability, is represented.
To estimate the hydraulic aperture of a fracture is a non-trivial task in terms of both modeling and experiments (Adler et al. 2012;Berkowitz 2002). Furthermore, for open fractures having anisotropic geometrical characteristics, the hydraulic aperture must be described by a tensor permeability K defining preferential directions for flow. This is often the case for rock fractures due to shearing (Auradou et al. 2005), or for filled fractures with anisotropy in the filling material.

Flow Across Fracture-Matrix Interfaces
In the general situation, the matrix is permeable, and there can be flow between fracture and matrix. The mass conservation equation is then the same as before for the fracture but with additional source terms representing flow into the fracture from the matrix. Assuming for simplicity that e V is constant in time and uniform in space, the equation becomes where v ⊥ M F denotes the flux across the interface M F , and the superscript on M F denotes the restriction of M F to the different sides of the fracture.
As the fracture is represented as a thin planar layer, flow in the direction tangential to the fracture is governed by Eq. (3), but an additional equation is needed to account for flow between fracture and matrix at each side of the fracture; that is, v ⊥ . Following Martin et al. (2005), by averaging, fluxes across the fracture-matrix interfaces at each side of the fracture are approximated by where n ± denotes the normal to ± M F pointing into M . Furthermore, κ ⊥ is the permeability in the normal direction, which represents conductance to flow across the fracture, p ± M is the pressure in the matrix at ± M F , and p F is the fracture pressure. The factor 2 comes from taking the pressure difference over half the fracture aperture. For the matrix domain, the fluxes (6) appear as Neumann-type boundary conditions; that is, The matrix pressure will thus experience a jump across the fracture surface, with a magnitude that depends strongly on whether the fracture is conductive or sealing.

Flow in Fracture-Fracture Intersections
The adaption of Eq. (3)

for flow in an intersection between two fractures is
where v | F is the flux along the one-dimensional interface, κ | is the permeability of the intersection, ∇ | is the one-dimensional gradient operator along the intersection, and g | is the component of the gravity acceleration vector along the intersection. For flow between a fracture and a fracture intersection, Eq. (6) can be adapted straightforwardly, considering the difference between the pressure at the boundary of the intersection and the pressure in the intersection ).

Modeling of Fractures as Lower-Dimensional Objects
Given the presented model, fractures can be treated as co-dimension one domains in the medium based on a domain-decomposition approach for intersections, fractures and matrix (Alboin et al. 2002;Angot et al. 2009;Martin et al. 2005). That is, for a three-dimensional domain, fractures can be considered as two-dimensional surfaces. Similarly, fracture intersections can be treated as co-dimension two domains (lines) in three dimensions ). This mixed-dimensional modeling of domains is useful in modeling as it avoids thin equi-dimensional domains representing the fracture.
In a co-dimension one representation of a fracture, its aperture can be treated as a parameter that possibly can vary, without the need to update the matrix domain. This is advantageous in coupling flow with dynamical chemical or mechanical processes that change the fracture aperture (Ucar et al. 2017).

Extensions to More Complex Flow Physics
The linear relation between flow and pressure gradient in Eq. (3) relies on the assumption of laminar flow (Whitaker 1986). If this assumption does not hold, because of high imposed pressure gradients relative to the resistance from fracture walls and gauge material, the result is reduced flow rates compared to the linear relation in Eq. (3). The standard extended model involves a Forchheimer correction term, which introduces a nonlinear coupling between pressure gradient and flow rates. Considering fractures as co-dimension one objects, the domain-decomposition framework presented above can be extended to include the Forchheimer correction in the fracture part of the domain (Frih et al. 2008). Furthermore, considering the development of boundary layers developing at the boundary between the matrix and fracture domains, the effective fracture permeability is shown to depend on the ratio of bulk and interface permeabilities, the fluid viscosity, and the fracture aperture (Vernerey 2012).

Conceptual Models for Representation of Fractured Porous Media
The dynamics of flow within a single fracture and the interaction with the surrounding matrix is in itself complex but feasible to treat with the models in the previous section. Additional challenges arise when these dynamics take place in a network of fractures, where the fractures together have decisive impact on the dynamics in the simulation domain. The fundamental choice in modeling and simulation of macroscale processes in fractured porous media is whether to represent fractures explicitly, or represent them implicitly by an effective continuum. This is not a binary choice, and several approaches are available that employ a combination of the two. In the following, different conceptual models, as illustrated in Fig. 2, are described.

Single-Continuum Models
In single-continuum models, the fractures are represented by adapting the permeability of the porous medium. The permeability may increase or decrease, and the orientation of the tensor may change depending on the properties of the fractures and the fracture network. Ultimately, one aims to identify an effective permeability K e of the fracture network together with the porous matrix such that where the angular brackets · denote the average over a suitable volume that contains fractures and porous matrix. In words, one expects that also at the larger scale the mean value of the velocity field is proportional to the mean value of the gradient of a potential. At this level, the approach is identical to classical upscaling of heterogeneous porous media. The averaging volume should ideally be corresponding to an REV, which as discussed may not exist due to the lack of scale separation in the fracture network. An effective permeability can still be computed for a given averaging volume, say a grid block, and used in simulations, but as a modeling concept this approach is less powerful than when an REV is available. To emphasize the difference, the term equivalent (grid-block) permeability (Durlofsky 1991) is used in the context of discretizations. An upscaling technique should account for inherent properties of the fractures themselves such as orientation, aperture, area, surface roughness and filling material as well as properties of the fracture ensemble such as fracture density and coordination number. Classical techniques for permeability upscaling, including methods based on numerical upscaling, asymptotic expansions and effective medium techniques, can all be applied to fractured media. Challenges pertinent to those techniques, as discussed by, e.g., Farmer (2002) or Fig. 2 Conceptual models of a fractured porous medium. The illustration provides a generalized picture of trade-offs in choosing between different models for fractured porous media and how they might relate considering a specific model problem for a given fractured medium. The placements of the different models relative to each other depend on the original fractured porous medium, to which extent scale separation exists, and on the applied modeling choices and upscaling procedures as well as the available information on the original medium. For example, if a DFM model represents all fractures of the original medium, it is a perfect model with the same placement as the original medium in the illustration. The same is true for the single-continuum model in the special situation where an REV exists Gerritsen and Durlofsky (2005), are amplified in fractured porous media. Furthermore, the merits of different upscaling approaches for a fractured medium differ extensively based on the amount and type of information that is required about the medium and the procedures involved. Oda (1985) devised an effective medium theory-based method that calculates the effective permeability of the fracture network based on fracture shape, size, aperture and orientation distributions of the fractures as well as matrix permeability. It does not account for the network properties, but it is still very popular due to its simplicity. A recent review on explicit expressions including experimental observations as well is provided by Liu et al. (2016). Even though it does not account for flow in the matrix, percolation theory has played a major role (Sahimi 1993). Also, the fractal nature of fracture distributions has been used to develop upscaling techniques (Liu et al. 2016). Further, effective medium theory has been adapted to cope with challenges posed by individual fractures (Saevik et al. 2013).
The single-continuum representation is desirable because of its low level of computational and model complexity. In comparison with non-fractured porous media models, no additional terms are required and the transfer of methods and simulation technology is straightforward. This comes at the cost of simplification. For single-phase steady-state flow, it is always possible to calculate an equivalent permeability for a given boundary value problem such that the average flux follows the upscaled Darcy equation (8). This equivalent permeability is, however, length scale, boundary condition and process specific. For more complex problems, even as simple as the evolution of pressure in compressible flow or tracer transport, the dynamics are not adequately represented by averaged flow alone. Thus, the boundary value problem to determine the equivalent permeability should be representative for the real flow field, which is not known a priori.

Multi-continuum Models
Multi-continuum approaches are families of methods that represent the fractured porous medium by several superimposed media with their own conservation equations and constitutive laws. The mass conservation equation for continuum α reads as where T α denotes the transfer into the continuum α from all other continua. The flux v represents the upscaled flux in continuum α, which is almost exclusively described by Darcy's law (2) with according upscaled permeabilities. The simplest approach consists of a fracture continuum α = f and a matrix continuum α = m and is referred to as dual-continuum model; see Fig. 2 for an illustration. An extension to a general number of continua can account for heterogeneities in properties and flow rates within the fracture continuum.
The underlying assumption of the multi-continuum approach is that it is possible to represent different domains in the pore-space and the dynamics within as well as the exchange between those domains based on continuous parameters and variables. Hence, the discussions and challenges for the single-continuum representation in Sect. 4.1.1 with respect to the existence of REVs apply to multi-continuum formulations as well. Due to the increased complexity of the equations, the issue is even more involved in the multi-continuum setting, and derivations of the models from a fine scale are available only under very specific assumptions on geometries and flow (Arbogast et al. 1990). More generally, multi-continuum models should be considered constitutive formulations that can still be useful, in particular when some continua provide volume while others provide permeability. Note that sealing fractures provide neither, and, hence, are usually not modeled through multi-continuum approaches.
The transfer term T α is the art and the challenge of multi-continuum modeling. It should represent the aggregated transfer behavior based on variables of the continuum formulation. In order to be predictive, its parameters should be derived from geometrical, topological and physical properties of the individual continua. The only physical constraint is mass conservation. For dual-continuum formulations, this implies that T m = −T f holds true. The first explicit form for the transfer term was proposed by Warren and Root (1963) and Barenblatt et al. (1960) as where the shape factor σ accounts for geometrical properties such as effective fracturematrix interface area and effective characteristic length over which the transfer takes place, and β accounts for physical properties such as matrix and fracture permeability and fluid viscosity. For simple regular and symmetric structures of the fractured porous medium and under additional constraints on the magnitude and heterogeneity of parameters and variables, Eq. (10) accurately represents transfer between fracture and matrix, with σ and β taken as constants. For some of those configurations, explicit formulas to determine σ and β have been developed. However, in many applications pressure differences are substantial and the fracture geometry is asymmetric and complex and, hence, the representation of the transfer through Eq. (10) is inaccurate. The common choice in extending the validity of Eq. (10) is to keep the linear pressure difference, and generalize the coefficient in front. No general approach is known; instead, the available expressions are strongly tied to assumptions on the flow regime and fracture network configuration. For processes with negligible buoyancy, linear density-pressure relationships and vanishing macroscopic flux in the matrix, upscaling techniques such as homogenization through asymptotic expansion (Arbogast et al. 1990) have been used to determine the transfer term. These yield transfer terms with memory functions that are non-local in time and therefore notoriously challenging to interpret or exploit in solutions and simulators. To simplify calculations, these have been expanded in series of exponentials yielding a series of first-order transfer terms with varying transfer rate coefficients (Tecklenburg et al. 2016). Also heterogeneous fracture spacings yield series of first-order transfer terms with varying coefficients due to varying shape factors (Haggerty and Gorelick 1995). Both are referred to as multi-rate transfer models. Another way of coping with the shortcomings of the transfer term in Eq. (10) is to expand the matrix continuum into a series of multiple interacting continua (MINC) with vanishing macroscopic flow (Pruess 1985).
For monotonic transfer processes, analytical solutions for the diffusion equations are available that produce explicitly time dependent shape factors (Rangel-German and Kovscek 2006) or explicit parameterizations of the transfer with respect to time (Lim and Aziz 1995;Zhou et al. 2017). While these approaches improve the accuracy, their implementation in simulators tends to be inelegant.

Explicit Representation of Fractures
Continuum models make no geometric distinction between M and F but instead employ different media that are coexisting in space. In contrast, conforming methods preserve a notion of M , F and M F as separate geometric objects. Depending on the specific approach being employed, this reduces or removes altogether the need to represent flow by upscaled quantities. As such, modeling within a framework with explicit representation of fractures is often conceptually simpler than the implicit counterpart, but at the cost of dealing with complex geometries.

Discrete Fracture Matrix Models
Discrete fracture matrix (DFM) models attempt to strike a balance between loss of accuracy by upscaling and geometric complexity. The fluid is located in explicitly represented fractures, defined as F as well as in the porous matrix in between, defining M . This allows for integrating fractures of a length much smaller than the size of the domain into M , providing secondary permeability, rather than forming part of the explicitly characterized fracture network. The DFM model is appropriate for handling fractures in a permeable medium that should be explicitly represented due to their structural and dominating impact on processes.
In principle, all fractures within the domain can be incorporated into F , and the governing equations are then given in Sect. 3. In practice, it is not feasible to retain an explicit representation of all fractures, mainly for reasons of computational complexity, as discussed in Sect. 5. DFM models therefore preserve some fractures in F , while others are upscaled, and replaced by averaged quantities within M . Hence, the effective equations will be those of a fractured porous medium, even if the host medium is considered impermeable. By resolving the interface between fractures and the rock, DFM models facilitate detailed modeling of processes on this interface, including capillary pressures in two-phase flow, and hydro-mechanical interaction between flow and the host medium.
Methods to decide which fractures to assign to F and M , respectively, are not well developed; common practice is to base the selection on fracture length (Lee et al. 2001). More advanced, and application dependent, selection criteria could include factors such as fracture connectivity and the volume of host material that has a fracture as its main connection to the global flow field could be included. The type of upscaled representation can also be varied: The traditional approach has been a single-continuum method (Karimi-Fard et al. 2004), but multi-continuum approaches can also be applied (Sandve et al. 2013;Jiang and Younis 2017).

Discrete Fracture Network Models
The network of fractures, disregarding the host medium, is commonly referred to as a discrete fracture network (DFN). By analogy, models that ignore flow in M , or, equivalently, consider M impermeable, are referred to as DFN models. In a DFN model, all fluid is assumed to be contained within the fracture network, which is represented by a lattice. As the material in between the fractures is considered impermeable, this is strictly speaking not a model of a fractured porous medium, but the model is still included herein for completeness. In particular, DFN flow models are appropriate for porous media where the entire porosity and permeability are due to fractures that can be explicitly represented. The model is commonly also used as a generalization to model fractures in low-permeable porous media.
Flow and transport are modeled in individual fractures i F based on Eq.
(3), with network interaction via the intersections I , as discussed in Sect. 3. This gives an accurate representation of dynamics in fractures, which can be of use in multi-physics couplings, such as reactive flows and mechanical deformation of fractures.

Numerical Modeling
The conceptual models presented in the previous section can all be used as starting points for simulation models. This section provides an overview of prevailing techniques and challenges in numerical modeling. The main topics are the preservation of heterogeneities between fracture and host medium, and the discretization of the interaction between the two. The full simulation model also requires discretization of the flow equations within the fractures and, in most cases, in the matrix. This is commented on when appropriate. In general, it is referred to Lie (2016) for an introduction to discretization methods for porous media flow. Quantitative comparisons between a number of different numerical approaches are presented for DFM models by Flemisch et al. (2018), DFN models by Fumagalli et al. (2018) and DFM and dual-continuum models by Moinfar et al. (2013).
Closely related to the conceptual models discussed in Sect. 4, which treat fractures implicitly (single-and multi-continuum models) or explicitly (DFN and DFM models), the numerical approaches are divided into three groups: (i) models with implicit fracture representation, (ii) models with explicit fracture representation with non-conforming meshes to

Implicit Fracture Representation
When fractures are represented implicitly, the PDE model contains no direct information of fracture geometry and is, hence, discretized with meshes that are non-conforming to the fractures.

Single-Continuum Approaches
Single-continuum models are by far the simplest approaches with respect to discretization: At any point, in space there is a single value for each parameter such as permeability and porosity. This leaves us with a standard simulation model for porous media flow.

Multi-continuum Approaches
Depending on which type of multi-continuum approach taken, these models have several coexisting values for physical parameters, each assigned to a separate medium. As for the single-continuum models, discretization of the individual media can be treated as standard (non-fractured) porous media flow, applying any number of numerical techniques. The main difference lies in the treatment of the interaction term: For linear approximations, this will not significantly alter the characteristics of the system, although the stiffness of the problem is increased and may cause convergence issues for linear and nonlinear solvers. More advanced interaction terms (nonlinear, non-local or time dependent) may substantially increase the complexity of the problem. For each included medium, the number of degrees of freedom increases significantly (the number of physical cells for single-phase flow, twice the number for two-component flow etc.). This can substantially increase the computational cost of the simulation, in particular due to increased time spent by linear and nonlinear solvers.

Explicit Fracture Representation with Non-conforming Meshes
The above approaches have in common that upscaling removes the distinction between individual fractures prior to the introduction of a numerical mesh. If the conceptual model takes a conforming approach to representing the fractures, these must either be explicitly resolved by the mesh, as discussed in Sect. 5.3, or upscaled as part of the discretization, as discussed here. The resulting simulation models typically have separate discrete variables for fractures and matrix. Difficulties that need to be addressed are associated with (i) the amount of geometric complexity contained within individual cells, (ii) discontinuities in matrix pressure across the fracture, caused by resistance to flow according to Eq. (6).

Embedded Discrete Fracture Methods
Embedded discrete fracture methods (EDFM), introduced by Li and Lee (2008), are based on the DFM conceptual model, and they are discretized with own degrees of freedom for M and F . In the separate discretization of flow within M and F , the latter commonly borrows techniques for conforming discretizations of DFM models discussed in Sect. 5.3.2; see, e.g., Moinfar et al. (2014). Flow between fracture and matrix within a matrix cell is taken as proportional to the pressure difference between fracture and matrix, with a proportionality constant that reflects the geometry of the fracture. The coupling structure of the discretized system is identical to that of dual-continuum models, as shown in Fig. 3, with the difference that the coupling terms are modeled in terms of discrete variables directly.

Extended Finite Element Methods
While the methods discussed so far can be implemented using several discretization schemes, the extended finite element method (XFEM) is tightly bound to finite element methods. Originally developed for fracture mechanics, the method has been applied to flow problems by, e.g., Fumagalli and Scotti (2013), Schwenck et al. (2015) and Flemisch et al. (2016). On non-conforming meshes, finite element approximations of the pressure are poor in cells that contain fractures, as the finite element basis functions cannot capture the pressure discontinuity over the fracture surface. XFEM enriches the finite element space by adding basis functions designed to capture the discontinuity. While the idea is simple, the number and design of the basis functions depends on the geometric configuration of the fractures within the cells, and implementation is difficult for fracture crossings, endings within a cell, etc. (Schwenck et al. 2015). These complications increase further in 3D, and to date, XFEM for flow problems has been used almost exclusively in 2D.

Non-conforming Meshing for DFN Models
A recent approach, developed by Berrone et al. (2013), applies non-conforming meshing to fracture intersections, mainly for DFN models. Flow in individual fractures is discretized independently, with fracture intersections treated as sink or source terms inside cells. Interaction between fractures is posed as a minimization of pressure and flux discontinuities over the intersections. This approach has no problems in dealing with extremely complex configurations of intersection lines and has been applied to highly complex fracture networks.

Explicit Fracture Representation: Conforming Mesh
Conforming methods avoid the computation of transfer functions between M and F by explicitly meshing F and I . This gives high resolution of near-fracture dynamics, which can be especially useful for multi-physics couplings. The price to be paid is in mesh construction, which can be extremely difficult, and computational cost, as the simulation model in general contains many more cells.

Meshing
A conforming mesh must by definition explicitly resolve M F by treating fractures as constraints in the mesh construction. By extension, the mesh should also conform to I . Unless the fracture geometry is highly regular, fitting a structured mesh with a Cartesian topology to the fracture surfaces is virtually impossible, and unstructured meshing algorithms mush be applied. Most algorithms apply simplex meshes (triangles in 2D, tetrahedrals in 3D), but Voronoi meshes have also been considered (Merland et al. 2014) To understand why simulations on conforming meshes can be challenging, recall that the performance of numerical methods in general is tied to the shape and number of cells: Most methods perform poorly on cells with small angles, often yielding low accuracy solutions for extreme cases. Moreover, computational cost in the best case scale linearly with the number of cells, and the effectiveness of linear and nonlinear solvers often depends on the ratio of the smallest to the largest cell, depending on whether effective preconditioners are available. In addition, the smallest cell size is bounded above by the size of the smallest fracture to be explicitly represented; in practice, the cells may be much smaller.
In the case of an equi-dimensional representation of the fractures, a very fine mesh is needed to resolve the fractures and avoid large aspect ratios of the fracture cells. With the mixed-dimensional representation of the fractured media introduced in Sect. 3.5, the dimension of the fracture cells is one lower than that of the matrix cells, and this issue is avoided. Nevertheless, meshing remains challenging for complex network geometries: Fig. 4 shows four different configurations known to cause problems in 2D. All of these configurations will either require highly refined meshes or cells with high aspect ratios; the fractures in Fig. 4b will also necessarily require cells with sharp angles. The complications are even more severe in 3D.
Meshing algorithms used in simulations can either be general purpose, drawing upon sophisticated software packages (Geuzaine and Remacle 2009;Si 2015;Shewchuk 1996), or dedicated to fractured porous media (Holm et al. 2006). Interestingly, the latter class includes approaches that allow slight modifications in the fracture geometry to ease the mesh construction (Karimi-Fard and Durlofsky 2016;Mallison et al. 2010), or even replacement of troublesome fractures in a stochastic setting (Hyman et al. 2014). Even when using such adaptations, meshing remains a bottleneck for conforming discretizations due to the inherent geometric complexity of fracture networks.

Discretization Approaches
Discretization of conformal models can be categorized as equi-dimensional, mixeddimensional or hybrid approaches. The model for the equi-dimensional approach is essentially that of a standard porous medium, with heterogeneous parameters in M and F , and, hence, this is not discussed further. Discretization of mixed-dimensional conforming methods is naturally treated as a domain-decomposition method, wherein the flow equations are discretized on M , F and I separately, leaving out M in DFN models, and then coupled by interface conditions on M F and F I . For the flow problem, this approach can be shown to be stable under very mild assumptions on the discretizations (Martin et al. 2005;Angot et al. 2009;Boon et al. 2018;Nordbotten et al. 2018). The hybrid approach to conforming discretizations is based on equi-dimensional modeling. However, fractures are treated as lower-dimensional objects (surfaces) in mesh construction, avoiding difficulties of meshing close parallel surfaces. In the discretization of flow equations, faces on the fracture surfaces are then expanded into equi-dimensional cells by a virtual extension perpendicular to the fracture surface. This is a convenient way to introduce a fracture model in an existing simulation model (Karimi-Fard et al. 2004).
The computational cost of a conforming method can in part be traced to cells in I , which typically have much smaller volume than cells in M and F . This results in badly conditioned linear systems in the flow problem, and time stepping constraints for explicit transport schemes. This can be circumvented by eliminating variables in I (Karimi-Fard et al. 2004;Sandve et al. 2012), with a varying loss of accuracy .
previous sections handle a range of scenarios based on the underlying structure of the fracture network.
In this section, aspects that may enter in selecting a specific model are discussed, before extensions to multi-phase flow and more general multi-physics couplings are considered.

Characteristics of Different Modeling and Discretization Approaches
A main remaining challenge in modeling and discretization of flow in fractured porous media is the appropriate selection of the conceptual model (as illustrated in Fig. 2) and corresponding governing equations as well as the necessary upscaling procedures and discretization approaches. This selection is closely linked to the specific modeling problem. In some cases, a continuum model, where the fractures are integrated in an effective porous medium, is sufficient. This is primarily the case when fractures are weakly connected and small compared to the characteristic length scales of interest. If fractures occur as dense connected networks, multi-continuum models can capture the different characteristic fluid velocities in the fractures and the surrounding domains. When the local geometry of the fracture network affects directly the flow patterns at the scale of interest, fractures must be represented directly, by DFM or DFN models. For these models, a lower-dimensional representation of fractures is beneficial from a modeling and, in particular, a meshing point of view. A combination of models, by restricting different models to different parts of a domain, or by using models with explicit representation of fractures as basis for upscaled continuum models, is also an alternative. This implies a balancing between discretization errors, directly associated with the mesh size and quality, and modeling errors caused by upscaling fractures. Little is known of how these errors relate. Furthermore, modeling choices are related to application-specific availability and quality of data as well as modeling purposes.
When aiming for simpler effective representations of fractured porous media, one has to account for the challenge that fractures introduce new scales compared to the standard microscale (pore) and macroscale (continuum scale) models of porous media. This results in a lack of scale separation, and a proper REV may not exist. Depending on how pronounced these features are, fractured porous media are more or less suited for upscaled representations. However, constraints on resources as well as available information often trigger a need for at least a partial effective representation of dynamics, often relying on engineering and physical intuition to balance errors and costs.
In general, the different options for conceptual model and discretization move the complexity, and the error, between modeling and discretization. This is best illustrated by comparing a non-conforming approach, be it multi-continuum or EDFM, and a conforming approach. From a discretization standpoint, non-conforming methods do not require complex meshes, and simulations are fast and simple compared to conforming methods. The price to be paid is in complex modeling and low resolution of fracture-matrix interaction. Fundamentally, these issues are caused by the inherent complexity of porous media flow dominated by fracture networks, and modeling of the flow will remain difficult.

Extensions to Multi-phase Flow
The focus of the previous sections has been the modeling of single-phase flow in fractured media. However, in many applications two or more fluid phases share the pore-space. The interplay of capillary forces, pore geometry, buoyancy and viscous forces then leads to nonlinearly coupled flow. Despite more than a century of research, the modeling of these phenomena across different length scales still forms an active research field and many challenges translate directly to fractured media. Providing a brief discussion on challenges and activities, the following discussion is grouped into three themes: (i) the conceptual and constitutive representation of multi-phase flow phenomena in a single fracture, (ii) the mathematical and numerical representation of conventional multi-phase Darcy terms in explicit representations of fractures and (iii) the representation of multi-phase flow phenomena in implicit representations of fractures.
The availability of experimental or numerical data for the constitutive relationships, capillary pressure and relative permeabilities, is surprisingly sparse compared to non-fractured media. However, since the first experimental works (Pruess and Tsang 1990;Persoff and Pruess 1995) and modeling work (Rossen and Kumar 1992), significant progress has been achieved (Bogdanov et al. 2003;Chen and Horne 2006). Characterization of flow patterns in a fracture, and the impact they have on constitutive functions and their upscaling, is currently an active field of research (Karpyn et al. 2007;Babadagli et al. 2015;Jones et al. 2017) impact they have on constitutive functions and their upscaling is currently an active field of research.
Assuming constitutive relations for relative permeabilities and capillary pressures are known, they can in principle be directly included in models with explicit fracture representations. However, differences in the composition of fracture and matrix may cause strong heterogeneities in the constitutive relations, manifested as capillary barriers and entry pressure effects, over the fracture-matrix interface. Depending on the specific form of relative permeabilities and capillary pressure functions, these heterogeneities are known to cause numerical challenges (Jaffré et al. 2011;Ahmed et al. 2017;Brenner et al. 2018).
The representation of multi-phase phenomena in continuum formulations is more challenging. As the transfer timescales are much slower than in single-phase flow, singlecontinuum representations are less likely to work because the local equilibrium assumption that underlies conventional relative permeabilities and capillary pressure relationships is not given. A numerical upscaling approach for parameters has, however, been outlined by Matthai and Nick (2009). The single-phase multi-continuum formulations have been expanded to incorporate multi-phase flow from early on (Kazemi et al. 1976;Gilman and Kazemi 1988) and continue to be developed (Lu et al. 2008;Tecklenburg et al. 2013;March et al. 2018). Comprehensive reviews are given by Lemonnier and Bourbiaux (2010a;2010b), Ramirez et al. (2009) andAl-Kobaisi et al. (2009). The fundamental challenge is that for a systematic development one would want to isolate different processes and develop transfer concepts individually while these processes are nonlinearly coupled.

Outlook: Multi-physics Couplings
This paper has reviewed the main modeling and discretization approaches for flow in porous media. However, what is often characteristic of fractured porous media is the strong mutual interaction between flow processes and the fractured structure itself. Flow processes can lead to deformation of the fractures due to mechanical, thermal and/or chemical processes, and altered fracture configurations again strongly impact flow. This is, for instance, the case considering membranes or thin layers in living tissues (Vernerey 2012). Another example is flow-induced seismicity, where flow reactivates and opens fractures that again severely affect flow patterns (Ucar et al. 2017).
The coupled thermal-hydraulic-mechanical-chemical (THMC) processes in the context of fractured porous media pose research challenges that have only partially been solved, despite strong research contributions over the last two decades (Rutqvist et al. 2002;Taron et al. 2009;Kolditz et al. 2016). For future research on flow in fractured media, inherent and remaining research challenges related to conceptual modeling and discretization approaches can be summarized in the following research questions: • How can disparate processes with different characteristic length and timescales be coupled in the mathematical models and discretization approaches presented herein? • How can macroscale representations of fractured porous media be amended to include not only changing fracture apertures but also propagating fractures? • In modeling of dissolution and precipitation processes, how should decisive phenomena such as clogging of fractures be modeled? • What are suitable upscaled models for interacting processes in fractured media?
• How can modeling be used for qualitative and quantitative insight when data on fracture configurations and information on constitutive relationships are poor?

Conclusion
This paper has reviewed mathematical and physical models and numerical approaches developed over the last decades to model flow in fractured porous media. While this is now a mature research field, providing modeling concepts and discretization approaches to handle a range of situations where networks of fractures dominate flow, substantial challenges remain in dealing with multi-phase flow and multi-physics couplings (e.g., multi-phase flow, mechanics, and chemistry).