1 Introduction

Naturally fractured reservoirs are often encountered in the exploitation of geothermal energy, hydrocarbon extraction, management of groundwater resources and carbon sequestration (Berkowitz 2002; Bonnet et al. 2001; Singhal and Gupta 2010; March et al. 2018). In all these applications, the transport of fluids is involved. Due to the geological complexity in naturally fractured reservoirs, numerical flow modelling methods are essential in fractured reservoir management. However, the highly heterogeneous nature of these reservoirs poses a significant challenge to numerical flow modelling. Direct representation of all fractures in a conventional reservoir simulator is simply too computationally expensive. In response, various fractured reservoir flow modelling methods have been proposed (Berre et al. 2018). These methods can be grouped into three categories: (1) continuum, (2) discrete fracture & matrix (DFM) and (3) hybrid methods. The main difference between the three groups of methods is how fractures are represented—implicitly or explicitly.

In continuum methods, fracture networks are represented implicitly as equivalent porous media. A widely used example of continuum methods is the dual-porosity (DP) method developed by Barenblatt et al. (1960). The DP method represents fractures as secondary porous media that interact with the primary porous media, the rock matrix. The two continua exchange fluids using transfer functions that capture the physics of cross-continuum flow (Lemonnier and Bourbiaux 2010a). Different transfer functions have been developed over time (Ramirez et al. 2009; Al-kobaisi et al. 2009). Here, we provide a brief overview of different transfer functions available and provide the equation for a transfer function that was used in our study. The first transfer function, introduced together with the DP method, was proposed by Warren and Root (1963). This transfer function related the transfer of a single-phase fluid between fractures and a matrix block to the difference in pressure between the two media. Kazemi et al. (1976) extended this transfer function to account for capillary forces, allowing the modelling of multiphase flow:

$$\begin{aligned} T_\alpha ^{\mathrm{DP}}=\sigma \frac{k_{\mathrm{m}} k_{{\mathrm{r}}\alpha }}{\mu _\alpha } \left( \psi _\alpha ^{{\rm f}}-\psi _\alpha ^{{\rm m}}\right) , \end{aligned}$$
(1)

where \(T_\alpha ^{\mathrm{DP}}\) is the fluid transfer from fracture to matrix, \(\psi _\alpha =p_\alpha +\rho _\alpha g d\) refers to flow potential, \(k_\mathrm{m}\) is an isotropic matrix permeability and \(k_{\mathrm{r}\alpha }\) is a relative permeability chosen based on upwind weighting. Finally, \(\sigma\) is a shape factor that depends on the geometry of the matrix blocks. While Eq. (1) uses flow potential as a driving force for fluid movement, the fact that flow potentials in the fracture and matrix are measured at the same point in space means that the gravity terms cancel each other out, thus negating the effects of gravity. Gilman and Kazemi (1983) then developed a transfer function that included matrix block saturation and height as inputs to explicitly account for gravity forces. While this formulation was able to include the effects of gravity, it assumed that gravity-driven flow was occurring on all faces of a matrix block, which is incorrect. To incorporate the directionality of gravity, Quandalle and Sabathier (1989) further developed a transfer function that explicitly accounted for fluid exchange mechanisms on each matrix block surface. Depending on the driving forces behind fracture–matrix fluid exchange, an appropriate transfer function has to be used. The transfer functions reviewed here all assume pseudo-steady state within matrix blocks—a state in which pressure and saturation are evenly distributed within the matrix block. Prior to reaching pseudo-steady-state flow, matrix flow initially goes through a transient state in which pressure and saturation distributions are higher near the matrix block boundaries and lower in the block centres. To account for transient matrix flow, a variation in the DP method known as the multiple interacting continua (MINC), introduced by Pruess (1985), should be used. The DP method also assumes that matrix blocks only exchange fluids with fractures, which is not always the case. When fracture apertures are small, capillary continuity may enable flow between matrix blocks. To account for inter-block flow, the dual-permeability (DK) method can be used (Lemonnier and Bourbiaux 2010b; Bourbiaux 2010).

A much simpler continuum approach is to lump the fractures with the matrix to form a new equivalent matrix, sometimes referred to as a pseudo-matrix (Rogers and Enachescu 2007). This single porosity (SP) approach enables conventional simulators to be used as is, but is only suitable for sparse and poorly connected fracture networks (Berre et al. 2018). The conversion of fracture networks into continua involves a process known as upscaling, which aims to derive a simplified representation of a problem that adequately reproduces the flow behaviour that would be observed in a fine scale representation of the same problem. A popular upscaling approach is to use steady-state flow fields to extract representative equivalent permeabilities. The steady-state flow fields can be numerically simulated or analytically solved for under suitable simplifications (Durlofsky 1991; Oda 1985; Sævik et al. 2013). The main benefit of continuum methods is that fracture network complexity is greatly simplified through the upscaling process, allowing for large-scale simulations to be performed. However, as we will show in this work, such simplifications are not always desirable because they oversimplify flow processes and lead to unreliable simulation outcomes.

In contrast to continuum methods, another family of techniques known as DFM methods explicitly represent the fracture networks in flow modelling. Many DFM methods exist, but they generally make the same simplifications to reduce computational costs (Flemisch et al. 2017). Firstly, the propensity for flow in fractures is captured using intrinsic fracture permeabilities. The cubic law relates the intrinsic fracture permeability \(k_\mathrm{f}\) and aperture a of a fracture through the relationship \(k_\mathrm{f}=a^2/12\) (Witherspoon et al. 1980). Secondly, fractures are modelled as lower-dimensional objects. Fractures are represented as lines in 2D domains and surfaces in 3D domains. These simplifications allow larger matrix grid cells to be used since grid cell sizes no longer have to be constrained by the space within the fractures. Larger grid cells have the advantage of reducing grid construction time and increasing the stability of explicit numerical schemes for solving fluid flow (Karimi-Fard et al. 2004). Closely related to DFM methods are discrete fracture network (DFN) methods. A DFN explicitly represents all fractures but assumes that matrix permeability is negligible (Long et al. 1982; Andersson and Dverstorp 1987; Hyman et al. 2015). Simulation of flow through the fracture network can be performed with DFN. Due to the impermeable matrix, isolated fractures in a DFN will not be able to contribute to fluid flow. In contrast, fractures in a DFM model may still communicate via the matrix even if they are not physically intersecting. We have generated a DFN for the purpose of this study, but did not use it for flow simulations. Instead, the DFN was used as an input to generate a DFM model.

One of the ways existing DFM methods differ from each other is grid conformity. Classical DFM methods require that fractures coincide with the interior boundaries of a matrix grid. Fractures are then represented as lower-dimensional cells that occupy the faces of matrix grid cells. As such, unstructured grids are required for the discretization of the matrix domain (Karimi-Fard et al. 2004; Paluszny et al. 2007). When using conforming grids, care must be taken to account for the change in capillary pressure and saturations across fracture–matrix interfaces. For DFM methods that employ finite element approaches, capillary pressures are defined for elements while pressure and saturations are calculated at cell nodes. As nodes are shared between fracture and matrix cells, it is necessary to increase the degrees of freedom at shared nodes to account for different saturations in the fractures and matrix (Reichenberger et al. 2006).

In the face of these issues, non-conforming DFM methods, which decouple the matrix grid construction from the fracture network, have been developed. The benefit of non-conforming grids is that they are much easier to construct. One such non-conforming DFM method is the eXtended finite element method (XFEM) which enriches the finite element solution space with basis functions that capture the effects of fractures. However, the XFEM approach assumes that fracture segments are through going within each finite element. As such, XFEM experiences stability issues at fracture tips (Schwenck et al. 2015). Another method which also uses the finite element method is the Lagrange multiplier method (LMM) introduced by Köppel et al. (2019). LMM does not enrich the finite element basis functions; instead, LMM represents fractures and the matrix in separate domains and allows them to interact via Lagrange multipliers that are to be solved for. Closely related to LMM is the embedded discrete fracture model (EDFM). Fluid exchange in EDFM is achieved through source/sink representations of fractures within the matrix grid. The advantage of EDFM over LMM is that EDFM is based on the finite volume method, making it compatible with commercial simulators such as Schlumberger’s ECLIPSE and Computer Modelling Group’s IMEX (Lee et al. 2001; Moinfar et al. 2012). While non-conforming DFM methods can be used with any matrix grid type, it draws its strength from its compatibility with structured matrix grids (Flemisch et al. 2017).

Another characteristic that differentiates DFM methods is the assumption of pressure continuity. For a conductive fracture which has an intrinsic permeability higher than the matrix permeability, the pressure field in the matrix is continuous across the fracture. However, for a sealed fracture that is less permeable than the matrix, the pressure field in the matrix is discontinuous across the fracture. Sealed fractures are not well modelled by all DFM methods. In particular, EDFM is not capable of capturing the discontinuity in pressure across low permeability fractures. This was shown by Flemisch et al. (2017) in their benchmarking of various DFM methods for a sealed fracture network. A recent improvement to EDFM was made by Tene et al. (2017) to account for sealing fractures; this variant of EDFM was named pEDFM. In our study, we only investigate fractures that are highly conductive in comparison with the matrix. As such, pressure continuity is not expected to impact our results.

The advantage of DFM methods is that they stay true to the fracture network geometry, allowing us to gain a better understanding of fluid flow dynamics in fractured reservoirs (Geiger and Matthäi 2014). However, the explicit representation of all fractures makes DFM simulations computationally intensive, especially for large-scale studies. Instead of a binary choice between implicit and explicit fracture representations, hybrid methods use a mixture of explicit and implicit fracture representations by simultaneously employing both continuum and DFM methods. The strength of hybrid models lies in their ability to capture long-range effects of large fractures, while simplifying the flow process in small fractures (Li and Lee 2008). The construction of hybrid methods involves two main decisions: partitioning and lumping. Partitioning is the process of splitting a fracture network into a set of small fractures that will be represented implicitly and another set of large fractures that will be explicitly represented. By using a suitable partitioning size, fractures smaller than the partitioning size will be modelled using continuum methods while larger fractures are modelled via DFM methods. Once partitioning is performed, implicitly represented fractures will have to be further subdivided such that each group of fractures can be incorporated into a continuum. We refer to this procedure as lumping. Implicitly represented fractures can be lumped into any number of continua. They can also be lumped together with the matrix to be upscaled into a new equivalent continuum, which is sometimes referred to as a pseudo-matrix. Depending on the partitioning and lumping choices, different hybrid methods can be used for flow simulation. The simplest hybrid methods are SP hybrid methods that combine the SP and DFM methods. The advantage of SP hybrid methods is that they can be applied using conventional reservoir simulators. However, they are only suitable for poorly connected fracture networks (Berre et al. 2018; Lee et al. 2001; Li and Lee 2008). DP hybrid methods which combine DP and DFM methods have been proposed by Jiang and Younis (2015), who also benchmarked the DP hybrid methods against fully resolved simulations. Their results showed that DP hybrid methods are suitable for implicitly representing well-connected fracture networks, while explicitly representing large-scale fractures. In this work, we name hybrid methods based on their component continuum and DFM methods (e.g. SP-DFM combines the SP and classical DFM methods while DP-EDFM combines the DP and EDFM methods). In practice, for full field simulations, both continuum and hybrid methods may be used since some or all of the fractures may be implicitly represented. For convenience, we refer to continuum and hybrid methods collectively as upscaled methods.

Since many upscaled methods are available for fractured reservoir simulations, it is often not obvious which method should be used for each application, and how the model should be constructed. Given a fractured reservoir, we have to first determine if any fractures have to be explicitly modelled. For fractures that will be implicitly represented, we also need to determine the number of continua needed, and which continuum each fracture should be lumped into. Depending on the partitioning and lumping choices made, a suitable upscaled flow model is then chosen and constructed. Partitioning and lumping choices are often made heuristically. For example, Lee et al. (2001) and Li and Lee (2008) proposed a SP-EDFM hybrid method that uses grid cell size to partition fractures. All fractures smaller than the grid cell size are lumped with the matrix to form a pseudo-matrix. On the other hand, Rogers and Enachescu (2007) studied flow in a chalk reservoir using the SP-DFM method by explicitly representing seismic scale fractures and lumping sub-seismic scale fractures with the matrix. However, there is no guarantee that these partitioning and lumping strategies will always work, especially if implicitly represented fractures are well connected (Berre et al. 2018). In terms of DP hybrid models, Jiang and Younis (2015) proposed and benchmarked the MINC-EDFM and DK-DFM methods. In their study, fractures were partitioned such that small-scale natural fractures are upscaled and large-scale natural or hydraulic fractures are explicitly represented. While benchmarking results have shown that this is an appropriate partitioning and lumping strategy, DP hybrid models may not always be necessary, especially if the small-scale fracture network is poorly connected. In such a case, SP hybrid models would be fit for purpose and simpler to implement. To further complicate matters, DP hybrid models require the selection of a transfer function that accounts for the correct matrix–fracture flow behaviour in NFRs; however, the flow behaviour cannot be known beforehand, making it difficult to make a selection (Berre et al. 2018; Lemonnier and Bourbiaux 2010a; Abushaikha and Gosselin 2013). In a more comprehensive approach, Bourbiaux (2010) proposed to choose between SP, DP, DK and hybrid methods using a heuristic evaluation scheme which compares fracture sizes to grid cell size, time scale of physical processes to simulation time steps and the connectivity of fractures. While useful, the evaluation scheme still relies on user defined parameters such as grid cell size and simulation time steps. Additionally, the time scale of physical processes may also be difficult to determine. These issues add an element of subjectivity in the application of the evaluation scheme.

Ideally, a rigorous evaluation would compare simulation results from different upscaled methods against a reference solution. Such a reference solution will have to be generated from a model that stays as close to reality as possible. In this regard, DFM methods are natural candidates as ’truth case’ generators. Unfortunately, DFM methods are unsuitable for full field simulations due to the prohibitive computational costs involved. Regardless, Geiger and Matthäi (2014) showed that DFM simulations performed on outcrop sourced fracture networks provide us with a deeper understanding of subsurface fluid flow. We thus hypothesize that our understanding of subsurface flow behaviour will benefit from simulation studies using small-scale DFM models representative of a naturally fractured reservoir; this information can then facilitate the selection of upscaled models (e.g. SP, DP, SP-EDFM or DP-EDFM) that are fit for large-scale fluid modelling. The flow characteristics observed can also guide the construction of the selected upscaled model.

It should also be noted that apart from using DFM simulations to gather insights regarding subsurface flow behaviour, there are also other analytical approaches that can be used. For example, dimensionless numbers such as the dual-porosity Damköhler number (DaDP) can be used. DaDP relates transfer and advective fluxes. The former describes fluid exchange rates between the matrix and fractures, while the latter describes fluid flow rates within the fracture network. DaDP is negative when advective fluxes are larger than transfer fluxes, zero when they are balanced and positive when advective fluxes are smaller than transfer fluxes. DaDP can be used to determine whether a DP model should be used. Calculation of the distribution of DaDP in a full field has been shown to take seconds (Spooner et al. 2019). One may also gain insights into the relative strengths of different transfer mechanisms between the matrix and fractures by analysing transfer rates. Lu et al. (2008) proposed to model matrix–fracture fluid transfers by superposing the effects of fluid expansion, viscous drive, diffusion, spontaneous imbibition and gravity drainage. The strength of each process is governed by transfer rates that depend on matrix block geometries and petrophysical properties. The relative magnitudes of these transfer rates allow us to identify the dominant drivers of fluid flow. Such analytical techniques are beyond the scope of this paper, but we advocate using such analytical techniques to corroborate the findings from DFM simulation studies.

In this paper, we propose a workflow that applies our hypothesis by using DFM simulations to qualitatively characterize flow behaviour in fractured reservoirs. By understanding the flow dynamics in a fractured reservoir, we can then select an upscaled flow modelling approach that is fit for purpose. Given a fractured reservoir, the proposed workflow begins with the creation of a DFM model that represents a small sector of the reservoir. Depending on the purpose of the final upscaled model, various simulation studies can be designed. For example, in carbon sequestration, March et al. (2018) highlighted the importance of designing CO\(_2\) injection rates to avoid early spill—a phenomenon where CO\(_2\) flows quickly through the fracture network with minimal displacement of brine in the matrix. Small-scale DFM models can then be used to investigate the sensitivity of CO\(_2\) sequestration to injection rate. In the case of hydrocarbon extraction from a fractured reservoir supported by an aquifer, it is important to prevent early water breakthrough due to high production rates. In this regard, small-scale DFM models can be used to study the flow behaviour of water at different production rates. By performing these simulations, the proposed workflow allows us to observe the effect of different fracture scales on fluid flow, how fractures interact and how matrix–fracture fluid exchange occur. These qualitative observations can then be used to select and construct upscaled flow models that are fit for purpose.

For our simulation studies, we have chosen to use EDFM, which is available in the MATLAB Reservoir Simulation Toolbox (MRST) (Lie et al. 2012; Shah et al. 2016). Flemisch et al. (2017) benchmarked EDFM against a full-dimensional fractured model for single-phase flow and found a good match between the two. However, to our knowledge, no such validation has been performed for more complex flow phenomena that involves capillarity or gravity. Nevertheless, EDFM has been used in many studies including that of Panfili et al. (2014) to simulate miscible sour gas injection into a fractured carbonate field, Siripatrachai et al. (2016) and Zhang et al. (2016) to understand the impact of capillary effects on hydrocarbon production in tight hydraulically fractured formations, Karvounis and Jenny (2016) for simulating flow in enhanced geothermal systems and Shakiba et al. (2018) to analyse production uncertainties in hydraulic fracture networks characterized by microseismic data. Since we are interested in a qualitative understanding of subsurface flow behaviour, and not a quantitative prediction of actual full reservoir behaviour, EDFM is adequate for our study. We chose EDFM for high-resolution simulations since it does not require grid conformity, thus simplifying the grid construction process. Moreover, as only conductive fractures are considered, we do not expect significant errors due to the assumption of pressure continuity in the EDFM method.

The rest of the paper is structured as follows: in Sect. 2, the proposed workflow is outlined in four steps, with further details provided for every step. The workflow is the main contribution of the paper. The rest of the paper then focuses on validating the workflow. Section 3 provides an overview of the validation procedure, which is to apply the workflow on a real gas reservoir and then perform error analysis on the chosen upscaled model. The first step in the application of the workflow is discussed in Sect. 4, where the features of the gas reservoir will be described. This section will show details of a sector scale EDFM model that captures the main features of the gas reservoir. The second step in the application of the workflow is discussed in Sect. 5, where a set of simulation studies is designed. Section 6 focuses on step 3 of the workflow. In this section, the simulation outputs will be analysed to improve our understanding of flow behaviour in the reservoir. Step 4 is discussed in Sect. 7, where the insights from the previous step are used to select an upscaled representation for the reservoir. Section 8 compares the selected upscaled model to reference solutions generated using EDFM. Additionally, other commonly used upscaled models are also studied. The comparison shows that EDFM studies on small-scale representative models allowed us to systematically select an upscaled model that adequately captured fluid flow behaviour and matched the reference EDFM model results. Finally, we summarize our findings in Sect. 9.

2 Flow Characterization Workflow

In this proposed workflow, we use DFM simulation studies to enhance our understanding of fluid flow behaviour in naturally fractured reservoirs. The obtained insights then allow us to decide how large-scale flow simulations can be performed. The workflow consists of four steps:

  1. 1.

    Construct a DFM model that represents a small sector of the reservoir being studied.

  2. 2.

    Determine a set of simulation studies to be performed on the representative model.

  3. 3.

    Analyse the simulation outputs to determine the physical processes influencing flow behaviour.

  4. 4.

    Choose an upscaled flow model that captures the highlighted physical processes.

In step 1, the DFM model created should capture the major features of a fracture network. While it is not possible to know the exact geometry of a subsurface fracture network, using information acquired from wellbores, seismic imaging and outcrops, we can stochastically generate discrete fracture networks (DFN) that honour the observed data. No fractures should be upscaled at this stage. However, we note that due to technological limitations, there will always be fracture scales that are un-observable during data acquisition. We assume that the conductivity of these fractures is similar to that of the surrounding rock. Thus, they can be treated as an integral part of the matrix. Apart from fracture network properties, the petrophysical properties used in the model should also be the same as those observed in the actual field. Since DFM methods are computationally intensive, the model will have to be restricted in size. Instead of a full field model, the constructed DFM model should be representative of a sector of the field. The DFM model can be on the scale of the drainage area of a water well, the distance between a geothermal well pair and the vertical distance between an oil production well and an aquifer. The reason for restricting the size of the DFM model is computational cost. If full field DFM simulations were feasible, then there would be no need for this proposed workflow. The computational cost of full field DFM simulations is the raison d’être of this proposed workflow.

Step 1 is often not straightforward due to the complexity of naturally fractured reservoirs. While it is not possible for us to cover all geological and physical possibilities, we will highlight some challenges that practitioners may encounter when applying the proposed workflow. For example, fracture network statistics such as density, average orientation, size ranges may be non-stationary, i.e. they vary spatially. Reasons for non-stationarity may be that different rock volumes have different mechanical properties and burial histories. In such cases, several DFM models should be considered in the workflow, each representing a zone (e.g. one rock volume) within which stationarity holds. The outcome of the workflow may then point to different upscaled models for each zone. In the case that a network’s fracture density is near the percolation threshold, stochastically generated DFNs may display a large variation in fracture connectivity. As such, it will be necessary to include multiple DFN realizations based on the same fracture network parameters in order to capture the potential variability of the fluid flow behaviour. Fracture networks may also display a high level of clustering. Under such circumstances, unstructured grids may be needed to sufficiently discretize tight spaces between fractures, while keeping grid cells reasonably coarse where fractures are not clustered. Finally, there may also be a lack of scale separation in fracture networks. In this case, it may not be possible to construct a sector model that adequately captures the behaviour of the field. For this situation, we suggest repeating the proposed workflow using increasingly larger models. The first model will be a DFM model. Subsequent models will be upscaled models. In each iteration, observed fluid flow behaviour can be used to determine how the next larger upscaled model should be constructed. The workflow can be repeated until the final upscaled model is on the scale of the full field. Testing all these ideas is beyond the scope of this paper. However, we hope that these suggestions will form the basis of further future research.

In step 2 of the proposed model, simulation studies are chosen and performed so that observations can be made about fluid flow behaviour. It is important to choose simulation studies that are relevant to the actual field application. For example, if we are investigating the displacement of gas by water from a bottom aquifer, the properties of the aquifer are important. As such, a DFM simulation study involving gas displacement by water injection will yield some insight into preferential fluid flow paths, but will not provide information relating the aquifer’s properties to fluid flow behaviour. The learnings obtained in a relevant simulation study will assist us in determining an upscaled flow model that is fit for large-scale flow simulations. An example of a simulation study relevant to aquifer driven gas production would be to investigate the amount of gas extracted given a range of production rates using a vertically oriented DFM model containing an aquifer model. Additionally, commonly encountered problems should also be considered. For example, given a sudden ambient temperature drop, demand for domestic heating increases. In this case, gas production may need to be increased to meet this rise in energy requirements. Such a scenario can also be included in step 2 to ensure that the chosen upscaled model accounts for the necessary physical processes if a similar study is required in the future for the full field.

After the simulations have been performed, the generated outputs are analysed in step 3. In particular, attention is paid to how different fracture scales influence flow patterns. Analyses should also be conducted to determine the driving forces behind fluid exchange between fractures and the matrix. Within the matrix, the distribution of saturation, pressure or solute concentration can also give an indication of the fluid flow regimes involved in fracture–matrix fluid transfer. These observations will enable us to make better informed partitioning and lumping choices in step 4. Additionally, if transfer functions are required for cross-continua flow, the analyses of fracture–matrix fluid exchange allow us to choose transfer functions that are fit for purpose. Once a partitioning and lumping strategy has been determined, a suitable upscaled flow modelling method is selected. For example, if all the fractures can be lumped with the matrix to form a pseudo-matrix, the SP method can be used. If the matrix and fractures cannot be lumped together, while the latter can be represented as a single continuum, then the DP method can be used. If some fractures have to be separately represented explicitly while the rest can be lumped together, then a DP hybrid method can be used.

We would like to further highlight that the proposed workflow is qualitative in nature. The DFM simulations are meant to improve our understanding of subsurface fluid flow behaviour in a reservoir of interest. This understanding can then be used to construct appropriate upscaled models.

3 Validation Methodology

To validate the proposed workflow, we applied it to a real fractured gas reservoir using the EDFM approach. Each step of the workflow was carried out as described in the previous section. To assess the performance of the full field upscaled model selected from the application of the proposed workflow, we would ideally compare it to a reference full field EDFM model. While running simulation for a full field upscaled model is possible, running a high-resolution EDFM simulation for a full field where all fractures are explicitly represented is prohibitively expensive. Since it is not possible to run a detailed full field EDFM simulation, we conduct two separate assessments on the small-scale EDFM model to study how well the selected upscaled model captures fluid flow behaviour in large-scale simulations. These studies are inspired by the work of Sablok and Aziz (2008), which considered upscaling errors as a combination of homogenization and discretization errors. Homogenization errors capture the impact of replacing highly heterogeneous permeability fields with less heterogeneous ones. On the other hand, discretization errors capture the effect of grid coarsening.

To calculate the homogenization error, we do not create upscaled models for comparison against the small-scale EDFM model. Instead, homogenized models which are as refined as the EDFM model are used. In our work, each homogenized model is created from the small-scale EDFM model in three steps. Firstly, the features to be homogenized are identified based on the lumping strategy used for hybrid or continuum modelling. Next, the identified features are upscaled to seek an equivalent permeability which captures the macroscopic flow behaviour of the features. The upscaling procedure is performed numerically using the method proposed by Durlofsky (1991). This method makes use of steady-state flow responses of incompressible fluid flow to calculate equivalent permeabilities from Darcy’s law. As such, we omit the effects of any aquifers, wells and gravity. We also temporarily substitute the original fluids in place with an idealized incompressible single-phase fluid. The simulation domain is the same as Fig. 1 but contains only features to be lumped.

In our model, the matrix is isotropic and the fractures are oriented either in the x- or y-directions. As such, we do not expect much crossflow to occur. Accordingly, the equivalent permeability is assumed to be a diagonal tensor. Moreover, since the fractures are equally long in both directions, the diagonal terms are also expected to be close to each other. To solve for the x-component of the equivalent permeability, a pressure differential is applied on opposing boundaries in the x-direction. For the remaining boundaries, a periodic boundary condition can be applied. Alternatively, since it is assumed that there are no off-diagonal terms in the equivalent permeability, a no-flow boundary condition can be used. In our work, we chose the latter approach since it is simpler to implement. The simulated flow field is then used to back-calculate an equivalent permeability using Darcy’s law. The y-component of the equivalent permeability can be found by repeating the process, but applying a pressure differential on opposing boundaries in the y-direction while flow is not permitted on the remaining boundaries.

The final step in the construction of a homogenized model is to replace the lumped features in the small-scale EDFM model with a homogeneous continuum. The permeability of the continuum is equal to the equivalent permeability previously calculated. The continuum will also be discretized at the same level of resolution as the small-scale EDFM model. Flow simulations can then be performed with the homogenized model (accounting for the actual fluids in place, gravity, aquifer and wells). For each production scenario, the results from flow simulations using homogenized models can be compared to their reference EDFM counterpart to calculate homogenization errors (Sect. 8.2).

The discretization error is calculated by comparing coarsened and homogenized models. The coarsened models are created from the homogenized models by re-discretizing the homogeneous continua using coarser grid cells. As was done to calculate homogenization errors, the coarsened models are subjected to flow simulations. The results of the flow simulations can then be compared to the results generated from homogenized models to compute errors resulting from grid coarsening. The calculation of discretization errors is discussed in Sect. 8.3.

Combining the findings from both studies allows us to assess the potential performance of an upscaled model for large-scale full field simulations. In addition to only comparing the upscaled model chosen using the proposed workflow in Sect. 2 against EDFM, we also compare two other popular upscaled models against EDFM to show that EDFM studies on fully resolved small-scale representative models enable us to choose an upscaled model with lower homogenization and discretization errors.

4 Workflow Step 1: Representative Model Construction

Our tests are run on a model representative of a small sector in a real naturally fractured gas reservoir with a bottom aquifer. The fracture network in the field is a result of geological folding, which gave rise to shear and tensile fractures that are, respectively, aligned and perpendicular to the direction of compression (Fig. 1). The tensile fractures are mainly concentrated at the crest of the fold. Borehole imaging results have shown that both fracture orientation sets show a degree of clustering which results in fracture corridors. The reservoir is dipped along the crest of the fold and has been producing under an aquifer drive. As the fractures are sub-vertical, the fracture network can be modelled in 2D. The 2D fracture network can then be oriented such that gravity points along one of its principal axes. This gas reservoir was chosen to test the proposed workflow due to various physical processes occurring during gas production. The presence of the aquifer results in two phase flow that involves multiple driving mechanisms (spontaneous imbibition, viscous and gravity drive). These processes occur on different time scales and will benefit from analysis using DFM methods like EDFM. Additionally, two fracture scales are present in the reservoir, which allows us to gain some insights into how the workflow performs given multiple length scales. The smaller fractures are highly dense and well connected, guaranteeing a separation of time scales in terms of fluid flow rates in the matrix and fracture network. The high fracture density also allows us to use only a single fracture network realization in our study as the variation in fracture connectivity is expected to be small.

Fig. 1
figure 1

EDFM model (right) is a representation of a small sector at the crest of a folded naturally fractured reservoir (left). The area of interest used to create the EDFM model is shown in red. The initial saturation field is shown in the EDFM model

Full size image

Based on proprietary field data, we used the method outlined in Priest (1993) to generate a simplified DFN (Fig. 1), which contains two orthogonal fracture sets. The fractures oriented along the x-axis represent the tensile fractures at the fold’s crest. The fractures oriented along the y-axis represent the shear fractures along the direction of compression. The fracture corridors are simplified and represented as large-scale high conductivity fractures. The remaining fractures, which we call diffuse fractures, are evenly dispersed in the domain. The parameters used for DFN generation are shown in Table 1. A minimum spacing was applied to the diffuse fractures to ensure that they are not artificially close to each other; this was done by imposing exclusion zones around fractures. New fractures are not allowed to intersect exclusion zones of existing fractures. The exclusion zone method is based on the observation of stress shadow zones around pre-existing fractures. The stress shadow zones prevent new fractures from forming near old fractures. The shapes of exclusion zones, for practical reasons, are chosen to be simple geometries that enclose fractures (Renshaw and Pollard 1994; Josnin et al. 2002). In our case, the exclusion zone used is rectangular in shape with dimensions 15.2 m \(\times\) 2.9 m. The fracture conductivities in Table 1 are related to the intrinsic fracture permeabilities, \(k_{\mathrm{f}}\), and apertures, a, by the relationship \(C_{\mathrm{f}}=k_{\mathrm{f}} a\) where the cubic law \(k_{\mathrm{f}}=a^2/12\) applies. As a result, the diffuse fractures have a uniform aperture of 61.9 \(\upmu\)m and the fracture corridors have a uniform equivalent aperture of 329 \(\upmu\)m. For the matrix, a porosity of 0.3 and permeability of 0.01 mD are used; these values are chosen to be on the same order of magnitude as field measurements of porosity and permeability. To account for the dip of the field, gravity was set to point in the x-direction. For simplicity, we do not adjust the strength of gravity based on the dip angle. In Sect. 6.1, we show that this simplification has minimal impact on the validity of our simulation studies because gravity does not contribute significantly to fracture–matrix fluid exchange. Next, an aquifer is specified at \(x=200\,\text {m}\) for which any fluid inflow is composed only of water. The aquifer pressure is fixed at \(p_{\mathrm{aq}}=100\) bar. For fluid production, a fixed rate horizontal well with 10 cm borehole radius is added to the EDFM model; the well is defined by the line connecting points (1.5 m, 40 m) and (1.5 m, 60 m). The well production rates used in this study are quoted in the percentage of gas in place (GIP) per year. GIP is calculated using the following equation:

$$\begin{aligned} \text {GIP}=\int \frac{S_{\mathrm{g}}}{B_{\mathrm{g}}}\phi \, \mathrm{d}V, \end{aligned}$$
(2)

where \(S_\mathrm{g}\), \(B_\mathrm{g}\) and \(\phi\) refer to gas saturation, gas formation volume factor and porosity. The integration is performed over both the matrix and fractures. Based on the current field production rate, the well is set to produce at 0.04 GIP/year.

Table 1 Parameters used to generate test model in Fig. 1
Full size table

The fluid phases flowing in this field are gas and water. The former is assumed to be an ideal gas while the latter is assumed to be incompressible. Their properties are shown in Table 2. Since there are two flowing media (fractures and matrix), two sets of relative permeability curves are generated for each medium. The fractures have straight line curves while the matrix relative permeabilities are modelled after Brooks and Corey (1964), which is as follows:

$$\begin{aligned} k_{\mathrm{rw}} = k_{\mathrm{rw}}^{\mathrm{max}} \left( \frac{S_{\mathrm{w}}-S_{\mathrm{wi}}}{1-S_{\mathrm{wi}}-S_{\mathrm{gi}}}\right) ^{n_\mathrm{w}}, \end{aligned}$$
(3)
$$\begin{aligned} k_{\mathrm{rg}} = k_{\mathrm{rg}}^{\mathrm{max}} \left( \frac{S_{\mathrm{g}}-S_{\mathrm{gi}}}{1-S_{\mathrm{wi}}-S_{\mathrm{gi}}}\right) ^{n_\mathrm{g}}, \end{aligned}$$
(4)

where \(k_{\mathrm{rw}}^{\mathrm{max}}\) and \(k_{\mathrm{rg}}^{\mathrm{max}}\) are the maximum relative permeabilities of water and gas. \(n_\mathrm{w}\) and \(n_\mathrm{g}\) are exponents that scale the relative permeability curves. Lastly, capillary forces are assumed to be negligible in the fractures since the apertures are relatively large (Firoozabadi and Hauge 1990; de la Porte et al. 2005). For the matrix, the relationship between capillary pressure and water saturation is again based on that of Brooks and Corey (1964), which is as follows

$$\begin{aligned} p_{\mathrm{c}}=p_{\mathrm{e}}\left( \frac{S_{\mathrm{w}}-S_{\mathrm{wi}}}{1-S_{\mathrm{wi}}-S_{\mathrm{gi}}}\right) ^{-n_\mathrm{p}}, \end{aligned}$$
(5)

where \(p_\mathrm{c}\), \(S_\mathrm{w}\), \(S_{\mathrm{wi}}\) and \(S_{\mathrm{gi}}\) are the capillary pressure, water saturation, irreducible water saturation and connate gas saturation. \(p_\mathrm{e}\) and \(n_\mathrm{p}\) are characteristics of the matrix; in this case, they have values of 1 bar and 1.5.

Table 2 Properties of different fluid phases
Full size table

For flow simulations, an EDFM test model is created using the inputs discussed above. A grid cell size of 1 m \(\times\) 1 m is applied to ensure that there is at least one full matrix grid cell between two diffuse fractures. Additionally, in Sect. 8.3, grid convergence is tested to show that the grid cell size chosen produces accurate results. Prior to conducting any simulation studies, a capillary transition zone is set up within the test model by imposing hydrostatic equilibrium (Fig. 1). Given the aquifer pressure and fluid densities, the hydrostatic pressure profile for gas and water above the aquifer is first calculated. Capillary pressure in the matrix is then calculated using the gas and water phase pressure differences. Finally, water saturation is calculated using Eq. (5). Since capillary pressure is negligible in the fractures, they are initially fully saturated with gas.

5 Workflow Step 2: Determine DFM Studies

Using the test model created in the previous section, we perform a series of studies that mirror some common investigations carried out for producing gas fields. The results from these simulation studies will allow us to better characterize fluid flow behaviour in the fractured gas reservoir (Fig. 1). For full field simulations, we can then fine tune our upscaling approach based on the observed flow characteristics. As a result, we can then arrive at an upscaled model that captures the small-scale physics.

Here, we outline the EDFM simulation studies that we performed (Table 3). In line with Geiger and Matthäi (2014), the purpose of these high-resolution EDFM studies is to better understand the flow behaviour that will be observed in a field. As such, it should be noted the studies in Table 3 are not the only studies that can be used when using high-resolution DFM simulations to characterize fluid flow behaviour. In particular, the types of studies that will be performed on the full field should dictate the types of high-resolution DFM studies to be conducted. While we are applying our methodology to a fractured gas reservoir, the simulation studies we have chosen to perform in this work are applicable to many types of fractured reservoir exploitations.

The goal of study 1 is to understand the contribution of viscous, gravity and capillary forces to hydrocarbon recovery. The results of the analysis will be used to determine if fractures can be lumped with the matrix when considering a hybrid or continuum modelling approach. If not, the observations will help to determine the appropriate transfer functions that incorporate the observed physical processes behind flow.

Table 3 EDFM simulation studies performed to characterize flow in a naturally fractured gas reservoir
Full size table

In study 2, we investigate the effect of stopping production using the fully resolved EDFM test model. The observed fluid flow dynamics are then used to guide the selection of an upscaled model. When extracting gas from a reservoir, produced formation water will have to be treated before disposal. Moreover, water may accumulate in the well and inhibit further gas production. A common remedial action taken when water breaks through at a well is to stop production. Allowing the reservoir to equilibrate will prevent water production when the gas production is eventually restarted. While it is unknown if early water breakthrough will happen in the naturally fractured gas field, using study 2, we take a pre-emptive approach by ensuring that the full field upscaled model is suitable for simulating well shut-ins if the need ever arises.

In study 3, the effect of production rates on flow behaviour is investigated. This study is often conducted for full field models to understand the effects of higher well rates. While high well rates are generally appealing, they can also lead to some undesirable effects. In our particular case, for the fractured gas reservoir, it is enticing to produce gas at a higher rate to increase profitability. Production supported by a bottom aquifer usually faces the risk of uneven water progression through the reservoir. Water progresses faster near the production well and slower further away. As a result, the gas–water interface resembles a cone, with the crest of the cone below production well. In the hydrocarbon industry, this phenomenon is referred to as water coning. Water coning is undesirable in the extraction of hydrocarbons as it leads to early water breakthrough. As such, in full field simulations, it is expected that various well rates will be tested to understand their impact on gas recovery at water breakthrough. In preparation for such investigations, study 3 focuses on simulating water coning phenomena using the small-scale EDFM test model. The simulation results will help to determine the type of upscaled model that should be used for full field studies.

6 Workflow Step 3: Analysis of Results from Simulation Studies

6.1 Study 1: Driving Forces

In this study, the driving forces behind hydrocarbon recovery from the naturally fractured gas field are investigated. Since the EDFM test model created in Sect. 4 readily includes viscosity, gravity and capillary forces, we perform a drawdown simulation using a rate of 0.04 GIP/year (Case 1a). The production rate is chosen such that the gas might be fully extracted from the model in 25 years. The simulator is allowed to march in time until water breaks through at the production well. The breakthrough time for water was found to be 10.75 years, with a recovery factor at breakthrough of 43%. The recovery factor is the volume of gas that has been produced, normalized by the GIP in the model. The evolution of the water saturation field is shown in Figs. 2 and 3, where a small amount of coning can be observed. Water from the aquifer preferentially flows through the fracture corridors and spreads into the diffuse fractures. Water is then drawn into the matrix through capillary and gravity forces; this slows down the progression of water in the fractures and results in a stable coned saturation front. It is also observed that as water imbibes the matrix from the fractures, water saturation in the matrix is high near the fractures and low in the centre of the matrix blocks. The difference in matrix water saturation between matrix block edges and centres suggests that matrix flow remains in transient state the entire time. The long transient matrix flow regime has implications on the choice of transfer functions that can be used in the DP method. Conventional DP transfer functions assume pseudo-steady-state flow. For short transient matrix flow regimes, this assumption is acceptable as matrix flow will be in a pseudo-steady-state regime for most of the time. With long transient matrix flow regimes, this assumption has been shown to result in an underestimation of matrix–fracture transfer rates. To account for transient matrix flow regimes, a MINC approach can be used. Alternatively, transfer functions that account for early time flow behaviour can be used (March et al. 2016).

Fig. 2
figure 2

Fracture water saturation maps at different stages of production from the test model at a base rate of 0.04 GIP/year (Case 1a). The recovery factors represent the fraction of GIP in the test model that has been produced

Full size image
Fig. 3
figure 3

Matrix water saturation maps at different stages of production from the test model at a base rate of 0.04 GIP/year (Case 1a). The recovery factors represent the fraction of GIP in the test model that has been produced

Full size image

To further understand the impact of capillary and gravity forces, the test model was modified to produce two other models: one with neither gravity nor capillary forces, and one with gravity included but without capillary forces. These models correspond to cases 1b and 1c. In case 1b, the model is initialized with constant saturation (\(S_{\mathrm{w}}=1-S_{\mathrm{gi}}\)) and pressure (\(p=p_{\mathrm{aq}}\)) fields. In case 1c, the model is initialized with a constant saturation (\(S_{\mathrm{w}}=1-S_{\mathrm{gi}}\)) field at hydrostatic equilibrium. As was done in case 1a, both models in cases 1b and 1c are produced at the same rate of 0.04 GIP/year. However, note that these models have different GIPs compared to the base case since there is no capillary transition zone.

The water saturation maps at breakthrough in cases 1b and 1c are shown in Fig. 4. With neither gravity nor capillary forces, water breakthrough in case 1b is immediate and no gas is produced; water progresses uninhibited to the production well through the fracture corridors. In case 1c, gravity facilitates the entry of water into the diffuse fractures and slows down the progression of the water saturation front. Gas is allowed the opportunity to decompress and be produced. However, despite the full strength of gravity being applied, very little gas in the matrix is displaced by water. As a result, the recovery at breakthrough is increased by only 12.5%. Recalling the results in Figs. 2 and 3, including capillary forces enables spontaneous imbibition of water into the matrix, increasing the recovery at breakthrough even more.

Fig. 4
figure 4

Water saturation maps (fractures in left column, matrix in right column) at breakthrough with a well production rate of 0.04 GIP/year. a Case 1b with neither capillary nor gravitational forces: breakthrough recovery at 0%, b Case 1c with capillary forces included: breakthrough recovery at 12.5%

Full size image

The findings in this study are summarized as follows:

  1. 1.

    Water preferentially flows through the fracture network, but its progression is delayed as water in the fractures imbibes the matrix blocks.

  2. 2.

    Water saturation in the matrix blocks is high near the fractures and low in block centres, indicating that the transient flow regime in the matrix blocks is long.

  3. 3.

    Gravity does not contribute significantly to matrix imbibition.

  4. 4.

    Capillary forces are the main driver for gas recovery from the matrix through spontaneous imbibition.

These findings have implications in terms of choosing an appropriate upscaled model for the fractured gas field. Point 1 shows that fracture–matrix fluid exchange needs to be accurately modelled to capture the dynamics of gas and water production in this field. Point 2 shows that the transient process of matrix–fracture fluid exchange may have to be modelled using the MINC method or a transfer function that accounts for transient flow. Points 3 and 4 show that if a DP model is chosen to simulate flow in this field, a transfer function that accounts for capillary forces is necessary while gravity may be neglected. A simple example of such a transfer function is that by Kazemi et al. (1976). However, this transfer function only accounts for pseudo-steady-state flow in the matrix. Alternatively, the transfer function by March et al. (2016) can be used instead.

6.2 Study 2: Well Shut-In

In this study, EDFM is used to gain an understanding of the field’s response to shut-in time after initial water breakthrough. During the shut-in duration, the fractured reservoir is expected to equilibrate. In the case of the gas field, water in the fracture network is expected to continue imbibing the matrix blocks until the fracture network is free of water. In this study, the EDFM test model is produced at a well rate of 0.04 GIP/year until water breaks through at 10.75 years. The well is then shut in for 5 years before being restarted again at 0.04 GIP/year. The simulator is then allowed to continue marching in time.

The water saturation maps at 2 and 5 years for the fractures and matrix after shutting in the well are shown in Fig. 5. It is observed that once the well is shut-in, the water saturation front in the fracture network settles and becomes more dispersed. The front recedes as water imbibes previously unexposed matrix blocks upon contact. However, despite several years of shut-in time, water saturation in the tight matrix remains high near fractures and low in the centre of the matrix blocks. Once matrix block edges are saturated with water, water imbibition from the fractures to the matrix is impeded, causing the fracture water saturation front to recede very slowly. Since there is still a significant amount of water in the fracture network after the end of the shut-in period, water breakthrough in the second phase of production occurs quickly after only 1 year since restarting the well.

Fig. 5
figure 5

Water saturation maps (fractures in left column, matrix in right column) at different times during the shut-in period. a 2 years, b 5 years

Full size image

The results obtained in this study have direct implications on the selection and construction of an upscaled model for the fractured gas field. In particular, the transient matrix flow behaviour needs to be accounted for in the upscaled model. March et al. (2016) showed that using transfer functions that do not account for transient flow in the matrix results in initially lower water imbibition rates. However, the matrix blocks reach full water saturation faster. In fractured systems with water wet or highly permeable matrices, these transfer functions are often applicable since the transient state is usually short and negligible. For the considered reservoir, however, shut-in periods will have to be simulated using upscaled models that account for the transient flow regime because the transient flow regime is long.

6.3 Study 3: Different Production Rates

It is expected that the upscaled model will be used to study the effect of production rate on water breakthrough time. In this study, the EDFM test model is used to simulate gas production under different well rates. The EDFM test model is produced using four different well rates: 0.04 GIP/year, 0.1 GIP/year, 0.3 GIP/year and 1 GIP/year. The lowest production rate was chosen such that the gas in place might be produced in 25 years. The other production rates are chosen through logarithmic step increases. The respective times and recoveries at breakthrough are shown in Table 4. It can be seen that recovery at breakthrough reduces as gas production rate increases (Fig. 7), implying that more gas has been bypassed due to coning effects.

Table 4 Breakthrough times and recoveries under different production rates
Full size table

At high production rates of 0.1 GIP/year and 0.3 GIP/year (Fig. 6a, b), coning effects become more severe than producing at 0.04 GIP/year (Figs. 2, 3), but water still manages to displace gas from a significant number of diffuse fractures. Water in the diffuse fractures then further imbibes the matrix. At the production rate of 1 GIP/year (Fig. 6c), coning becomes severe. Water bypasses most of the diffuse fractures and progresses directly to the well through the fracture corridors. As a result, there is barely any water imbibition into the matrix, resulting in lower recoveries.

Fig. 6
figure 6

Water saturation maps (fractures in left column, matrix in right column) at different well production rates. a 0.1 GIP/year, b 0.3 GIP/year and c 1 GIP/year

Full size image

Through this study, we showed that water from the aquifer preferentially flows through the fracture corridors. Water then enters the diffuse fracture network through the fracture corridors. Upon contact with the matrix, water imbibes the matrix from the diffuse fractures. At low production rates, there is enough time for water to imbibe the matrix through the diffuse fractures. On the other hand, at high production rates, the water saturation front reaches the production well through the fracture corridors before water is able to displace gas from the diffuse fractures. This observation suggests that, for the fractured gas field, the selected upscaled model should separately represent fracture corridors and diffuse fractures. To accurately study the dependency of gas recovery on production rates, a DP-EDFM model in which the fracture corridors are explicitly represented is likely more suitable than a DP model which lumps fracture corridors and diffuse fractures together.

7 Workflow Step 4: Selection of Upscaled Representation

The observations made in the three simulation studies allow us to make an informed decision as to how the fractured gas reservoir in question should be represented in an upscaled model.

Firstly, study 3 highlighted that fracture corridors and diffuse fractures have to be separately represented in the upscaled model. The separate representation is important if full field sensitivity studies are expected to be performed to explore the impact of different production rates. As such, the fracture network should be partitioned such that the fracture corridors are explicitly represented while the diffuse fractures are implicitly represented.

Studies 1 and 3 also highlighted that water imbibition into the matrix occurs through the interface between diffuse fractures and the matrix. Hence, the diffuse fractures and matrix cannot be lumped together into a pseudo-matrix since such a representation would not reproduce the observed flow path. Explicitly representing all the diffuse fractures in a full field model would be unrealistic. Instead, they should be lumped and represented implicitly in a continuum.

Study 2 showed that when the production well is shut in, water redistributes itself in the fracture network to seek hydrostatic equilibrium. Consequently, water comes into contact with more matrix blocks (through the diffuse fractures) and imbibes them. More importantly, it was also observed that fluid flow in the matrix remained in transient state for the entire shut-in duration. In order for the upscaled full field model to reproduce this phenomenon, simulation approaches that account for transient matrix flow need to be used.

Finally, study 1 showed that the exchange of fluid between diffuse fractures and the matrix is mainly driven by capillary pressure and to a small extent by gravity. As such, the transfer function used to model flow between diffuse fractures and the matrix need only account for capillary forces in order to be fit for purpose. The simplest transfer function that accounts for capillary forces is that by Kazemi et al. (1976).

Given the preceding analysis, we concluded that the DP-EDFM method should be used for full field flow simulations in the fractured gas reservoir. In this model, the fracture corridors should be explicitly represented, while the diffuse fractures should be upscaled into a continuum. Flow between the diffuse fractures and matrix should be modelled using transfer functions that account for capillary pressure and transient matrix flow. An example of such a transfer function is that by March et al. (2016). Alternatively, the MINC-EDFM method can also be used to account for the transient matrix flow behaviour.

8 Error Quantification: Comparison of Upscaled Models to EDFM

In this section, we seek to study the potential impacts of model selection on simulation outcomes by comparing EDFM results to those generated using DP-EDFM, DP and SP-EDFM models. As discussed in Sect. 7, DP-EDFM was identified as an appropriate approach to simulate flow in the fractured gas field. The DP method was included in the study to illustrate the effect of lumping all fractures together, regardless of size, into one continuum. The SP-EDFM method was included to show the effect of lumping diffuse fractures together with the matrix. The comparison between EDFM and the upscaled models is performed through two studies as outlined in Sect. 3. The first study investigates errors arising from homogenization by comparing EDFM and the homogenized models. The second study investigates discretization errors arising from grid coarsening by comparing coarsened models against homogenized models. The results of both studies, in conjunction, allow us to assess whether or not an upscaled model would be suitable for large-scale simulations. In both studies, models constructed are subject to production under the following rates: 0.04 GIP/year, 0.1 GIP/year, 0.3 GIP/year and 1 GIP/year. We then use the recovery and saturation fields at breakthrough predicted by the models for comparison.

8.1 Construction of Homogenized and Coarsened Models

Here, we describe in detail the procedures used to construct the homogenized and coarsened DP, SP-EDFM and DP-EDFM models. For each model, the properties of every continuum are calculated. Gravity is set to point in the x-direction—similar to the reference EDFM model. A production well is specified at the same location as shown in Fig. 1. An aquifer is positioned at \(x=200\) m with a constant pressure of \(p_{\mathrm{aq}}=100\) bar. These models will then be used to generate flow responses that can be compared to the fully resolved EDFM solutions previously discussed in Sect. 5.

8.1.1 DP Models

In the homogenized DP model, the fracture corridors and diffuse fractures are lumped together and homogenized using the procedure outlined in Sect. 3. Since the fracture network is isotropic, the equivalent permeabilities calculated for both the x- and y-directions are similar. The resulting averaged equivalent permeability of the fracture continuum is 32.27 mD.

To complete the homogenized DP model, petrophysical properties and a transfer function have to be specified. The fracture continuum relative permeability curves used are the same as shown in Table 2, and capillary pressure remains negligible. The matrix continuum retains the same properties as the matrix in our EDFM test model. As discussed in Sect. 7, a transfer function which accounts for transient matrix flow should be used. However, such a transfer function was not available in MRST at the time of this work. As a compromise, we chose to use the transfer function by Kazemi et al. (1976), which accounts for capillary pressure but assumes pseudo-steady-state matrix flow. To calculate the shape factor needed to parametrize the transfer function, the size of the matrix block is required. The fracture spacing is often used as the matrix block size. However, since both fracture corridors and diffuse fractures (which have different fracture spacings) were lumped into the fracture continuum, an equivalent block size should be used instead. We found that a block size of 3.3 m allows the homogenized DP model to match the EDFM reference model reasonably in terms of recovery at breakthrough (Fig. 7). The block size of 3.3 m is exactly the average fracture spacing of the diffuse fractures (\(1/P_{21}\)). Due to the large number of diffuse fractures in comparison to the fracture corridors, it is not surprising that using the diffuse fracture spacing as the equivalent matrix block size allows the homogenized DP model to match reference EDFM solutions.

The homogenized DP model can be compared to the fully resolved EDFM model to quantify the homogenization error arising from lumping all fractures into a continuum. To understand the effect of grid coarsening, coarsened DP models have to be compared to the homogenized DP model. As such, five coarsened DP models were created using different grid cell sizes: 2 m, 4 m, 8 m, 16 m and 32 m.

8.1.2 SP-EDFM Models

In the homogenized SP-EDFM model, fracture corridors are explicitly modelled, while the diffuse fractures are lumped together with the matrix to form a pseudo-matrix. The homogenized pseudo-matrix is constructed using the procedure outlined in Sect. 3. As the diffuse fracture network and the matrix permeability are isotropic, the components of the equivalent permeability in both directions are similar. The final equivalent permeability of the pseudo-matrix is 1.869 mD (two orders of magnitude higher than the matrix permeability). In the SP-EDFM approach, the explicitly modelled fracture corridors communicate with the pseudo-matrix using the EDFM approach of establishing matrix–fracture connections based on the intersections between matrix and fracture cells.

Since the matrix and diffuse fractures have different relative permeabilities and capillary pressures, equivalent versions of these properties will need to be calculated for the pseudo-matrix. This procedure is known as pseudoization, and the results are known as pseudorelative permeabilities and pseudocapillary pressures. These pseudoized properties represent the bulk behaviour of the pseudo-matrix under multiphase flow. While pseudoization is a commonly used approach, this technique has been shown to be an unreliable tool for upscaling high-resolution models (Barker and Thibeau 1997). As such, in our case, we argue that since the volume ratio of matrix to diffuse fractures is large, we simply let the pseudo-matrix inherit the petrophysical properties (relative permeability and capillary pressure) of the matrix. This simplistic treatment will later be shown to result in a reasonably accurate homogenized SP-EDFM model when compared to the fully resolved EDFM model. Similar to the coarsened DP models, five coarsened SP-EDFM models were created using different grid sizes (2 m, 4 m, 8 m, 16 m and 32 m) for studying discretization errors.

8.1.3 DP-EDFM Models

The homogenized DP-EDFM model represents the fracture corridors explicitly while the diffuse fractures are independently upscaled into a fracture continuum (no lumping with the matrix). The homogenized fracture continuum is constructed using the procedure outlined in Sect. 3. As the diffuse fracture network is isotropic, both the x- and y-directional components of the equivalent permeability are similar. In the homogenized DP-EDFM model, the average of the two values is used for the isotropic equivalent permeability. The final calculated fracture continuum permeability is 1.801 mD. The calculated fracture continuum permeability is very close to but lower than the pseudo-matrix permeability used in the homogenized SP-EDFM model. This is not surprising since the matrix permeability is very tight at 0.01 mD. Including the matrix in the upscaling process would thus only increase the equivalent permeability marginally.

In terms of petrophysical properties, since all diffuse fractures have the same capillary pressure and relative permeability curves, these curves are also applied to the diffuse fracture continuum (zero capillary pressure and straight line relative permeabilities with zero residual saturations). Fluid flow between the explicitly modelled fracture corridors and the fracture continuum (which accounts for the small fractures) is accounted for using the EDFM approach by Moinfar et al. (2013). In the original EDFM approach, connections are established between the matrix and explicit fractures. Instead, for DP-EDFM, we replace all matrix properties in the original EDFM formulation with properties of the fracture continuum. This small modification allows us to establish fracture corridor-to-fracture continuum connections based on the intersections between fracture corridor and fracture continuum grid cells. In our case, direct interaction between the explicitly modelled fracture corridors and the matrix is neglected since the previous case studies have shown that transfer from the fracture network to the matrix mainly occurs through the diffuse fractures. The interaction between the fracture continuum and matrix is based on the DP approach (Warren and Root 1963). As was done for the homogenized DP model, fluid exchange between the fracture continuum and matrix is modelled using the transfer function by Kazemi et al. (1976) with a block size of 3.3 m (calculated from the inverse of diffuse fracture density). No calibration of the block size is performed since fracture corridors were not lumped into the fracture continuum. Additionally, five coarsened DP-EDFM models were created using different grid sizes (2 m, 4 m, 8 m, 16 m and 32 m) for studying errors arising from grid coarsening.

8.1.4 Coarsened EDFM Models

In addition to creating coarsened DP, SP-EDFM and DP-EDFM models, the same is done for the EDFM model in order to investigate the performance of EDFM as grid cell sizes increase. Five coarsened EDFM models were created with grid cell sizes of 2 m, 4 m, 8 m, 16 m and 32 m. Production well and aquifer locations are similar to the reference EDFM model in Sect. 4.

8.2 Effect of Homogenization

To investigate the effect of homogenization on model performance, simulation results generated by homogenized models are compared to the reference EDFM solutions which are shown in Sect. 6.3. The gas recoveries at breakthrough predicted by the homogenized DP, SP-EDFM and DP-EDFM models are compared against EDFM results (Fig. 7). The well bottom hole pressures are also compared (Fig. 8). Overall, the homogenized DP model performs well since the equivalent block size was calibrated to match the EDFM results. The homogenized SP-EDFM and DP-EDFM models, respectively, over- and under-predict recovery at breakthrough compared to the EDFM results. In most cases, DP-EDFM models also predict larger well pressure drops compared to the reference EDFM results. On the other hand, SP-EDFM models under-predict the pressure drop in the well.

Fig. 7
figure 7

Recoveries at breakthrough predicted by homogenized DP, SP-EDFM and DP-EDFM models at the same level of grid refinement. The reference solution is produced using EDFM

Full size image
Fig. 8
figure 8

Well bottom hole pressure evolution at different well rates predicted by homogenized DP, SP-EDFM and DP-EDFM models at the same level of grid refinement. The reference solution is produced using EDFM. a 0.04 GIP/year, b 0.1 GIP/year, c 0.3 GIP/year, and d 1 GIP/year

Full size image

The comparisons of recovery at breakthrough and well bottom hole pressure may give the impression that DP models are adequate for flow modelling in this field. However, the homogenized DP model does not produce saturation fields that match the reference EDFM results at all since it lumps all fractures together in a singular representation. Figure 9 shows the saturation field at breakthrough for a production rate of 0.3 GIP/year, which do not match the saturation profile shown in Fig. 6b. On the contrary, by explicitly modelling the fracture corridors, the homogenized SP-EDFM and DP-EDFM models produce saturation fields that better resemble the reference EDFM results. Figures 10 and 11 show the saturation profiles at breakthrough for the homogenized SP-EDFM and DP-EDFM models with a production rate of 0.3 GIP/year. The saturation profiles are able to capture the saturation profile simulated by the EDFM model (Fig. 6b). Between the two hybrid approaches, DP-EDFM best captures the overall saturation distribution in the fracture corridors, diffuse fractures and matrix.

Fig. 9
figure 9

Water saturation maps (left to right: fractures, matrix) at breakthrough for the homogenized DP model with an equivalent matrix block size of 3.3 m. The well production rate is 0.3 GIP/year

Full size image
Fig. 10
figure 10

Water saturation maps (left to right: fracture corridors, pseudo-matrix) at breakthrough for the homogenized SP-EDFM. The well production rate is 0.3 GIP/year

Full size image
Fig. 11
figure 11

Water saturation maps (left to right: fracture corridors, diffuse fracture continuum, matrix) at breakthrough for the homogenized DP-EDFM model. The well production rate is 0.3 GIP/year

Full size image

Despite producing saturation profiles that better match EDFM results, SP-EDFM and DP-EDFM were not able to match the recovery at breakthrough and well bottom hole pressure mainly due to the way fluid transfer between diffuse fractures and the matrix is modelled. In the homogenized SP-EDFM model, the matrix and diffuse fractures have been lumped together. Without separately representing the diffuse fractures and the matrix, water that enters the diffuse fractures immediately imbibes the matrix. As a result, the matrix is fully saturated with water behind the water saturation front. The fast imbibition of water into the matrix then impedes the movement of the water saturation front in the fracture corridors. Thus, the water saturation front takes a longer time to reach the production well and increases recovery at breakthrough. Accordingly, well bottom hole pressure drops slowly since the fracture corridors contain more gas than predicted by the EDFM model at any given point in time. On the other hand, in the homogenized DP-EDFM model, fluid exchange between diffuse fractures and the matrix is modelled using a transfer function that assumes pseudo-steady state in the matrix. This transfer function allows the homogenized DP-EDFM model to better capture the difference between advective and transfer fluxes. As a result, the matrix is not fully saturated behind the water saturation front. However, as shown in the studies in Sect. 5, matrix flow remains in a transient state for the entirety of the simulation period. Thus, the pseudo-steady-state assumption is not applicable. Without accounting for the transient matrix fluid flow, less water imbibes the matrix in the early time regime (March et al. 2016). Water then moves through the fracture network less inhibited in the homogenized DP-EDFM model compared to the EDFM model. Consequently, the water saturation front reaches the production well earlier. Hence, the well pressure reduces faster as the fracture corridors contain less gas than predicted by the EDFM model at any given point in time.

To study the errors arising from homogenization, one first needs to define an error measure. Several possibilities include comparing the recoveries at breakthrough, bottom hole pressures, pressure distributions or saturation distributions. The error measure chosen should correspond to the reason for running such flow simulations. In this test case, the main concern with gas production from an aquifer driven reservoir is to prevent the production of aquifer water. In the event that water does break through at the production well, the next course of action is often to shut in the production well and drill an infill well in a location where water production is unlikely in the short term. Since the time it takes for water to break through in the new well is directly linked to the saturation distribution at water breakthrough in the first well, it is our opinion that the comparison of saturation distributions should take precedence over other forms of error measures. To determine errors based on saturation distributions, ideally one would calculate the difference between saturation fields generated from EDFM and upscaled models. However, such a comparison would not be straightforward because EDFM, DP, SP-EDFM and DP-EDFM grids may have the same average grid cell size, but will have different number of grid cells. Take for example, the number of grid cells in an EDFM model will be the number of matrix grid cells plus the number of discrete fracture grid cells. The number of grid cells in a DP model will be the number of matrix grid cells plus the number of fracture continuum grid cells (which is the same as the number of matrix grid cells). Comparing EDFM versus DP, EDFM versus SP-EDFM, EDFM versus DP-EDFM will all require different types of saturation field conversions. The way the conversions are performed is subjective. Homogenization errors measured this way are thus difficult to compare. Alternatively, it is likely that the continuing the simulation beyond water breakthrough might further accentuate the inadequacy of the DP model due to the poor saturation profiles predicted. Such a study will help quantify the homogenization errors introduced through DP modelling. We have not performed such simulations since, in reality, gas wells are often shut in upon water intake. Instead, we propose a homogenization error measure inspired by the performance of a new production well. The error measure involves picking a location for a new well and then measuring the shortest distance from the well to the water saturation front at breakthrough, \(d_{\mathrm{ws}}\). More precisely, \(d_{\mathrm{ws}}\) is the solution of the following problem.

$$\begin{aligned} \begin{array}{lll} \text {maximize}&{} \quad d \\ \text {subject to}&{} \quad S_{{\rm w}}(x)=S_{{\rm wi}}(x), \quad \forall x \in \{x\in \Omega \,|\, D(x,x_{{\rm well}})\le d\}, \\ &{}\quad d\ge 0, \end{array} \end{aligned}$$

where \(\Omega\) refers to the simulation domain, \(x_{\mathrm{well}}\in \Omega\) refers to the location of a new well and D(x, y) refers to the Euclidean distance between points x and y. Note that \(d_{\mathrm{ws}}\) is a function of \(x_{\mathrm{well}}\) as well as the model used. With the same \(x_{\mathrm{well}}\), different models will yield different \(d_{\mathrm{ws}}\) values. In this study, \(x_{\mathrm{well}}\) is chosen as the solution of the following problem, applied to the EDFM model.

$$\begin{aligned}&\underset{x_{\mathrm{well}}}{\text {maximize}} \quad A(\{x\in \Omega \,|\, D(x,x_{\mathrm{well}})\le d_{\mathrm{ws}}(x_{\mathrm{well}})\}) \\&\text {subject to}\quad x_{\mathrm{well}}\in \Omega , \end{aligned}$$

where \(A(\cdot )\) is a function that takes in a subdomain of \(\Omega\) and outputs the area of the subdomain. This approach to selecting \(x_{\mathrm{well}}\) is analogous to picking a new well position to maximize gas recovery from the well before water breaks through. As an example, in the case of a production rate of 0.04 GIP/year, the new well location and distances to the saturation front for each upscaled model are illustrated in Fig. 12. The new well location is indicated by a red dot, while the red circles have radii equal to \(d_{\mathrm{ws}}\). Note that the well location was chosen by maximizing the area within the red circle in Fig. 12a. \(d_{\mathrm{ws}}\) can be interpreted as a measure of how fast water might break through at the new well. The smaller the value of \(d_{\mathrm{ws}}\), the sooner the breakthrough, making the new in-fill well less economical. New well locations and \(d_{\mathrm{ws}}\) values are determined for all model and production rate combinations. For each production rate, the \(d_{\mathrm{ws}}\) values for the homogenized models are then compared to their EDFM counterpart (Fig. 13). We note that the homogenization errors for the SP-EDFM and DP-EDFM models do not change monotonically with respect to the production rate. This behaviour is due to the explicit representation of the fracture corridors. Figure 14 shows the \(d_{\mathrm{ws}}\) measurements for the homogenized DP-EDFM model at different production rates. For a production rate of 0.3 GIP/year (Fig. 14c), the location on the water saturation front that is nearest to the new well is along the y-directional fracture corridor at \(x=50\) m. The homogenization error measure we have used is insensitive to the location of the water saturation front in this fracture corridor. A more robust homogenization error measure would be to run an actual gas production simulation at the new well location and use the resulting breakthrough time for error calculations.

Fig. 12
figure 12

Potential location of a new well (red dot). Circles are defined by the distance from the well to the nearest water saturation front. Water distributions shown are at breakthrough with a production rate of 0.04 GIP/year. a EDFM model, b DP model, c SP-EDFM model, d DP-EDFM model

Full size image
Fig. 13
figure 13

Change in \(d_{\mathrm{ws}}\) for each upscaled model with respect to EDFM results

Full size image
Fig. 14
figure 14

Change in \(d_{\mathrm{ws}}\) for the homogenized DP-EDFM model for different production rates. Saturation fields shown are for diffuse fractures. Production rates are a 0.04 GIP/year, b 0.1 GIP/year, c 0.3 GIP/year, d 1 GIP/year

Full size image

Nevertheless, the proposed homogenization error measure is useful for comparing the performance of the homogenized models relative to the small-scale EDFM model. The results (Fig. 13) show that homogenization error is the largest in the homogenized DP model (between 27 and 83%). On the other hand, both homogenized hybrid models fare better (SP-EDFM below 33% and DP-EDFM below 24%). The improved performance of the homogenized hybrid models (especially the DP-EDFM model) relative to the homogenized DP model shows that the EDFM simulations performed on a small-scale representative model allow us to systematically select an upscaled model with a lower homogenization error.

The results produced by the three homogenized models show that none of the homogenized models are a perfect replacement of the fully resolved model. Among the three, the homogenized DP-EDFM model performs the best as it better captures the shape of the water saturation front as well as the water saturation in the matrix behind the water saturation front. However, the homogenized DP-EDFM model does not perform well in terms of predicting the correct recoveries at breakthrough. Homogenization errors in the DP-EDFM model may be reduced by using a more accurate transfer function that accounts for the transient matrix flow observed in the high-resolution EDFM simulations. The homogenization errors of various continuum and hybrid methods tested here show that in order to make effective engineering decisions, it is important to ensure that an appropriate upscaled model is chosen. We also showed that by using a small-scale EDFM model to characterize flow behaviour in a real naturally fractured reservoir, we were able to better select and construct upscaled simulation models that capture the flow behaviour in the reservoir.

8.3 Effect of Grid Coarsening

The homogenization error studied in Sect. 8.2 only focuses on the effect of substituting heterogeneous permeability fields with homogeneous ones. As discussed in Sect. 3, upscaling error is a combination of both homogenization and discretization errors. The latter is the error arising from grid coarsening. All the simulations in the previous study of homogenization errors have been performed using 1-m grid cell sizes. However, continuum and hybrid models are usually used with much coarser grids. As such, we are also interested in the performance of EDFM, DP, SP-EDFM and DP-EDFM models at coarser grid resolutions. Using the coarsened models generated in Sect. 8.1, production simulations were performed at well rates of 0.04 GIP/year and 0.3 GIP/year. We then determine recoveries at breakthrough and calculate their percentage change relative to the homogenized models (which have 1-m grid cell sizes) (Fig. 15).

Fig. 15
figure 15

Performance of coarsened EDFM, DP, SP-EDFM and DP-EDFM models. Top row corresponds to well rate of 0.04 GIP/year. Bottom row corresponds to well rate of 0.3 GIP/year. a, c Predicted recoveries at breakthrough, b, d percentage change relative to corresponding models with 1-m grid cell size

Full size image

The change in recovery at breakthrough in comparison with the fine grid results shows that EDFM is the most sensitive to grid cell size. The maximum discretization error, equal to the percentage change in recovery at breakthrough, is 45.16% for the 32-m model producing at 0.3 GIP/year. As smaller grid sizes are used, the recovery at breakthrough begins to converge. The DP method is the least sensitive to grid cell size and produces a maximum discretization error of only 3.57%. The SP-EDFM and DP-EDFM methods have maximum discretization errors of 17.95% and 10.53%, respectively. The likely reason for the low discretization error for the DP-EDFM models in comparison with the SP-EDFM model is the use of a transfer function. As grid cell size is coarsened, the transfer function may continue to regulate flow from the diffuse fractures into the matrix. Water then breaks through at the well at around the same time despite the numerical dispersion introduced by grid coarsening. At large grid sizes, numerical dispersion eventually becomes significant enough to cause earlier water breakthrough. On the other hand, the SP-EDFM model assumes that the diffuse fractures and the matrix are always at the same pressure and saturation. As the grid is coarsened, fluid transfer into the matrix may increasingly be overestimated, leading to a delay in water breakthrough. At large grid sizes, numerical dispersion in the explicit fractures eventually becomes significant and begins to shorten the breakthrough time. The DP, SP-EDFM and DP-EDFM models are all more resilient to grid coarsening as their properties were determined through upscaling (Sect. 8.1). The low discretization errors of the coarsened hybrid and continuum models in comparison with the coarsened EDFM models highlight the suitability of continuum or hybrid models for full field simulations using coarse grid cells. Between the simulation methods considered here, although the DP method has the lowest discretization error, it is unsuitable for large-scale simulations due to its large homogenization error. The DP-EDFM method, which has reasonable homogenization and discretization errors, is the most suitable method for large-scale flow simulations. This result reaffirms that the EDFM simulations on small-scale representative models are effective in helping practitioners systematically select full field simulation methods that honour the flow characteristics in the reservoir.

9 Summary and Conclusion

Various simulation methods exist for modelling flow in naturally fractured reservoirs. The selection and construction of full field simulation models are often a subjective process that involves trial and error. In this paper, we suggested that DFM simulation studies can be used to reveal flow characteristics in naturally fractured reservoirs. The improved understanding of flow behaviour allows us to systematically select a simulation method for large-scale flow modelling. This approach was applied to a real fractured gas field using the EDFM method. We first generated an EDFM test model that represents a small sector of the gas field. The EDFM test model was subjected to three different simulation studies (Table 3). The chosen simulation studies are based on commonly studied phenomena that are applicable to gas production. The simulation studies allowed us to improve our understanding of the flow dynamics in the fractured gas reservoir. Using the observed flow behaviours as a guideline, the DP-EDFM method was identified as an appropriate technique for full field simulations.

The applicability of the chosen upscaled method (DP-EDFM) for modelling flow in the fractured gas field was then assessed by comparing it against the EDFM method. Additionally, we also compared the SP-EDFM and DP methods against EDFM. As it is not feasible to run full field simulations with the EDFM method, we cannot compare full field simulation results generated by the upscaled approaches against EDFM. Instead, we adopted the approach of Sablok and Aziz (2008) who separated upscaling errors into homogenization and discretization errors. In our work, the homogenization error was measured by comparing hybrid and continuum models against EDFM models at a similar level of grid refinement. The discretization error is measured by comparing each of the fine grid models against coarse grid versions of themselves. The results showed that DP-EDFM produces lower homogenization and discretization errors when used to model flow in the fractured gas field.

Through this validation procedure, we ascertained that DFM studies on fully resolved small-scale models can reveal useful information about flow behaviour in a naturally fractured reservoir. This information can be used to guide the selection and construction of full field simulation models. By ensuring the full field upscaled models are able to reproduce the observed flow behaviour, we can then ensure that the upscaled model is fit for purpose. As a final note, we would like to re-emphasize that while we have used DFM to characterize flow behaviour in a fractured gas field, the selected simulation studies in Table 3 are not unique to the hydrocarbon industry. As such, the approach can similarly be used in other areas of application like carbon sequestration and geothermal energy production.