Dual‐Porosity Flow Diagnostics for Spontaneous Imbibition in Naturally Fractured Reservoirs

We extend single‐porosity flow diagnostics to dual‐porosity systems using a novel retardation factor R to account for the effect of fracture‐matrix transfer on breakthrough times and displacement efficiency during two‐phase flow in fractured reservoirs. R is based on an analytical solution for capillary‐driven fluid exchange between the fractures and rock matrix. By linearizing R the time‐of‐flight τ⋆ is adjusted to include fracture‐matrix transfer and derive new metrics, the dynamic Lorenz coefficient Lc⋆ to quantify the dynamic heterogeneity, and the dual‐porosity sweep efficiency Ev⋆ to estimate how efficiently the injected fluid displaces the reservoir fluid over time. We have tested different formulations of R across three case studies with increasing complexity to analyze the applicability and limitations of dual‐porosity flow diagnostics. This analysis reveals that as long as flow in the fractures is faster than fracture‐matrix transfer, dual‐porosity flow diagnostics provide useful approximations when assessing displacement efficiencies and identifying the wells that are at most and least likely to experience early breakthrough. We show that Lc⋆ and Ev⋆ can be combined with stochastic optimization algorithms to improve the displacement efficiency in a 3D reservoir case study. Since a single dual‐porosity flow diagnostics calculation requires less than 1 min while a full‐physics simulation takes 2 h, we can now quickly screen a large parameter space to identify scenarios that need to be studied in more detail using full‐physics simulations. Hence, our new dual‐porosity flow diagnostics complement and accelerate state‐of‐the‐art uncertainty quantification and optimization workflows for fractured reservoirs.

Flow diagnostics are a well-established technique for screening the ensemble and accelerating uncertainty quantification studies. Flow diagnostics perform computationally efficient numerical tests on each ensemble member to obtain qualitative and visual information about the approximate flow behaviors in a reservoir in the form of fluid flow paths, swept and drained reservoir volumes, and other quantitative metrics of the heterogeneity such as Lorenz coefficients or breakthrough times (Datta-Gupta & King, 2007;Krogstad et al., 2017;Møyner et al., 2014;Shahvali et al., 2012;Shook & Mitchell, 2009;Watson et al., 2020). Flow diagnostics work particularly well in situations where the reservoir dynamics are dominated by pressure gradients between wells and by geological heterogeneity, and when the nonlinear effects of multi-phase flow are small. They achieve faster computation times compared to full-physics simulations through simplified reservoir physics that consider an incompressible, steady-state system which usually contains only a single phase. Flow diagnostics utilize the time-of-flight, which describes the travel time it takes an arbitrary particle to travel a certain distance through the reservoir (Batycky et al., 1997;Crane & Blunt, 1999;Dagan, 1984;Datta-Gupta & King, 2007;Di Donato & Blunt, 2004;Lie, 2019;Pollock, 1988). Because flow diagnostics are significantly faster than full-physics simulations, they allow more exploration of the parameter space. Since flow diagnostics also capture the reservoir dynamics, we can apply clustering (e.g., Arnold et al., 2016;Scheidt & Caers, 2009) to group models based on flow to decide on the subset of models for full-physics simulation.
Flow diagnostics are well-established for single-porosity reservoirs, that is, reservoirs where fractures are absent or only enhance the reservoir permeability. However, reservoirs that contain well-connected fracture networks exhibit a dual-porosity, or mobile-immobile, behavior with two distinct time-scales and flow behaviors: fast fluid flow in the fractures which is dominated by advection and slow fluid flow in the rock matrix which is dominated by diffusion (in the absence of a second fluid phase) or capillary and gravitational forces (in the presence of a second fluid phase) (e.g., Al-Kobiasi et al., 2009;Berkowitz et al., 2001;Haggerty et al., 2001;Kazemi & Gilman, 1993;Neretnieks, 1983;Ramirez et al., 2009;Warren & Root, 1963).
The aim of this study is to build upon our earlier work (Spooner et al., 2019) and develop a dual-porosity flow diagnostics framework for naturally fractured reservoirs with well-connected fracture networks. To achieve this aim, we make use of extensive research from the groundwater literature which quantifies the exchange of solutes between mobile and immobile domains (e.g., Berkowitz et al., 2001;Brusseau et al., 1989;Cvetkovic & Dagan, 1994;Dagan, 1984;Haggerty et al., 2001;Maghrebi et al., 2013;Sanchez-Vila et al., 2007;Willmann et al., 2008) as well as relatively new semi-analytical solutions that allow us to model the exchange of fluid phases between fractures and rock matrix due to capillary forces (i.e., during spontaneous imbibition) and graviational forces (i.e., during gravity drainage) (March et al., 2016(March et al., , 2018Schmid & Geiger, 2012Schmid et al., 2016). This approach allows us to reformulate the time-of-flight and obtain new metrics that enable the computationally efficient approximation of the rate at which injected fluids (e.g., water in a hydrocarbon reservoir or CO 2 in a saline aquifer) migrate through fractured reservoirs. The availability of such metrics is crucial for the subsequent comparison and ranking of the individual ensemble members that capture the inherent geological uncertainties in a naturally fractured reservoir model.
The study is organized as follows: In Section 2, we develop the methodology by first providing an overview of the time-of-flight calculations and then discuss how these calculations can be extended to dual-porosity systems under certain assumptions; in this section we also introduce existing and new metrics to quantify flow processes in dual-porosity systems using flow diagnostics. In Section 3, we present three applications of increasing complexity for the dual-porosity flow diagnostics: we start with a simple one-dimensional model and a simple two-dimensional model to illustrate the key features and limitations of the dual-porosity flow diagnostics, before applying the flow diagnostics to a full 3D model that is based on real subsurface data. We close with a simple optimization study that demonstrate how the flow diagnostics could be used to accelerate calculations that help us to identify improved well controls in the 3D reservoir model.

Methodology
The aim of this section is to introduce the fundamental concepts of dual-porosity flow diagnostics by expanding the classical time-of-flight calculations for single-porosity media to dual-porosity systems. We develop the methodology by introducing a retardation factor R, which quantifies the interactions between the SPOONER ET AL. 10.1029/2020WR027775 2 of 21 fractures and rock matrix. By including R in the time-of-flight calculations, we can then estimate how the breakthrough of the injected fluid traveling through the mobile fractures is delayed due to its interactions with the immobile rock matrix. More specifically, R models the capillary-driven exchange of two immiscible and incompressible fluid phases between the fracture and matrix. This retardation process is linearized by comparing, in each grid block, the rate at which the fluids are exchanged between the two domains with the rate at which the fluids are transported in the fracture.

Brief Review of the Time-of-Flight Calculation
The time-of-flight τ is defined as the time needed by an imaginary particle to travel a certain distance s along a streamline with total length L (Pollock, 1988 where  u is the Darcy velocity along the streamline and φ is the porosity. The velocity field in the reservoir is obtained by solving the standard steady-state pressure equation and Darcy's law in the absence of gravitational forces where k is the permeability tensor, μ the fluid viscosity, q a source term, and p the fluid pressure.
The forward time-of-flight τ f is defined as 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; similarly, the backward time-of-flight τ b is the time to travel from point  x to a producer. The residence time ˆ is the sum of τ f and τ b .
Injected particles flowing through a reservoir take different paths to the producer; some particles may take a shorter or faster path while others take a slower or longer path. τ is therefore a distribution of travel times with the average τ commonly used to assess the swept and unswept regions of the reservoir (Møyner et al., 2014;Shahvali et al., 2012). The first arrival time is τ corresponding to the fastest moving particle, which is a useful measure to assess and predict breakthrough times for injected fluids.
Solving Equation 2 in discretized form (e.g., using a finite difference method) and with the appropriate boundary conditions yields the velocity field in the reservoir based on which the streamlines can be traced (Batycky et al., 1997;Datta-Gupta & King, 2007). Streamlines therefore characterize the impact of heterogeneity in porosity and permeability for a given set of flow boundary conditions, but they do not capture changes in flow due to changes in well rates, well locations, or of all reservoir mechanisms such as viscous cross-flow or gravity slumping. It is important to note that while the velocity  u varies along a streamline, the velocity is locally constant within the a given grid block i that the streamline traverses at a distance s i (Figure 1).
Traditionally, τ is computed along streamlines (Batycky et al., 1997;Datta-Gupta & King, 2007) and represents a spatial coordinate, that is, τ varies as a function of s ≤ L. It is also possible to compute τ directly on the grid of the reservoir model without the additional step of routing streamlines through the model, which SPOONER ET AL.
10.1029/2020WR027775 3 of 21 Figure 1. Illustration of an imaginary particle (denoted by the black dot) traveling along a streamline with length L between an injector and producer. The streamline passes through a series of grid blocks i, j, and k, and the particle travels the corresponding distances s i , s j , and s k , respectively, when traversing these grid blocks. The porosity φ and Darcy velocity  u varies along the streamline but is locally constant and uniform in each grid block. Likewise, concentration c, saturation S, retardation factor R, and Damköhler number  vary along the streamline but are locally constant within each grid block. The time t the particle spends in gridblock i is given by , and the time-of-flight τ that the particle experiences by traveling a certain distance s ≤ L along the streamline is given by Equation 1. The breakthrough is delayed if the particle interacts with the rock, for example diffuses from a mobile domain (indicated by the streamline path) into an immobile domain (e.g., the rock matrix, indicated by the dashed line), because the particle experiences retardation R i ≥ 1 in grid block i. As a result, the increased travel time t ⋆ through grid block i is given by and the increased time-of-flight τ ⋆ , which the retarded particle experiences when traveling a distance s ≤ L, is given by Equation 7. yields volume-averaged values of τ, that is, the average values of τ for all streamline segments that cross a given grid block. The grid-based values for τ hence present the propagation of a piston-like displacement front in a reservoir and can be obtained by solving the following equation using the finite volume method (Eikemo et al., 2009;Natvig et al., 2007) with τ f | inflow = 0 and τ b | outflow = 0.

Calculating the Time-of-Flight for Mobile-Immobile Domains
Having reviewed the fundamental concepts of the time-of-flight calculation, we now expand it to a dual-porosity system by introducing the concept of the retardation factor R to this calculation. It is not uncommon that particles that flow through a geological formation are retarded because they absorb to or react with the mineral surface, or are temporarily trapped in regions with much slower transport velocities. This exchange of particles between the mobile and immobile domains causes a retention of the particles, that is, the average breakthrough time increases. A wide range of models have been proposed to quantify the rate at which particles move through a reservoir when there is additional interaction between the mobile and immobile domains in a geological formation (e.g., Berkowitz et al., 2001;Brusseau et al., 1989;Cvetkovic & Dagan, 1994;Dagan, 1984;Di Donato & Blunt, 2004;Haggerty et al., 2001;Sanchez-Vila et al., 2007;Willmann et al., 2008).
To develop a methodology that allows us to estimate the delay in τ caused by the retardation of fluids in a dual-porosity system due to the interactions between the mobile fractures and immobile rock matrix driven by capillary and/or gravitational forces, we first start by using simple single-phase transport analogy before we arrive at the two-phase flow process. In this analogy the injected wetting phase can be imagined as an adsorbing solute where the capillary-or gravity-driven exchange between fractures and matrix is akin to an adsorption process. Assuming incompressible single-phase flow in 1D along a streamline, the general equation for the kinetic transport of a solute of concentration c is given by where v = u/φ and x ≡ s is the distance along a streamline in 1D. In the following, we assume that the dispersion coefficient D is zero. The adsorption term    Γ t describes the kinetics of adsorption onto the rock, where Γ is the mass of solute or volume of solute on the rock and ɛ is a factor depending on the units. Equation 4 is closed by using the appropriate model for Γ. At equilibrium conditions, such a model could be the well-known Freundlich isotherm Γ eq = K ⋅ c 1/n where K is the Freundlich constant and 1/n quantifies the intensity of the adsorption. A possible kinetic equation where α 1 and α 2 are coefficients to quantify the rate of adsorption.
inserting Equation 5 with D = 0 into Equation 4, and rearranging yields Recall that c, v, and R can vary spatially, that is, along a streamline or from grid block to grid block, but are always locally uniform within a given grid block ( Figure 1). The retarded time-of-flight τ ⋆ that an imaginary particle experiences while it flows through the mobile fractures and interacts with the immobile rock matrix is given by Cvetkovic and Dagan (1994), Maghrebi et al. (2013), and Villermaux (1974

Computing the Matrix-Fracture Transfer
Having expanded the time-of-flight calculation to account for the fracture-matrix interactions that an imaginary particle might experience, we now need to calculate τ ⋆ for situations where two immiscible fluids flow in a dual-porosity system. We consider incompressible and immiscible two-phase fluid flow where capillary and gravitational forces in the mobile domain are negligibly small. In analogy to Equation 4 with D = 0, the change in water saturation S w along a streamline in the mobile domain of a dual-porosity system is given by Di Donato and Blunt (2004) where v t = v w + v nw is the total fluid velocity of the wetting phase w and nonwetting phase nw and    w is the fractional flow function in the fracture, with k r being the relative permeability. Superscripts m and f denote matrix and fracture properties, respectively. Since k r for a fracture is usually approximated by k rw = S w and k rnw = S nw , we obtain piston-like displacement with v ≈ v t if μ w ≈ μ nw . Note that the approximations obtained by the flow diagnostics become less accurate for cases where μ w and μ nw as will be discussed in the example applications below. With these definitions, we obtain the following identities between Equations 4 and 9: As in the previous sections, S, v t , and  can vary along a streamline or from grid block to grid block but are locally uniform within a given grid block.
The transfer function  defines the exchange of fluids between the mobile and immobile domain, respectively the fracture and rock matrix, as a function of fluid and matrix properties. Dropping the index w for simplicity,  is given by Di Donato et al. (2007) where β is a coefficient describing the rate of fluid exchange due to spontaneous imbibition (March et al., 2016;Schmid et al., 2016;Schmid & Geiger, 2012 and S m⋆ is the maximum saturation the wetting phase can reach in the matrix; for spontaneous imbibition into a water-wet rock matrix we have S m⋆ = 1 − S nwr with r denoting the residual saturation. Equation 10 has the same functional form as the widely used transfer function of Kazemi and Gilman (1993). In contrast to these existing transfer functions, however, β can be computed semi-analytically for a given set of matrix properties (i.e., relative permeability, capillary pressure, permeability, and shape of a matrix block) and fluid properties (i.e., fluid viscosity) and therefore varies from grid block to grid block as well. Note that since we neglect gravitational forces in Equation 9, we focus only on spontaneous imbibition, that is, capillary-driven exchange between fractures and matrix. In this case, β is given by Geiger (2012, 2013) and Schmid et al. (2016)   where ℓ is a characteristic length, φ m ℓ can be understood as normalized pore volume in 1D, and A quantifies the amount of water imbibing into the rock matrix until a given time t. For a detailed derivation of β see SPOONER ET AL.
10.1029/2020WR027775 5 of 21 Schmid et al. (2011), Geiger (2012, 2013), Schmid et al. (2016). Using β, the recovery RF at time t due to spontaneous imbibition as a function of the ultimate recovery RF ∞ is given by Aronofsky et al. (1958) For completeness, we note that β can also be computed directly for a given set of matrix and fluid properties in cases where gravitational forces dominate (March et al., 2018) and Equation 9 includes the appropriate gravity term.

Retardation Factor for Spontaneous Imbibition in Dual-Porosity Systems
Next, we need to find a linear form of R that allows us to calculate τ ⋆ during capillary-driven fluid exchange between the fracture and rock matrix, using the locally uniform properties in each grid block. We seek a linear form of R to avoid any costly nonlinear iterations that, while possibly providing greater accuracy for the flow diagnostics, would also come at increased computational cost. Considering the analogy between Equations 4 and 9, we hypothesize that R for two-phase flow in a dual-porosity system must depend on β.
In a two-phase dual-porosity system, the retardation of the injected fluid is related to the strength of the transfer between the fractures and the matrix. In regions of the reservoir where the rate of transfer is fast, the saturation fronts are slowed down to a greater extent. The rate of transfer is affected by the surface area of the fractures, the porosity and permeability of the matrix, the saturation and viscosity of the fluids, and the wettability (e.g., Andersen, 2019;Andersen et al., 2014;March et al., 2016March et al., , 2018Rangel-German & Kovscek, 2002;Schmid et al., 2016;Schmid & Geiger, 2012. This dependency is captured by β. In order to obtain a linear form of R, we further hypothesize that the functional form of R depends on the Damköhler number, which is uniform in a grid block. We start the development of R for two-phase flow by using the concept of the dual-porosity Damköhler number . In each grid block,  quantifies the local balance between the advective fluid flow and fracture-matrix transfer (Andersen, 2019;Andersen et al., 2014;Rangel-German & Kovscek, 2002;Spooner et al., 2019) as The average of  within the drained regions of each producer, denoted by , quantifies the transfer-flow regime affecting individual wells. Note that it is convenient to express  and  in Log 10 because the values can vary over orders of magnitude depending on the wettability and permeability of the matrix or the density and connectivity of the fracture network. The following flow regimes exist: : Matrix-fracture transfer is slow compared to advective flow with possible risk of early breakthrough; the reservoir behaves closer to a single-porosity system where fluid flow is dominated by the fracture network.
: Matrix-fracture transfer is on a comparable timescale to advective flow where flow in fractures and fracture-matrix transfer are balanced. •  1 : Matrix-fracture transfer is fast compared to advective flow with limited risk of early breakthrough; the separation of time-scales in the mobile fractures and immobile rock matrix is no longer guaranteed and the fundamental assumption inherent to the dual-porosity concept may no longer hold.
Using the concept of , a simple model for the corresponding retardation factor in a given grid block can be defined as Equation 14 states that if β is large,  is large which results in an increased R, that is, travel times are longer because τ is increased. Conversely, when β is small,  is small, and R tends to 1, meaning that τ is influenced only by the fracture network and the reservoir behaves like a single-porosity system. Note that we set the lower bound of R  to 1 so that τ ⋆ = τ but do not assume an upper bound for R  , which can become much larger than 1 if  1  .
SPOONER ET AL.
10.1029/2020WR027775 6 of 21 As discussed above, adsorption is often modeled using isotherms. The Aranofsky model (Aronofsky et al., 1958) in Equation 12 provides such an isotherm and a corresponding retardation factor can be defined as where t ⋆ is the residence time of the fluid within a grid block. By definition the maximum value for R Aronofsky is 2. It is intuitive to re-scale R Aronofsky by accounting for the total volume of fluid that can be recovered from the rock matrix by considering an appropriate pore volume multiplier where PV is the volume of the mobile phase in the matrix PV = (1 − S nwr − S wr )φ m .
Note that once the velocity field has been obtained for the reservoir by solving Equation 1, R can be computed for each grid block considering the locally uniform values for β. u t , φ, , and t ⋆ . Therefore, Equation 7 remains linear subject to the aforementioned assumptions.

Dynamic Metrics for Dual-Porosity Systems
The final step of our methodological development is to use the concept of τ ⋆ and R developed above to expand some well-known metrics that approximate the flow behavior in a reservoir, for example, assess the influence regions around each well, estimate the risk of early breakthrough, or identify regions with poor displacement efficiency, that is, poorly swept regions in the reservoir. A range of dynamic metrics have been introduced in the past which account for more than just the geological heterogeneity, that is, the distribution of the static properties such as permeability and porosity; they quantify the geometry and heterogeneity of the flow behavior itself, respond to changes in the well placements and well controls, and can be readily used to rank individual models of a larger model ensemble (Krogstad et al., 2017;Møyner et al., 2014;Shahvali et al., 2012;Shook & Mitchell, 2009;Watson et al., 2020).
A classical way to quantify heterogeneity in a reservoir is the Lorenz coefficient, which compares the distribution of permeability and porosity along a well path to approximate the flow capacity (i.e., permeability times thickness k ⋅ h of a given layer) for a given storage capacity (i.e., porosity times thickness φ ⋅ h of the same layer) (Lake & Jensen, 1991). This measure, however, does not account for fluid flow between geological layers, aerial heterogeneity within a geological layer, or the effect of well placement and well controls on fluid flow. Shook and Mitchell (2009) defined a dynamic Lorenz coefficient L c to overcome this limitation. L c replaces φ ⋅ h with a swept pore volumes Φ and k ⋅ h with a flow capacity F that is based on the actual volumetric flow which is a function of both, geological heterogeneity and well controls, respectively well location. For a single-porosity system, F and Φ can be computed functions of τ by ordering the local pore volume of each grid block located at position  x based on its local value of ˆ. The functional forms for F and Φ are hence given by Møyner et al. (2014) and L c is defined as L c has been shown to be an excellent proxy for the displacement of one fluid phase by another during immiscible and incompressible two-phase flow (Krogstad et al., 2017;Møyner et al., 2014;Shook & Mitchell, 2009;Watson et al., 2020). It ranges from 0 to 1, with L c = 0 indicating perfectly homogeneous displacement during two-phase flow and L c = 1 indicating indefinitely heterogeneous displacement.
SPOONER ET AL.
10.1029/2020WR027775 7 of 21 Since the definitions for F and Φ depend on ordering the pore volumes of each grid block based on the local value of ˆ, a natural extension for a dual-porosity system is to order F and Φ based on ˆ. Since τ ⋆ is more heterogeneous than τ due to the additional fracture-matrix interactions, L c will also differ. We therefore define the dynamic Lorenz coefficient for dual-porosity systems c L  to indicate this additional heterogeneity.
Another important quantity to characterize immiscible two-phase flow displacement is the displacement (or sweep) efficiency E v , which compares the pore volume contacted by the injected fluid at time t to the total reservoir volume. For a single-porosity system, E v can be estimated as (Shook & Mitchell, 2009).
While E v can approximate the displacement in the fractures, it does not account for the fracture-matrix interactions. We therefore propose a sweep efficiency v E  that accounts for the fracture-matrix transfer processes in dual-porosity systems. The fracture pore volume that has been contacted by the injected fluid Φ f at a given time t can be computed using Equation 17 for any fractures where τ ⋆ ≤ t. The fracture-matrix transfer that occurs in the fractures which have been contacted by the injected fluid is an inherently time-dependent process which has been linearized in our approach, as discussed above. Fractures where τ ⋆ ≪ t have experienced fracture-matrix transfer for a much longer time than the fractures where τ ⋆ ≈ t. For any fractures where τ ⋆ > t, fracture-matrix transfer has not started yet because these fractures have not yet been in contact with the injected fluid. The corresponding matrix pore volume that was contacted due to fracture-matrix interactions Φ m can hence be approximated by Using Φ f , Φ m , and the total reservoir pore volume Φ T , which is independent of time, we define the dual-porosity sweep efficiency as

Application Examples
In this section, we present three application examples of increasing complexity to test the applicability and limitations of the proposed flow diagnostics. The new dual-porosity flow diagnostics were implemented and the corresponding simulations for all application examples were performed using the Matlab Reservoir Simulation Toolbox (MRST) (K. A. Lie, 2019; K. A. Lie et al., 2012). Unless noted otherwise, all full-physics simulations that provide the reference solutions were computed using the unified fracture simulation framework in MRST (March et al., 2021).

1D Analysis of Dual-Porosity Flow Diagnostics
We consider a simple 1D dual-porosity model (Figure 2) to test if the proposed dual-porosity flow diagnostics can approximate the time it takes a wetting phase to travel from the injector to the producer. We compare the breakthrough times predicted from flow diagnostics using the three different retardation factors (Equations 14-16) to full-physics simulations that compute the advance of the wetting phase and corresponding fracture-matrix fluid exchange fully transiently. The full-physics simulations ignore compressibility, heterogeneity, gravity and capillary forces, and use the transfer function given by Equation 10 Table 1 for further parameters.
The 1D model is homogeneous and its parameters are detailed in Table 1. Fracture relative permeabilities are linear, that is, k r = S. Fluid viscosities are equal and constant, with μ = 1cP. We consider five different cases of fracture-matrix transfer rates. They range from Case 1 with σβ = 6 × 10 −9 m −2 s −1 to Case 5 with σβ = 6 × 10 −5 m −2 s −1 . σ is a shape factor, quantifying the fracture-matrix surface area and approximating the shape of the matrix block. For Case 1, the transfer rate is typical for mixed-wet rock under spontaneous imbibition while for Case 5, the transfer rate is typical for strongly water-wet rock. Note that these values were fixed rather than computing β from Equation 11 for a given set of fluid and rock properties. Figure 3 shows the evolution of the saturation profiles in the fractures and matrix. These results were obtained from full-physics simulations and help us to analyze what processes can causes differences between the flow diagnostics approximations, which are calculated using a linearized steady-state solution, and the full-physics simulations that capture the nonlinear and time-dependent interplay of fracture-matrix transfer. For Case 1 (slow transfer), the injected wetting phase progresses quickly through the fractures. The saturation gradient between the fractures and matrix behind the saturation front is high since little transfer has occurred. This means that the matrix behind the saturation front is still capable of exchanging fluid with the fractures. In this region, the injected fluids are hence still subject to delayed movement due to fracture-matrix transfer. Similar behavior can be observed for Case 2. For both cases the saturation gradient remains high even after the fraction of the produced wetting phase reaches 50% at the production well. When the fracture-matrix transfer is slow, transfer can occur behind the saturation front as long as S m < S m⋆ (Equation 10). For Case 3, the matrix saturation reaches values similar to the fracture saturation early on, that is, at the time of breakthrough, indicating fast equilibrium between the fracture and matrix and fracture-matrix transfer behind the saturation front stops quickly since S m ≈ S m⋆ . Although not shown here, similar behavior is observed for the cases with fastest fracture-matrix transfer (Cases 4 and 5). The change in transfer behavior is also apparent in the  values for each case, which increases from −1.95 (Case 1) to 2.05 (Case 5). For Case 3,  0.05  , indicating a near perfect balance between fracture flow and fracture-matrix transfer. Figure 4 compares the different breakthrough times, calculated after the fraction of the injected wetting phase has reached 0.5% at the producer. We compare the different breakthrough times rather than the saturation evolution in fracture and matrix because the flow diagnostics simulations are steady-state and linear in nature, that is, they do not preserve the history of fracture-matrix fluid exchange. It is therefore not surprising that we observe differences between the estimated breakthrough times from the dual-porosity flow diagnostics and those estimate from full-physics simulations. Flow diagnostics using R Aronofsky underpredict the simulated breakthrough time in all cases. The bounded form of R Aronofsky limits the transfer for the cases where fracture-matrix transfer is fast, which is required to capture the s-shape behavior in the breakthrough time ( Figure 4). However, R Aronofsky does not capture the magnitude in the breakthrough time accurately. Flow diagnostics that use R AronofskyPV account for additional storage and increase R, but overpredict the simulated breakthrough for all but the lower transfer rate cases (Cases 1 to 3) and vastly overpredict the breakthrough times for Cases 4 and 5 because R AronofskyPV can reach much larger values than R Aronofsky if β is large. Since R  can also reach very large values for cases where β is large, the breakthrough times are also overpredicted for Cases 4 and 5.
Recall that in the previous section we assumed piston-like displacement in the fractures because the relative permeability is assumed to be linear and μ w ≈ μ nw . If there is a strong contrast in viscosities, or the relative permeability is nonlinear, saturation gradients will occur in the fracture, which will further complicate the temporal evolution of fracture-matrix transfer and likely lead to further differences between the simplified flow diagnostics approximations and full-physics simulations.
This 1D analysis shows how variations in fracture-matrix transfer can lead to situations where R changes transiently and a linearized version of R is an over-simplification. It is possible to capture such nonlinear SPOONER ET AL.   (Bourbiaux, 2010). Instead, the matrix becomes the more mobile domain and flow is dominated by the nonlinear diffusion processes occurring during spontaneous imbibition, requiring single-porosity flow diagnostics with appropriately averaged fracture and matrix permeabilities (e.g., Wong et al. 2020) or a dual-permeability approach.
We also note that these complex nonlinear and time-dependent interactions between fracture flow and fracture-matrix transfer render the full-physics simulations more sensitive to grid resolutions while the SPOONER ET AL.

10.1029/2020WR027775
10 of 21  (Table 1) before breakthrough and at a later time. The bar charts show the normalized saturation difference between the matrix and fractures for each case. Note that the behavior for Cases 4 and 5 are very similar to that of Case 3 and are hence not shown.
steady-state nature of the flow-diagnostics simulations exhibit hardly any sensitivity to the grid resolution (Spooner, 2021).

Dual-Porosity Flow Diagnostics in 2D With Geological Heterogeneity
The previous example shows that using R to correct τ for matrix-fracture transfer can approximate the breakthrough time in a dual-porosity system as long as  1  . We now consider a more complex 2D example where the matrix properties vary spatially. We previously have shown that flow behaviors, swept and drained reservoir regions, and individual well performance can be assessed using  (Spooner et al., 2019). We now extend this work to evaluate how τ ⋆ can be used to estimate breakthrough times.
The parameters for the 2D dual-porosity model are detailed in Table 2. The model has four different rock types, based on different relative permeability curves ( Figure 5) which have been distributed across the model ( Figure 6). β for each rock type was computed using Equation 11. To accentuate the differences between flow rates in the fracture and matrix, the resulting values for β have also been multiplied by a factor of 0.1 and 10, which is identical to decreasing or increasing the matrix permeability by a factor 10, respectively. The fracture porosity and permeability are uniform, fracture relative permeabilities are linear, fluid viscosities are equal, and the wells are spaced evenly and operate at the same drawdown pressure and injection pressures. Hence the flow field in the fractures is symmetric and the only heterogeneity arises from the different rates at which fractures and matrix exchange fluids, leading to spatial variations in  and significant differences in τ and τ ⋆ (Figure 6). Figure 7 displays the comparison of the breakthrough times for each well using full-physics simulations and the breakthrough times approximated using flow diagnostics with R  , R Aronofksy , and R AronofksyPV . The results show that the best and poorest performing wells in terms of breakthrough time are always correctly identified, independent of which formulation of R is used. However, as in the 1D model, breakthrough times predicted from full-physics simulations and flow diagnostics diverge if  1  and flow in the fractures is slower than fracture-matrix transfer, that is, when the fundamental assumption underpinning the dual-porosity model no longer holds (Bourbiaux, 2010). This behavior is particularly apparent for Wells 1 and 4, which are located in the regions where the rock matrix is water-wet and hence fracture-matrix transfer is highest ( Figure 6). In particular, R AronofksyPV overpredicts the breakthrough times in these situations because, like in the 1D model, it estimates the highest R values.
As discussed for the 1D model, in cases where β is large and  1  , the rates of the fracture-matrix transfer continue to change even once the injected fluid has broken through; sometimes the transfer rates are strong enough to completely impede the advancement of the saturation front until the matrix blocks are fully depleted. The linearized formulation of R and the steady-state nature of the flow diagnostics calculations lack the memory of the saturation history, meaning there is no memory of the spatio-temporal changes in saturation in the rock matrix (Tecklenburg et al., 2016). We reiterate, however, that the dual-porosity flow diagnostics accurately predict the order of the wells at which they experience breakthrough and as long as the dual-porosity assumption holds, that is, for situations where  1  , breakthrough times can be approximated with reasonable accuracy.
To summarize, expanding τ to τ ⋆ using R provides reasonable quantitative estimates of breakthrough times within the following limits: SPOONER ET AL.      . Note that τ in (d) is nonuniform and exhibits the symmetric shape typical to a five-spot well pattern in a homogeneous porous media but appears to have the same value due to the scaling used for τ ⋆ (e). For further model parameters see Table 2.
• Fast transfer: When  1  , we cannot assume a constant retardation factor as the approach overpredicts breakthrough times. However, in this situation the fractures no longer dominate the flow behavior so the dual-porosity model is not appropriate and dual-permeability or single-porosity models with average fracture and matrix properties should be applied. • Intermediate transfer: When    1 1  , fracture-matrix transfer is still possible behind the saturation front but we can assume a constant retardation factor. R Aronofsky predicts breakthrough times more accurately compared to R  because R Aronofsky does not allow for transfer once the matrix has been fully depleted. • Slow transfer: When  1  , R  tends to provide slightly more accurate breakthrough time predictions compared to R Aronofsky and in extreme cases when τ ⋆ ≈ τ, a single-porosity will approximate breakthrough times with sufficient accuracy.
Furthermore, the extended diagnostics can offer qualitative, visual information to inform us of regions where the injected fluid is likely to travel slower or faster, and allow us to assess which wells are most and least likely to experience early breakthrough (Spooner et al., 2019). Such insights are not only useful to understand the limitations of dual-porosity flow diagnostics; they can also help to improve the design of full-physics reservoir simulation model, indicating areas in the reservoir where flow and transport can be modeled by a dual-porosity system and areas where other model concepts (e.g., dual-permeability or single-porosity concepts) are more appropriate.

3D Field Application: Teapot Dome Reservoir
We now present a case study which applies the dual-porosity flow diagnostics to a real fractured reservoir, the Teapot Dome reservoir in Wyoming, USA. We previously analyzed how single-porosity diagnostics, , and  can be used to characterize flow behaviors and compared this analysis to full-physics simulation results (Spooner et al., 2019). We now extend this work by including an analysis of breakthrough times using the dual-porosity flow diagnostics and investigate how flow diagnostics can be deployed to identify possible well controls that allow us to manage the reservoir more effectively. In contrast to the 1D and 2D models, the full-physics simulations are now carried out using a commercial simulator (IMEX by CMG Ltd) with the default Gilman and Kazemi dual-porosity model option (Kazemi & Gilman, 1993

Overview of the Teapot Dome Reservoir
The Teapot Dome reservoir is an elongated North-South trending anticline located at the South-West edge of the Powder River Basin in Wyoming. The Pennsylvanian Tensleep formation forms the main reservoir. The formation consists of aeolian sands, interbedded with sabkha and marine dolomites (Ouenes et al., 2010). The Teapot Dome is considered to be predominantly a Type II fractured reservoir according to Nelson's classification (Nelson, 2001), meaning that the reservoir is a dual-porosity system where fractures provide the essential permeability and the matrix only contributes to flow. The Discrete Fracture Network (DFN) and simulations models used in this study were built by Ahmed Elfeel et al. (2013) using input data described by Cooper et al. (2006) and Schwartz (2006). The final simulation model contains 267,294 active grid blocks. There are 10 producers and three water injectors, all of which are bottom hole pressure (BHP) controlled. The initial average pressure of the reservoir is 2,290 psi, the pressure is set to 2,300 psi at each injector and to 2,100 psi at each producer. A wetting phase is injected to displace the nonwetting resident fluid. A full-physics simulation takes around 2 hours to simulate 3,500 days of production on a standard workstation with four core CPU with 32 GB RAM and 4 GHz processor speed.
The matrix has average permeability of 17 mD, an average porosity of 9%, and is water-wet. There is only one rock type, corresponding to the most water-wet relative permeability and capillary pressure curve in Figure 5. Fracture relative permeabilities are assumed to be linear and fluid viscosities are equal. The fracture permeability varies over nearly three orders of magnitude, leading to a wide range of τ ⋆ values which are generally lower in the vicinity of the wells (Figure 8). Early water breakthrough is common at the Teapot Dome Reservoir, as evident in the distribution of  and  (Figure 8). However, areas that are poorly swept exist away from the wells, which is reflected by large τ ⋆ values and regions where  is positive. The variance in  is low, which indicates flow behavior is influenced by the wells and not by the differences in matrix-fracture transfer across the reservoir. The relatively high values of  and  can be expected for a water-wet reservoir. Given the above discussions, flow diagnostics will likely overpredict breakthrough times.
In addition to the possible discrepancies between the breakthrough times predicted using dual-porosity flow diagnostics and full-physics simulations caused by the relatively high  values, flow diagnostics also do not account for gravitational effects. The location and shape of saturation fronts after 10 years of water injection modeled with flow diagnostics and full-physics simulations show relatively good agreement, although differences exist (Figure 8). For example, the shape of the saturation fronts on either side of the crest of the reservoir is clearly affected by gravitational effects in the full-physics simulations. Figure 9 compares the breakthrough times from full-physics simulations with the breakthrough times approximated using single-and dual-porosity flow diagnostics. We include the single-porosity flow diagnostics results here to show how the additional fracture-matrix transfer processes, which are not captured by traditional flow diagnostics, influence our assessment of possible well behaviors. Breakthrough times computed using R Aronofsky are close to the single-porosity breakthrough times, and slightly underpredict the breakthrough times for some wells (10_46, 12_76, 4_73) but provide good approximations in general. The 1D simulations have shown that R  and R AronofskyPV can vastly overpredict breakthrough times in cases where  1  , which can also be observed for certain wells (e.g., 2_76 and 6_72). Well 6_72 is also located in the more steeply dipping flank of the reservoir where gravitational forces are larger, which likely exacerbates the differences between flow diagnostics and full-physics simulations. For larger inter-well spacing, the lack of time dependency of the retardation factor discussed in the 1D example above leads to an increased error in the predicted breakthrough times as the overestimation occurs in a much larger area of the reservoir, that is, in every grid cell located between injector and producer. This is particularly apparent at well 12_76 where the breakthrough time from dual-porosity flow diagnostics based on R AronofskyPV exceeds 4,000 days and based on R  exceeds 30,000 days, compared to a breakthrough time of 1,318 days from full-physics simulations. Figure 10 shows the ranking of the individual wells in terms of their breakthrough times. Despite differences in predicting the exact breakthrough times, R Aronofsky , and R AronofskyPV correctly identify the best and worst performing wells in terms of breakthrough times, that is, the wells that are most likely to experience SPOONER ET AL.

10.1029/2020WR027775
14 of 21 breakthrough very early or very late. Where the ranking of the wells differs, it is often due to relatively small (i.e., order of 10s of days) differences in breakthrough times which, from an operational point of view, would likely be negligible. While the ranking of wells based on R Aronofsky and R AronofskyPV is similar, R  provides a very different ranking in some cases. In particular, the worst performing well is incorrectly identified with R  . However, in cases where R Aronofsky and R AronofskyPV do not correctly identify the ranking, R  does (e.g., wells 4_73 and 13_67-1).
In general, the single-porosity flow diagnostics underpredict and the dual-porosity flow diagnostics overpredict breakthrough times. The range between the single-and dual-porosity breakthrough times can therefore be used to estimate the lower and upper bounds at which breakthrough is expected to occur. Figure 11 SPOONER ET AL.  shows the breakthrough times for single-and dual-porosity flow diagnostics using R AronofskyPV . For 90% of the wells, the real breakthrough times lie between the lower and upper bounds of the single-and dual-porosity breakthrough times with only the simulated response for well 6_72 falling just outside the range. This result is unsurprising given the well's location down-dip and its distance from an injector, that is, the fact that gravitational forces dominate in the vicinity of this well.
Although our dual-porosity flow diagnostics do not account for the time-dependent and dynamic changes in fracture-matrix transfer, comparison with the full-physics simulations result shows that the approach provides useful approximations to the breakthrough time for the Teapot Dome reservoir, especially in the regions of the reservoir where  1  and the concept of the dual-porosity model can be applied. The breakthrough times estimated using single-porosity flow diagnostics provide a useful indication of the lower range of breakthrough times that can be expected and hence provide complementary information to the dual-porosity flow diagnostics which tend to overestimate breakthrough times. Figure 12 shows the F-Φ curves based on τ and τ ⋆ for each retardation factor as well Ev ⋆ . F-Φ curves for τ and τ ⋆ based on R Aronofsky and R AronofskyPV are similar and hence L c and c L  are close. In contrast, the F-Φ curve based on R  and the resulting c L  coefficient indicate increased heterogeneity. Recall that the breakthrough times for R  for wells 6_72 and 12_76 are significantly overestimated (Figure 9), which is caused by the difference in τ ⋆ across the reservoir that manifests itself as an increased dynamic heterogeneity indicated by the F-Φ curve and larger c L  coefficient.

Flow Diagnostic Metrics for the Teapot Dome Reservoir
Ev ⋆ curves differ depending if τ ⋆ is calculated based on R  , R Aronofsky , or R AronofskyPV (Figure 12). τ ⋆ calculated using R  is larger overall, as evident in the overpredicted breakthrough times (Figure 9). Ev ⋆ based on R  is therefore lower than Ev ⋆ based on R Aronofsky and R AronofskyPV . Ev ⋆ based on R Aronofsky is highest because R Aronofsky is by definition limited to a maximum value of 2 (Equation 15), which yields lower τ ⋆ values across the reservoir. Lower τ ⋆ values mean that the injected fluids travel faster through the reservoir and contact a greater volume for any given time.
There is an important difference between Ev ⋆ and c L  : Ev ⋆ quantifies the speed at which the reservoir volume is contacted while c L  quantifies the heterogeneity in the shape of the connected flow paths as determined by τ ⋆ . Since R Aronofsky and R AronofskyPV predict different magnitudes of τ ⋆ but do not change the relative distribution, the corresponding c L  values are similar. However, Ev ⋆ curves differ because of the different time scales, given SPOONER ET AL.

10.1029/2020WR027775
16 of 21  by τ ⋆ , at which fluids are predicted to travel through the reservoir. Dual-porosity flow diagnostic metrics such as Ev ⋆ and c L  should therefore be used in combination to assess the dynamic heterogeneity in fractured geological formations.

Optimizing Displacement Efficiency in the Teapot Dome Reservoir
We now investigate how we can use Ev ⋆ and c L  to guide the optimization of the well controls to increase the efficiency at which the injected fluid displaces the resident fluid, aiming to increase the cumulative volume of resident fluid produced. Such optimization is possible because Ev ⋆ and c L  are a function of τ ⋆ , which depends on both, the reservoir geology and the well controls. The main challenge with all optimization algorithms is the need to run a very large number of simulations to explore the parameter space adequately, which is computationally expensive. Technologies such as proxy models (e.g., Josset et al., 2015;Oladyshkin et al., 2011;Schulte et al., 2020) or streamline simulations (e.g., Datta-Gupta & King, 2007) can be deployed to reduce the CPU time needed for optimization studies. Here, we focus on using flow diagnostics and the resulting metrics to guide optimization studies and reduce computational overheads, following the work of Møyner et al. (2014).
We consider seven optimization exercises which examine the effectiveness of flow diagnostic metrics for optimizing the well controls to enhance the displacement efficiency. The constraints for these optimization studies aim to change the pressures at the wells without restricting the fluid injection rates and are therefore expected to lead to the well-known results that the more fluids are injected, the more fluids can be produced. We also do not consider an ensemble of reservoir models to optimize the reservoir dynamics while considering geological uncertainties. Note that this optimization exercise is only a proof-of-concept to demonstrate that the flow diagnostics framework can be used as a "pre-processor" for optimization studies but real field applications would require more complex optimization constraints in terms of well controls and possibly the consideration of uncertainties inherent to the geological model. The objective for each of the seven optimization tests is described as follows:  Minimizing c L  aims to reduce the dynamic heterogeneity. Maximizing of Ev ⋆ computed at 1 year of continuous injection aims to increase the total contacted volume of the reservoir and therefore increase the production of the resident fluid. To optimize flow in the reservoir, we applied a particle swarm optimization (PSO) algorithm to minimize c L  or maximize Ev ⋆ , respectively. PSO is one example of stochastic optimization algorithms (Kennedy & Eberhart, 1995;Mohamed et al., 2010). The PSO algorithm samples the distributions of BHP controls for the injectors and producers. The pressures are allowed to vary from 2,400 to 5,000 psi for the injectors and from 50 psi to 2,100 psi for the producers, using uniform prior distributions. PSO was set to run for a maximum of 2,000 iterations for each optimization exercise. Figure 13 shows the optimized BHP values for each well as determined in each optimization exercise. The optimal well controls that emerge from optimizing L c and c L  based on R Aronofsky and R AronofskyPV are similar for 9 out of 13 of the wells, indicating an increase or decrease in the same wells. The optima that emerge from optimizing c L  based on R  differ. This is to be expected, considering the differences in the F-Φ curves discussed earlier (Figure 12). The optimal well controls that emerge from optimizing Ev ⋆ exhibit a very similar behavior, causing an increase in injection pressures and decrease in producer pressures at the same wells. Figure 14 shows the results from full-physics simulations that were carried out using the optimal well controls as identified by minimizing c L  or maximizing Ev ⋆ , respectively. As expected, simulations with improved well controls, which lead to higher injection rates, cause the volumes of the produced resident fluid to increase. Optimal well controls that are based on dual-porosity metrics (Exercises 2-7) show better displacement efficiency, that is, improved production forecasts, compared to well controls based on the single-porosity metric (Exercise 1). Optimal well controls based on minimizing c L  (Exercises 2-4) yield higher production forecasts than the optimal well controls based on minimizing L c (Exercise 1), with the production forecasts based on R AronofskyPV being close to but still higher than those based on L c . The highest production forecasts are obtained by optimizing the well-controls based on maximizing Ev ⋆ , and forecasts are largely independent of the way how R is being approximated. The differences in optimized production forecasts arising from minimizing c L  or maximizing Ev ⋆ , respectively, reflect the differences in how these two metrics respond to the variations in τ ⋆ across the reservoir. By optimizing metrics based on τ ⋆ , not just τ alone, we can improve the well controls to better account for the fracture-matrix processes that influence fluid flow in dual-porosity systems.

Summary for the Teapot Dome Reservoir
The Teapot Dome reservoir is a water-wet reservoir with comparatively fast matrix-fracture transfer rates, and hence significant parts of the reservoir have  1  . This reservoir is therefore at the limits of the applicability of the dual-porosity flow diagnostics, meaning that breakthrough SPOONER ET AL.  times are overpredicted at individual wells. Nevertheless, the dual-porosity flow diagnostics still provide useful insights into the reservoir dynamics and hence help to improve our understanding of how geological heterogeneity impacts fluid flow.
Most notably, the overall ranking of the wells in terms of breakthrough time is captured correctly. Furthermore, τ and τ ⋆ provide reliable lower and upper boundaries for estimating breakthrough curves. The corresponding metrics, c L  and Ev ⋆ , can be used to guide the identification of optimized well controls that improve displacement efficiency and fluid production in the reservoir. Similar to the discussions of Watson et al. (2020), we also observe that care must be taken when selecting flow diagnostic metrics to assess the reservoir dynamics because different metrics emphasize different aspects of the complex flow behaviors occurring in fractured formations.
For this realistic 3D model, flow diagnostics have been computationally very efficient. Each optimization run with 2,000 iterations took ∼7 hours using four cores on a standard workstation type PC (32 GB RAM and 4 GHz processor speed). This is equivalent to 50 seconds per flow diagnostics run, which includes setting up the model and exporting results after each run. Recall that in contrast a single full-physics simulation run takes around 2 hours on the same workstation, rendering an optimization with 2,000 iterations computationally very challenging if not intractable. Dual-porosity flow diagnostics therefore provide a quick method to investigate the approximate flow behaviors in fractured geological formations and estimate the dynamic response for a reservoir model, allowing the many permutations of the design parameters (e.g., well controls) that need to be considered in a detailed optimization study that employs full-physics simulations.

Conclusions
Flow diagnostics can approximate the reservoir dynamics in a matter of seconds, providing visual and quantitative information than can be used to interpret reservoir flow behavior and quantify dynamic heterogeneity. Flow diagnostics can be applied to modeling flow processes in a variety of subsurface systems (e.g., CO 2 storage, geothermal energy, oil recovery, and hydrology) where geological uncertainty is prevalent and fast screening methods are required to consider the range of scenarios required to capture uncertainty.
In this study, we extended traditional flow diagnostics to account for fracture-matrix interactions in dual-porosity or mobile-immobile systems. This extension is based on a linearized retardation factor R, which approximates the time-of-flight τ ⋆ in a fractured geological formation while estimating the delay an injected wetting phase experiences when it interacts with the rock matrix due to capillary forces and spontaneous imbibition. We proposed extended metrics based on τ ⋆ , namely the dynamic Lorenz coefficient c L  and the dual-porosity sweep efficiency Ev ⋆ , which are proxies for the dynamic heterogeneity in the reservoir and the efficiency at which the injected fluid contacts and displaces the resident fluid as a function of time, respectively. Care must be taken, however, when selecting metrics because different metrics emphasize different aspects of the reservoir dynamics.
We proposed and tested three formulations for R. No single formulation provided robust approximations to the breakthrough time in all cases because R is assumed to be linear and constant in time. Fracture-matrix interactions becomes complex and time-dependent in cases where flow in the fractures is equal to, or slower than, fracture-matrix transfer, that is, when the average Damköhler number  1  . In these situations, a dual-porosity model, which assumes a separation of time-scales between fracture and matrix, is no longer applicable. However, in many practical situations τ ⋆ , c L  , and Ev ⋆ provide reasonable approximations of the reservoir performance, including well-by-well breakthrough times. The single-porosity time-of-flight τ and τ ⋆ furthermore provide convenient estimates of the range of breakthrough times that can be expected for each well.
The applicability of these metrics was further demonstrated in a proof-of-concept optimization study where improved well controls were identified for a real fractured reservoir based on minimizing c L  and maximizing Ev ⋆ , respectively. Each flow diagnostics calculation took less than 1 minute compared to over 2 hours run time for a full-physics simulation. This significant reduction of CPU time therefore enables more efficient optimization studies where dual-porosity flow diagnostics can be used to explore large design spaces to identify possible well controls, which are then subsequently refined using additional full-physics SPOONER ET AL.
simulations. Flow diagnostics are therefore not intended as a replacement to full-physics simulation, but instead offer a complementary technology to accelerate optimization and uncertainty studies in fractured geological formations.

Data Availability Statement
All flow diagnostics and 1D and 2D simulations were carried out using SINTEF's open source Matlab Reservoir Simulation Toolbox MRST (https://www.sintef.no/projectweb/mrst/). Earlier versions of this study have greatly benefited from the careful and constructive comments by the reviewers as well as the examiners of Victoria's PhD thesis.