## Abstract

Reliable production forecasting for fractured carbonate reservoirs is a challenge. Natural fractures, adverse wettability and complex matrix heterogeneity are all uncertain and can all negatively impact upon recovery. Ideally, we should consider different reservoir concepts encapsulated in a large ensemble of reservoir models to quantify the impact of these and other geological uncertainties on reservoir performance. However, the computational cost of considering many scenarios can be significant, especially for dual porosity/permeability models, rendering robust uncertainty quantification impractical for most asset teams.

Flow diagnostics provide a complement to full-physics simulations for comparing models. Flow diagnostics approximate the dynamic response of the reservoir in seconds. In this paper we describe the extension of flow diagnostics to dual porosity models for naturally fractured reservoirs. Our new diagnostic tools link the advective time of flight in the fractures to the transfer from the matrix, identifying regions where transfer and flux are not in balance leading to poor matrix oil sweep and early breakthrough. Our new diagnostics tools have been applied to a real field case and are shown to compare well with full-physics simulation results.

Carbonate reservoirs contain much of the world's remaining hydrocarbon resources (Burchette 2012), but recovery is typically low. Natural fractures are ubiquitous in carbonates; they provide high permeability pathways and reduce sweep efficiency. Carbonates also tend to be mixed to oil wet which means significant volumes of oil can be left behind in the rock matrix. Fractures, matrix heterogeneity and wettability, and the interactions of these properties, not only reduce the recovery, but these are also primary uncertainties. These features are multi-scale, often sparsely sampled and their distributions can vary significantly well to well and across the reservoir, leading to significant uncertainty in the static and dynamic characterization of the reservoir.

Ideally, we should consider different reservoir concepts encapsulated within a large and diverse ensemble of possible models to provide a robust production forecast capturing the impact of uncertainties (cf. Scheidt & Caers 2002). Full-physics reservoir simulation is still the most common tool used to compare and rank reservoir models. However, for large complex models the simulation time can extend to many hours or even days. Due to time constraints, typically only a small subset of scenarios and sensitivities are considered, which may collapse into a single base case, from which development decisions are made. A single base case model cannot robustly quantify uncertainty and provides little confidence in the decision-making process (cf. Ringrose & Bentley 2015).

Flow diagnostics offer a rapid alternative for screening many models quickly by providing an approximation to the dynamic response of a static model (e.g. Black Oil) (Møyner *et al.* 2015). Flow diagnostics calculate and visualize key dynamic properties of the reservoir and provide preliminary estimates of: the time of flight, drained and swept reservoir volumes, time to breakthrough, sweet spots and well-allocation factors.

The time of flight is the most important quantity obtained in flow diagnostics, it is a useful metric for assessing swept and unswept regions in the reservoir (Møyner *et al.* 2015). The time of flight field is obtained by solving the standard pressure diffusion equation at steady state conditions once and then calculating the fluid velocities for each grid block in a post processing step. The forward time of flight (τ_{f}) is the time it takes an imaginary particle injected at an inflow boundary, such as an injector, to travel to a point, x, in the reservoir and the backward time of flight is the time to travel from point x to a producer. The total time of flight, or residence time, is the sum of the forward and backward time of flight. The time of flight is defined as:*et al.* 2015).

The time of flight is usually obtained from computationally efficient streamline simulations (Batycky *et al.* 1997; Datta-Gupta & King 2007) but more recent tools can compute the time of flight directly on structured and unstructured reservoir model grids (Shahvali *et al.* 2012; Møyner *et al.* 2015; Zhang *et al.* 2017).

Computing tracer partitioning in the reservoir provides further useful insights into reservoir dynamics in the form of swept and drained reservoir regions using a simple majority principle (Møyner *et al.* 2015). The tracer partitioning is determined from a numerical simulation where conservative tracers are injected at the injectors. Under steady-state conditions the tracer concentration is computed as:

Flow diagnostics require an initial solution of the steady state pressure field (with the given well controls) to compute the time of flight, basic quantities and dynamic metrics. This means that the run time for flow diagnostics is almost negligible. Using fast linear solvers, the run time is of the order of seconds, even for multi-million cell models (Stueben *et al.* 2007).

The short run times of flow diagnostics allow many static reservoir models to be screened, compared and ranked, based on their approximated dynamic response. This not only offers a robust and fast method to quantify the impact of geological uncertainty in production forecasting, but also allows us to select a small number of representative reservoir models that capture the full range of uncertainties for calibration, history matching, performance forecasting and robust optimization workflows using traditional full-physics reservoir simulations. Flow diagnostics hence complement, but not replace, traditional reservoir simulation workflows.

This paper is structured in two parts. Firstly, we introduce our development of flow diagnostics for dual porosity models. Here two new metrics will be introduced and compared with full-physics simulations. In the second part of this paper we present a case study of flow diagnostics performed for a real fractured reservoir.

## New flow diagnostics methods for fractured reservoirs

The existing flow diagnostics discussed so far are for single porosity systems only. Single porosity models do not capture important processes for flow and recovery in fractured reservoirs. Dual porosity models are commonly used for modelling flow and recovery for fractured reservoirs (e.g. Bourbiaux 2010; Lemonnier & Bourbiaux 2010*a*, 2010*b*). Fractured rock is represented by a discrete network of matrix blocks, between which there is no flow, which provide the majority of the storage and a flowing medium with little storage for the fractures (Barenblatt *et al.* 1960; Warren & Root 1963). The two separate porous media may exchange fluids via a transfer function, whose form depends on the recovery process. Figure 1 illustrates the dual porosity concept. The matrix–fracture exchange is a primary driver for oil recovery and needs to be accounted for in our diagnostics.

The work presented in this paper extends our earlier work (Spooner *et al.* 2018*a*, 2018*b*) with thorough validation of the flow diagnostics response with full-physics simulation results for a field case. Our new dual porosity flow diagnostics can quickly assess how the matrix–fracture exchange affects reservoir performance. We are particularly concerned with identifying regions of the matrix that are poorly swept either due to poor sweep in the fractures or by poor transfer potential and well performance. We introduce two new metrics for dual porosity models that quantify the effects of transfer and flow. Firstly, we introduce physics-based transfer models which are used to capture the matrix fracture exchange.

### Physics-based transfer models

The transfer functions in most commercial simulators are based on Darcy flow. In these formulations the transfer is computed as a function of the pressure difference between the fracture and matrix (Warren & Root 1963, Kazemi *et al.* 1976; Gilman & Kazemi 1988). Classical transfer functions have limitations, including the lack of modelling of transient effects in the matrix because they assume pseudo-steady-state conditions (Lu *et al.* 2006; Ramirez *et al.* 2007; Abushaikha & Gosselin 2008; Al-Kobaisi *et al.* 2009). Also, the transfer functions take no account of sub-grid-block heterogeneity, which can affect transfer (Di Donato *et al.* 2007; Maier *et al.* 2013).

With classical transfer functions, the transfer must be computed at every time step. Therefore, they are not suitable for flow diagnostics where we only solve the governing equations once to obtain the basic dynamic quantities. To overcome this limitation, we employ physics-based transfer functions which provide the transfer rates without the need for transient simulation. Physics-based transfer models provide a transfer rate as a function of the saturation difference between the fractures and matrix (Kazemi *et al.* 1992). The transfer rate constant, herein referred to as β, models the speed at which fluids are exchanged between fractures and matrix for a given recovery process (e.g. imbibition or gravity drainage). For any given petrophysical and fluid data, β can be computed directly based on analytical solutions which capture both early- and late-time behaviour during fracture–matrix fluid exchange, specifically imbibition (Schmid & Geiger 2012, 2013; Schmid *et al.* 2016) and gravity drainage (March *et al.* 2016, 2017, 2018). Figure 1 illustrates the difference sin formulation between classical and physics-based transfer models. We note that these transfer functions are not limited to dual porosity models but can also be applied, in principle, to dual porosity–dual permeability models, i.e. cases where fluid flow in the matrix is not negligible. However, extending the work to dual porosity–dual permeability models is beyond the scope of this paper.

In the discussion, we assume that recovery from the matrix is dominated by imbibition in our flow calculations, but we note that other recovery mechanisms such as gravity drainage can be included through the appropriate choice of β. For imbibition, β can be used in Aronofsky's model (Aronofsky *et al.* 1958) to estimate the recovery R at time t as a function of the ultimate recovery R_{∞} and β*et al.* 1976) rather than the characteristic matrix block length in our implementation of the transfer rate formulation as*et al.* 2017, 2018), have been implemented in the open source Matlab Reservoir Simulation Toolbox (MRST) (Lie 2014).

### Metric 1: dual porosity Damköhler numbers

Damköhler numbers (Da) are dimensionless numbers used in chemical engineering to relate the chemical reaction timescale to the timescale of flow or convection (Fogler 1999). The Da can be used as a metric for conversion efficiency. In a reaction cell, the most efficient conversion is achieved if the reactant is produced and removed at a rate which maintains the optimum concentration of reactant. A Da of 1 occurs when the reaction rate and transport rate are equal.

The Da principle is applied to the dual porosity model by taking each grid cell to be equivalent to a reaction vessel producing oil into the fracture and an advective flux which removes the oil from the grid cell location. We define the dual porosity Damköhler number (DaDP) as logarithm (base 10) of the ratio between the transfer flux from the rock matrix and the advective transport flux within the flowing fractures

#### Comparison with simulation for a simple case

The following simple 2D example illustrates the new dual porosity flow diagnostics method and how it compares with full-physics simulations. We take an inverted 5-spot box model consisting of 441 cells (21 × 21 × 1), with a water injector in the centre and a producer at each corner. The model has homogenous porosity and permeability in the matrix and the fracture distribution is assumed to be homogenous with a constant shape factor of 1 m^{−2}. The model features four different rock types based on wettability. Two realizations of this model are considered; each has the same number of cells of each rock type but these are distributed differently in each case. For Case A the rock types are correlated and grouped together and for Case B the rock types are distributed randomly (Fig. 2). A water flood was simulated for both realizations of the 5-spot model for 5000 days in MRST. Due to the 2D nature of the reservoir, gravity effects are negligible and imbibition dominates fracture–matrix fluid exchange. All model parameters and rock fluid data are given in Appendix 1.

The transfer rate coefficient β was computed for each rock type via equation (4). The transfer rate coefficient for each rock type is given in Table 1.

Having computed the values for β, we can then compute basic flow diagnostic quantities and the first dual porosity metric, DaDP. Figure 3 shows the DaDP values for each case and the average DaDP values in each drained region. Figure 3 also shows the simulated saturation of oil in the matrix and the fractures for each case.

First, we examine the DaDP values. In each case, the DaDP values are positive for cells with water wet rocks, and negative where the rock is mixed wet (Fig. 3). The average advective flux in our 5-spot model is 0.0016 m^{−2} s^{−1}. In terms of orders of magnitude, we have for the mixed-wet rocks:^{−7} s^{−1} (two orders of magnitude greater than for the mixed wet rock types), hence the DaDP values are all positive in cells of rock types 2, 3 and 4.

The distribution of the rock types also impacts the average DaDP values. For Case A the avDaDP values for three wells (W1,W3 and W4) are positive, indicating that on average the DaDP regime is optimal. However, for W3, which is located in a region with predominantly mixed wet rock, the avDaDP value is negative. In contrast, for Case B the avDaDP values are all positive with a lower variance between the average DaDP values compared to Case A. We will shortly highlight the impact of this result on well performance.

Next, we compare the DaDP values with the simulated saturation of oil in the matrix. In the cells with negative DaDP, the saturation of oil remaining in the matrix after 1410 days of production is still high due to little transfer. This is most clear for Case A where we see large regions of oil remaining in the matrix.

A key observation from the fracture saturations is the striking difference in the shape of the flood fronts (Fig. 3d). For Case B we observe a uniform front but, for Case A the flood front progresses quickly through the regions with negative DaDP values. In this case, the transfer from the matrix with very low β values is too low to replenish the fractures adequately to delay water breakthrough.

For Case B the random nature of the high quality and low quality rock types essentially homogenizes the transfer with respect to the inter-well distance. This essentially averages the transfer resulting in a system that behaves like a single-rate system despite containing multiple rock types with different transfer rates.

Differences in the flood fronts significantly affects the predicted water breakthrough (Fig. 4). For Case, A the water cut profiles are very different for the individual producers and the range of water breakthrough times are large (>1200 days). The values of avDaDP correspond well to the simulated water breakthrough. Breakthrough occurs first at Well 2, which has a negative avDaDP. Breakthrough at Well 3 is early compared to the more highly performing wells. The avDaDP for Well 3 has a small positive value (0.08). Although this value is positive we should note that the transfer flux used to compute this value is the maximum possible value assuming the fracture block is fully water saturated. However, actual transfer rates may be much lower depending on the actual saturation in the fractures. Therefore, small positive DaDP values can also indicate poor performance compared to other wells.

Wells 1 and 4 both have positive DaDP values greater than one. For these wells, the water breakthrough occurs much later at over 1000 days. In each of these wells’ locations, the rock matrix is dominantly water wet with much greater transfer rates allowing adequate replenishment of the fractures for longer.

For Case B the water cut profiles for all four wells are similar (Fig. 4). Here the avDaDP values are all positive and the variance is low (0.05) compared to Case A (0.65). The variance in avDaDP is therefore also useful in determining if wells in a field are likely to behave similarly or not.

This example has shown that very low transfer rates can significantly impact flood fronts in the fractures and water breakthrough as well as matrix oil drainage. Our DaDP metric quantifies ‘low transfer’ in the context of the flow regime. Hence DaDP and avDaDP can be used to assess flood patterns, drained and swept regions of the matrix, and possible well performance.

### Metric 2: dual porosity F–Γ curves

The second dual porosity flow diagnostic we introduce is the F–Γ curve (Spooner *et al.* 2018*b*). This metric is based on the same principles as the dynamic F–Φ curve proposed by Shook & Mitchell (2009). Dynamic F–Φ curves offer an extension to the familiar static F–Φ curves which are based on porosity thickness and permeability thickness (Shook & Mitchell 2009). Dynamic F–Φ curves are instead based on swept pore volumes computed from the velocity field and hence capture the reservoir dynamics (Shook & Mitchell 2009). Consequently, they respond to any factors affecting the flow such as stagnant regions and well placement. F and Φ are defined based on Shook & Mitchell (2009) as_{j(i)} is the advective flux and Φ_{j(i)} is the pore volume in the j^{th} (i^{th}) cell. The dynamic Lorentz coefficient (L_{c}) can then be computed from the F–Φ data as

We have used the principle of the dynamic F–Φ curves to compare the heterogeneity of the transfer rate compared with the heterogeneity of the flow field in the fractures. We compare the cumulative transfer flux (Γ) to F. We define the cumulative transfer flux Γ as

When the transfer rate distribution is proportional to the distribution of pore volumes the F–Γ curve will match the F–Φ curve. This is shown for a simple realization of our 5-spot model where only a single rock type is present (Fig. 5), i.e. this example shows single-rate behaviour with a uniform β. This uniform transfer behaviour is favourable since we can optimize the sweep efficiency based on the single porosity flow diagnostics alone.

For Case A the F–Γ curve (shown in yellow) is closer to ideal but further from the F–Φ curve than Case B. The similarity between the F–Γ and F–Φ curves can be quantified by the root mean square error (RMSE). The RMSE for Case B is low. This result agrees with the simulated behaviour we showed previously in that Case B behaved like a single-rate or uniform system, whereas Case A showed non-uniform transfer behaviour across the field with regions of different wettability impacting the sweep fronts.

We have introduced two new metrics for dual porosity systems, DaDP values and dynamic F–Γ curves. The DaDP value combines the advective flux and physics-based transfer to identify regions where flow and transfer are not in balance, which can lead to poor matrix sweep and poor sweep fronts. The avDaDP metric provides a quantification of the DaDP regime affecting the drained volume for each well. Negative avDaDP can indicate that a well may experience early water breakthrough. The variance in the avDaDP metric indicates whether we can expect similar or widely different behaviour between the wells in a field.

Our second metric, the F–Γ curve, and the resulting Lorentz coefficient can be compared with the dynamic F–Φ curve to assess and quantify whether the system is likely to show multi-rate or single-rate behaviour.

Next, we will apply our newly extended flow diagnostic tools to a real field case study. We will assess which of the flow diagnostic metrics agree well with the simulation results and which offer potential to be used for screening and ranking models.

## Teapot Dome case study

We now present a case study comparing the flow diagnostic response to the full physics simulation results for a real fractured reservoir. The Teapot Dome is a naturally fractured reservoir located in Wyoming, USA, and it has been operated by Stranded Oil Resources Corporation since 2015 (US Department of Energy 2015). A previous study noted that despite relatively good data about the fractures in the Teapot Dome there is still a high degree of uncertainty in the distribution of fracture intensity between the wells (Ahmed Elfeel *et al.* 2013). The correlation of the fracture intensity between the wells affects the flow paths at the field scale and hence the flooding patterns and well connectivity. The simulation run-time for the model is around 2 hours, hence simulating the number of models required to quantify the impact of uncertainty would require a very large time investment. We demonstrate the application of dual porosity flow diagnostics using the Teapot Dome model and validate the flow diagnostics response with full-physics simulations.

### Background

The Teapot Dome is an elongated north–south-trending anticline located at the SW edge of the Powder River Basin. The geological and simulation models used in our study were originally built to examine the impact of DFN upscaling methods on history matching and production forecasting (Ahmed Elfeel *et al.* 2013). The Pennsylvanian Tensleep formation forms the main reservoir in the Teapot Dome. The formation consists of aeolian sands, interbedded with sabkha and marine dolomites (Ouenes *et al.* 2010). The reservoir can be divided into three main units. The uppermost unit (Unit 1) is fractured sandstone, Unit 2 consists of fractured dolomites and Unit 3 is a sandstone (Fig. 6). The rock matrix is hard, tight and water wet. The mean matrix porosity is 7% and the horizontal matrix permeability is 14 mD. A single relative permeability and capillary pressure curve (provided in Appendix 2) describes the wettability of the rock matrix.

The fracture data is for the Teapot Dome reservoir include borehole images, seismic data and outcrop data. Outcrop fracture analogues were described by Cooper *et al.* (2006). They found fracture sets that were parallel, perpendicular and oblique to the fold hinge. Schwartz (2006) analysed borehole images and provided subsurface strike dips and aperture data. The DFN for the model used in this study was built by Ahmed Elfeel *et al.* (2013) using input data described by Cooper *et al.* (2006) and Schwartz (2006). The fracture intensity is constrained to well log data and co-kriged with the grid curvature. The fracture intensity is populated via sequential indicator simulation on the grid with a correlation length of 2450 m in both the major and minor directions. The correlation length was highlighted as an uncertainty in previous work (Ahmed Elfeel *et al.* 2013).

The simulation model contains 267 294 active dual porosity cells. There are ten producers and three water injectors, all of which are bottom hole pressure (BHP) controlled. We have used CMG's IMEX Black Oil reservoir simulator to simulate the production for a water flood. For the base case model, the simulation time for 3500 days production is around 2 hours on a standard workstation.

Teapot Dome experiences steep water production profiles (Fig. 7). In terms of time to water breakthrough the best performing well is 12_76. This well and other wells where breakthrough occurs later are furthest from the water injectors. A cluster of wells (5_62, 1_43-2, 2_63 and 13_67-1) show breakthrough after 1 year of production. The cumulative oil production varies significantly between wells (Fig. 7). It is not clear simply by examining the static models why there should be such large differences between the wells' productivity.

### Teapot Dome flow diagnostics

The flow diagnostics response and the simulation results for the base case are compared below. The run time for the flow diagnostic is *c.* 6 s on a standard workstation. Note that we have used an aggregation-based algebraic multigrid method solver AGMG (Notay 2010) to accelerate the incompressible pressure solution in MRST.

#### Computing transfer rates

The first step is to compute the physics-based transfer rate coefficients. Under water flooding we consider spontaneous imbibition as the dominant recovery process since the matrix is water wet and used equation (4) to compute the transfer rate coefficients. This solution of equation (4) requires an iteration for each unique set of input values (porosity, permeability, wettability data), each iteration takes around 10 s.

The matrix wettability for the Teapot Dome is described by a single rock type, but the porosity and permeability are continuous which means there are thousands of possible combinations of porosity and permeability values. Computing the transfer rates for so many unique sets of input parameters to equation (4) would take an impractical amount of time. Instead, we truncate the porosity and permeability to discrete values that are chosen to ensure good representation of the original distribution. The shape factor is applied as a multiplier after the iterative process used to solve equation (4); therefore, we can use the original shape factor values directly. The resulting β values are shown for the base case in Figure 8. The highest values of β are predicted to be in the fractured dolomites (Unit 2). Higher values are also observed in the anticlines and synclines in the upper sandstone unit.

As a first diagnostic test we can assess the distribution of the transfer rate coefficients and their impact on the ideal recovery predicted via the Aronofsky model (Equation 3). Figure 8b shows a combined histogram and recovery curves for the transfer rate distribution. The Aronofsky model (Equation 3) was used to compute the normalized recovery for different orders of magnitude of β. The line thickness and intensity indicate the proportion of cells with a transfer rate of that order. In this case most of the cells have transfer rates of the order of 10^{−8} to 10^{−6} m^{−2} s^{−2}. Under the assumptions of the Aronofsky model this represents a block recovery of 3– 30 years to maximum recovery. This implies that as long as the fractures are well swept, recovery from the matrix is on typical field life time scales. However, the actual recovery from the blocks will be slower as the Aronofsky model does not consider saturation differences.

#### Time of flight

The time of flight is computed in the flowing fractures. For realistic reservoir models it is often useful to view the log (base 10) of the total time of flight to exaggerate the swept and unswept regions. Figure 9a shows the time of flight distribution for the Teapot Dome. The regions coloured in bright yellow, such as the downthrown fault block and the outskirts of the field, have the highest residence time, i.e. these regions take much longer to sweep. The highest residence times are 300 years or more, so these regions would normally be considered unswept.

#### Dynamic F–Φ and F–Γ plots

Having computed the time of flight we can examine the F–Φ and F–Γ plots (Fig. 9). The F–Φ curve lies skewed towards the Y axis, indicating that overall displacement is heterogeneous. However, this includes cells that are poorly swept in the downthrown fault block. Improved sweep of this region will move the F–Φ curve towards more homogeneous displacement.

The similarity between the F–Γ and F–Φ curves indicates that the swept volumes are in balance with the distribution of transfer. As shown for the simple 5-spot example, when the transfer and sweep are well balanced the system behaves as a single-rate system. In a single-rate system the overall transfer effects are uniform and the system behaves as though a single averaged transfer rate is acting. In these cases, the advective flow dominates the dynamic response. Therefore, for the Teapot Dome reservoir, the single porosity flow diagnostics are expected to provide adequate approximations for the dynamic behaviour.

#### Dual porosity Damköhler number

DaDP is high in poorly swept regions, such as the downthrown block, and lower closer to the wells where the flow is faster (Fig. 10a). All of the average DaDP values for each drained volume are positive (Fig. 10b). Higher values of DaDP (>3) are associated with the unswept regions and are not included in avDaDP. The variance between the avDaDP values is low indicating that in terms of flux and transfer the wells are likely to behave fairly similarly.

#### Drained and swept volumes

We noted previously that the wells closer to the injector breakthrough earlier than the wells furthest from the injectors (Fig. 7). We can approximate the breakthrough time using the tracer thresholding utility within the main flow diagnostics toolbox. This allows us to visualize swept volumes and the location of the flood fronts for each injector at a given time.

Figure 11 shows the flooded reservoir volumes for each injector after 365 days of water flooding. There are some notable differences in the flood fronts between the flow diagnostics and full-physics simulation. In the simulation the flood fronts preferentially travel eastwards down dip whereas the flow diagnostics flood fronts spread more uniformly. The difference can be attributed to the fact that flow diagnostics do not account for gravity in the advective flow. The greater the impact of gravity (i.e. the greater the reservoir thickness and the larger the vertical permeability and fluid density contrast) the greater the discrepancy between diagnostic results and simulation.

The Teapot Dome is a relatively thin reservoir (36 ft thickness), therefore we can use flow diagnostics with the caveat that we expect some degree of error due to gravity effects. The differences impact our method of averaging DaDP within each drained volume and the time at which the diagnostics predicts breakthrough.

We have approximated the water breakthrough time for each well by examining the location of the flood fronts at different times and compared this with the simulated breakthrough time. We have grouped the results by wells that breakthrough within six-month intervals. The ranking determined by the flow diagnostics and simulation show good agreement (Fig. 12a). Both methods identify the same most poorly performing wells and the same highly performing wells in terms of breakthrough time.

We can also use the tracer thresholding utility to estimate drained volumes at a given time. We have computed the drained regions at the end of production for the base case model and computed the reservoir pore volume (fracture and matrix volume) of these regions. We have ranked the wells by drained volume and compared this to the cumulative oil production ranking from the simulation. The overall ranking for drained volume and cumulative oil production shows excellent agreement (Fig. 12b).

The diagnostics are effective in determining the best and worst performing wells in terms of oil production. There are some differences in the ranking for the mid performing wells (5_62, 2_63, 1_43-2 13_76-1 and 12_76). We attribute these differences to the differences in the predicted flood fronts (Fig. 11). The flow diagnostic calculation does not account for gravitational effects in the of the flow field. This simplification in the physics manifests as differences in water breakthrough times and drained volumes, particularly for wells that are close to injectors. Wells 1_43-2 and 5_62 show these different behaviours clearly. These wells are at a similar distance from injector 3_53 but located either side of the crest with well 5_62 on the more steeply dipping side. The flow diagnostics indicate that the breakthrough is earlier and cumulative production is lower for well 1_43-2 than for well 5_62. The simulation accounts for the stronger gravity effect for the higher dip for well 5_62 and hence predicts earlier breakthrough with a higher oil recovery.

The dual porosity metrics indicate that the dynamic behaviour of the Teapot Dome system is dominated by the advective flux in the fractures. The transfer rates throughout the field are fast enough to support adequate replenishment of the fractures and there are no large areas with low transfer rates that could lead to earlier breakthrough than expected. The lack of strong transfer effects means that water breakthrough is linked more strongly to the proximity of the producers to the injectors. This behaviour is precisely what we would expect for a reservoir with simple and favourable wettability.

Despite some differences caused by modelling simplified physics, the flow diagnostics yields an overall similar ranking for both the water breakthrough times and produced oil compared to the full-physics simulation. Flow diagnostics not only agrees well with the full-physics simulation results but importantly provides these results in a fraction of the time.

## Conclusion

Flow diagnostics offer a fast approach for comparing geological models based on their approximated dynamic response. In a similar manner to streamline simulation, flow diagnostics sacrifice some level of physical detail in exchange for fast run times. The key benefit is that flow diagnostics provide key information about the reservoir dynamics under different production mechanisms in seconds.

Flow diagnostics are not a replacement for reservoir simulation but offer a natural pre-processing step that complements modern reservoir simulation, uncertainty quantification, and robust optimization workflows, especially for naturally fractured reservoirs where fracture conductivities are notoriously difficult to constrain but are a first-order control on reservoir performance. The rapid run time makes a multi-realization approach to forecasting possible. A diverse subset of geological models can be chosen based on the flow diagnostic response. These models can be subsequently history matched in a Bayesian framework to available production data to produce a set of geologically diverse but plausible reservoir models for forecasting with more robust uncertainty envelopes. Flow diagnostics therefore allows the user to consider a broad range of uncertainties whilst reducing the number of models required for forecasting using computationally costly, but physically more accurate, full-physics simulations.

Our new dual porosity flow diagnostics identify and quantify the interplay of fracture–matrix transfer and advective flow, providing insight into the interaction between the rock matrix and the fractures, which is a key uncertainty in naturally fractured formations. The DaDP metric offers new understanding of the competition between fracture flow and fracture–matrix transfer processes in a dual porosity system. DaDP hence provides a quantitative measure to explain poor matrix sweep, identifies sites for infill drilling and optimization of well controls. The F–Γ curve can be combined with the dynamic F–Φ curve to assess the relative importance of both the range and distribution of transfer rates.

We have applied dual porosity flow diagnostics to a real fractured reservoir, the Teapot Dome in Wyoming, USA. The dual porosity flow diagnostics show that, due to the simple rock type distribution, the overall impact of the transfer is fairly uniform across the field. Hence there are no regions indicated to be poorly swept or at risk of early breakthrough due to poor replenishment from the matrix. Instead advective flow is shown to dominate the behaviour in the system. The oil production and breakthrough time in each well can be predicted using the single porosity flow diagnostics.

Flow diagnostics offer geoscientists a tool to quickly assess the reservoir dynamics directly from the geomodel. This provides a unique opportunity to understand the flow behaviour of the system early and explore the impact of different geological scenarios prior to full flow simulation, helping geoscientists and engineers select the most important features and cases to consider for forecasting.

## Acknowledgements

The authors would like to thank Energi Simulation and EPSRC for their financial support. All data for this paper are properly cited and referred to in the reference list. Simulations were carried out with the open source code MRST from SINTEF available from https://sintef.no/projectweb/mrst/, and the related Matlab scripts are available from http://carbonates.hw.ac.uk

## Funding

This work was funded by Energi Simulation and The Engineering and Physical Sciences Research Council is a British Research Council.

## Appendix 1: Supplementary information for the 5-spot model

**Table A1.** *Parameters for the 5-spot box model*

Parameter | Value | |
---|---|---|

Dimensions | 2000 m × 2000 m × 20 m | |

Grid cells | 21 × 21 × 1 | |

Matrix properties | ||

Porosity (−) | 0.2 | |

Permeability (mD) | 10 | |

Fracture properties | ||

Porosity (−) | 0.01 | |

Permeability (mD) | 500 |

## Appendix 2: Supplementary information for the Teapot Dome

- Received October 2, 2018.
- Revision received May 9, 2019.
- Accepted May 9, 2019.

- © 2019 The Author(s). Published by The Geological Society of London for GSL and EAGE. All rights reserved