Optimization of waterflooding performance by using finite volume-based flow diagnostics simulation

From a visual point of view, volumetric information about reservoir portioning and communication such as sweep, flow patterns, and drainage zones are longer better interpreted and pictured when presented by an average volumetric flux calculation. To this hand, finite volume discretization can be used to substitute streamline simulation-based finite difference to assess flow diagnostic information. Herein, we use finite volume-based flow diagnostics to optimize waterflooding. In particular, we discretize in finite volume the flow equation from Darcy’s law single-phase incompressible flow steady state combining with two auxiliary flow equations, time of flight and stationary tracers using the two-point flux approximation to describe fluid particles motion and flow lines. In addition, with the estimation of dynamic heterogeneity, we compute the Lorenz coefficient to highlight the reservoir flow and storage capacity characterization. To optimize waterflooding rates, we first, use an objective function the equalized Lorenz coefficient got through the evaluation of average travel time in cells to increase sweep efficiency and decrease the dynamic heterogeneity coefficient. Second, following the same target, we use the flow diagnostic interactive tools to study the volumetric sweep displacement front and harmonize the flooding breakthrough. In this work, our conceptual approach is to see the reservoir initially filled with oil; then, optimizing the Lorenz coefficient leads us to an oil recovery improvement. To be pragmatic, we apply our waterflooding performance optimization model on two case studies, the ninth SPE comparative solution project, a reexamination of black-oil (synthetic case) and ZHNBA Chinese oilfield (real field dataset).


Introduction
The road toward strategies for recovery improvement of hydrocarbon resources compulsory goes through reservoir rock characterization and the study of fluids flow in porous media. Herein, reservoir modeling is taken as an investigation of the interactive motion of fluid particles, depicting the physics and mathematical equations involved. This means disassembling the partial differential equations describing the fluid behaviors in the reservoir and then looking for numerical solutions and implementation. Thereby, the flow diagnostics workflow module from MATLAB Reservoir Simulation Toolbox (MRST) is the perfect modeling tool that can be used as a waterflooding project processor to provide a reservoir model to establish connections and quickly procure a qualitative picture of the flow patterns. Flow diagnostics are the result of flow field analysis properties, and they are the study of well responses and fluid displacement front distribution to understand the flow paths and communication patterns in reservoir models. Howsoever, flow diagnostics-based finite volume discretization can be used to measure and optimize the dynamic heterogeneity, and increase well volumetric sweep efficiency. Flow diagnostic with the interactive tools can be used to manage the flooding breakthrough, and forecast an optimal injection well rates configuration.
The study and use of flow diagnostic tools have been the purpose of several types of research. Authors have mostly involved flow diagnostics to develop strategies for hydrocarbon resources optimal recovery. However, we have distinguished two main approaches, the streamline-based finite difference, and finite volume-based flow diagnostics.
Streamline-based flow diagnostics has received attention these last years; it is a simulation method based on the finite-difference discretization. This started with modeling reservoir geostatistical realizations based on the study of time flight equations by Idrobo et al. (2000), and then Thiele and Batycky (2003) with an automatic method for waterflooding and well rates optimization. Then, later the same year, Ates et al. (2003) ranked and up-scaled geostatistical reservoir models, and then Akhil Datta- Gupta and King (2007) stated the features of streamline technology and complement to conventional finite-difference simulation. Batycky et al. (2005) revisited class flood surveillance methods applied to injection production data, and Park and Datta-Gupta (2011) worked on the workflow for waterflood rate optimization using streamline-based flood sufficiency maps. Lastly, Haegland et al. (2019) simulated fluid flow and transport in porous media.
Streamline-finite difference-based method has been efficient for these above purposes, but it has limitations in terms of computational complexity in its extensibility to irregular computational domains like complex geological structures such as faults. However, reservoir drainage and volumetric sweep information can be well-interpreted and pictured when presented by an average volumetric flux calculation.
Recently, finite volume discretization has been explored as a mathematical ground for reservoir simulation methods. In the waterflooding framework, we can quote Shook and Mitchell (2009) for their work on the new method for estimating heterogeneity in the earth model using time of flight and volumetric flow rate information. Then, Shahvali et al. (2011) proposed finite volume method as an alternative to streamline for obtaining flow diagnostic information. Finally, Møyner et al. (2015) proxied flow diagnostics for reservoir optimal management workflows.
This work is purposed to decrease the impact of reservoir rock dynamic heterogeneity known as one of the unfavorable factors affecting waterflooding overall sweep efficiency. The dynamic heterogeneity can be represented as the imbalance permeability distribution in the reservoir leading to an early breakthrough in high permeability zones and stagnant region in low permeability zones. Our work will be to find ways to equalize the measure of dynamic heterogeneity by computing its values using grid cells volumetric average approximation. The resulting target at the end of this work will be: • Having dynamic heterogeneity volumetric grid cells average values (optimal values) • Harmonizing the flooding breakthrough • Optimizing the water injection well rates • Optimizing the waterflooding performance

Mathematical model
Here, we simply present the governing equations sustaining this work. We first start by developing the single-phase incompressible steady-state flow equation from Darcy's law.
where p is the pressure of the fluid, is the density of the fluid, g is the acceleration from gravity, Z the depth is a vector function of (x, y, z) which has the g direction. However, K is shown to be the absolute permeability tensor measured in Darcy's units the equivalence of length squared in the mathematics of reservoir simulation, Ewing (1983). In most reservoir simulation using Darcy's law, the permeability ( K ) is considered as a special diagonal tensor.
with K x , K y and K z as the permeability in the (x, y, z) directions. For K x = K y = K z , isotropic medium; otherwise, anisotropic. The mathematical model for the incompressible singlephase flow is derived by the fusion of the conservation of mass, Darcy's law on steady state and the state equation on a control volume ( v ) shown by Lie 2016. The mass conservation implies that the mass accumulated in v equals the mass flow across the boundary of v plus the mass injected into v via wells with , the fraction of volume, v available for flow, the fluid density per unit volume, v the computational domain with normal vector ⃗ n , ⃗ u the macroscopic Darcy's velocity and q the mass flow rate or fluid sources or sinks (inflow, outflow) per volume unit at certain locations.
The mass conservation equation is given as: By using the Gauss divergence theorem, we get the following transformation: This equation is valid for any volume ( v ) and also must satisfy the continuity equation To continue our process, we have to introduce constitutive equations which are the link between two different states (pressure, volume, temperature) of the system at specific physical given conditions. Herein, in our work, we can take an example, Darcy's law equation, which links the velocity ( ⃗ u ) to the fluid pressure ( p ). Herein, a variation in density will cause a pressure and temperature variation; in thermodynamics, those variations can be considered as a change in volume ( v ) for a fixed number of particles.
T and p are the indication that the above process takes place at constant pressure and temperature. Also, we can notice that for a fixed number of particles, v can only be constant, d v = dv ; then, the last equation can be rewritten as follows: Herein, c f is the isothermal compressibility and f is the thermal expansion coefficient. Many subsurfaces keep the density almost constant which makes the conduction to maintain the temperature constant that situation drives us to the following simplification.
c f , the isothermal compressibility, is a nonnegative factor and depends on pressure and temperature, c f (p, T). Reintroducing Darcy's law and rock and fluid compressibility in the continuity in Eq. (5), we get a parabolic equation described as follows: Herein, c t = c r + c f , the total compressibility, the parabolic equation is nonlinear because and c t depend on p . In the case of incompressible flow, we have to make the assumption that c t = 0 , and then we get an elliptic equation with variable coefficients which represent the single-phase incompressible flow equation.
When bringing the fluid potential = (p − g z) , the single-phase incompressible flow equation can become the generalized Poisson equation: −∇K∇ = q or when there is no source, neither sink can become Laplace equation: −∇K∇ = 0 , all for steady-state flow.
Secondly, we introduce the flow transport auxiliary equations: the time of flight and tracer equations. We designate, as time of flight, the time a particle takes to travel from a source or a sink to a particular point in the reservoir; it is the time a fluid particle takes to travel a distance (a) in the interstitial velocity ⃗ u . The time of flight is an auxiliary equation which describes the fluid flow lines in the reservoir, and it is an important equation for the visualization and analysis and can be written as follows: Here, τ is the time of flight; therefore, using directional derivative, we can write: From this last expression (12), we can write the time of flight equation as follows: We also introduce a set of tracer equations that characterized neutral particles that passively flow with the fluid without changing the flow properties. The reservoir tracer concentrations are expressed by the following continuity equation: Then, by simulating the evolution of non-diffusive tracers that keep constant values during the fluid expansion or compression, we get the communication between patterns. From this point, we set at certain inflow boundary part tracer concentration that equals one, and we compute the solution from the non-conservative equation steady-state conditions: The resulting tracer distribution from the last equation defines the reservoir portion which will be influenced by a coming inflow boundary flow. By reversing the flow field sign, we can identify the reservoir portion that will be influenced by an outflow boundary or sink. Going further by stretching this method to all parts of source flow, we can compute the instantaneous flow field.

TPFA of the single-phase incompressible flow equation steady state
We start the discretization by rewriting the steady-state equation by its integral form which is the simple form of Eq. (3), where and have been removed because they are dependent on time.
The next step of our discretization is to compute the flux which goes through each face of the cells. Here, we are talking about our two neighbor cells ( i , c ). We assume that both grid cells are supposed to match, and then we assign the face i,c as half face of the grid cell i that associates with the normal vector ⃗ n i,c , and we can write as follows: Since the grids are assumed matching, the half face i,c has his twin i,c and the area A i,c = A c.i , but ⃗ n i,c = −⃗ n i,c opposite normal vector. Now, using Darcy's law and the mid-point approximation, we can derive the flux as: where ⃗ x i,c is the centroid on i,c . Using the finite volume method makes us assume that the pressure is linear or constant inside each cell because finite volume discretization is the approximation of the average quantities inside cells; then, we get: We can see in this last equation the introduction of T i,c which is one-side transmissibility or half transmissibility that links the flux across a cell precisely between the cell center and the edge. Finally, let us apply the flux continuity through all faces as follows.
u i,c = −u c,i = u ic also the continuity of the pressure of edge faces i,c = c,i = ic . Then, we get: When we remove the interface ic , we get the two-point flux approximation (TPFA) Here, T ic is the transmissibility between both cells ( i , c ). TPFA is a method that uses the average pressure of two cells through an interface to approximate the flux across those cells. To conclude, we clearly see that TPFA satisfies the equation ∇ ⋅ ⃗ u = q as follows:

TPFA of time of flight and tracer equations
By combining the time of flight and tracer equations, we can write the flow transport equation (Lie 2016): . We can then start to discretize our steady-state flow transport equation by applying the Gauss divergence theorem over a grid cell i .
In the TPFA section, we described the discretization of an interface flux ic between i and c with v ≡ 1 . In this case, we can find an approximation value of v ic on ic to write the interface flux as follows: For the reason of scheme stability avoiding oscillations, we need to define the upwind value of v ic Then, we can end the discretization by writing the discrete form of Eq. (24) By expressing (v) for each face f using Eq. (17) with f f ∈ S f , we get:

Flow diagnostics impact of dynamic heterogeneity on flooding sweep
One of the best ways to characterize a reservoir rock capacity to store and transmit fluids goes through the deep investigation of its intrinsic properties. Heterogeneity developed by Shook and Mitchell 2009 is the variety of constituent particles encountered in porous rock. In other words, it is the property that defines the reservoir rock permeability and porosity. Studies on the measure of heterogeneity have converged on two models. The first static heterogeneity is all about the rock storage and transmissibility distribution of fluid. Therefore, high static heterogeneity may conduct to a rock good capacity to store and transmit fluid. On the other hand, the dynamic heterogeneity focuses on the measurement of interconnected void space. It describes fluid motion and flow-path connections, thus, a high dynamic heterogeneity may cause a large residence time (time a particle takes to travel from an injector well to a producer well) value. The measure of dynamic heterogeneity is described by the reservoir flow capacity and storage capacity diagram precisely, and it is expressed by the computation of the residence time (sum of backward and forward time of flight). This method has appeared for decades in reservoir engineering, but recently used and upgraded by (Lie 2016). The theory consists to consider the reservoir as a collection of N permeable stream-tube layers. Each layer has no volumetric communication with others; for a flow rate q i , a volume V i , we can define a residence time i = V i ∕q i . Layers are arranged to make their residence time being ascending, 1 ≤ 2 ≤ 3 ⋯ ≤ N , and we assume the displacement within every layer piston types. We can describe the normalized storage capacity i and flow capacity F i as follows.
where i represents the portion of stream-tube layers which have broken through at the residence time i , and F i describes the corresponding volumetric flow portion, and it represents the fraction of the produced injected fluid. We can plot these both data i and F i to get the flow capacity and storage capacity F-ϕ diagram (Fig. 1).
The left plot represents the volume V i , plotted as a function of the flow rate q i ; the blue curve describes a heterogeneous fluid motion and the green shows a homogeneous displacement. On the other hand, the right plot presents the F-ϕ diagrams, where the flow rate and volume have been normalized.

Lorenz coefficient
The characterization of fluid motion between injection and production wells highly implies a perfect knowledge of reservoir dynamic heterogeneity. Used to approximate how much reservoir oil volume can be displaced and extracted from the oilfield. However, the Lorenz coefficient is one of the most popular measures of dynamic heterogeneity, calculated from the F-ϕ diagram. It is the difference in flow capacity from which a piston-type displacement has been applied, defined as follows (Fig. 2).
Graphically, the Lorenz coefficient is two times the area bounded by the F (ϕ) curve (blue) and the F = ϕ line (green) and has a value from (0) to (1). A Lorenz coefficient of 0 falls along the F = ϕ line; green line describes a homogenous fluid motion with an equal volumetric displacement front, and the breakthroughs occur at the same residence times 100% sweep. A value of (1) describes an infinitely heterogeneous fluid motion blue line on the F-ϕ diagram, and the breakthroughs occur at different residence times and lead to stagnant regions.
Flow diagnostics on case studies: SPE9 (Killough 1995) and ZHNBA China oilfield SPE9 case study

SPE9 reservoir model setup
Here, the reservoir setup consists of eight wells in total, six producers and two injectors forming a repetitive well   five spots scheme. All wells are controlled by flow rate, for the producers we set 70 m 3 ∕day and injectors 280 m 3 ∕day . We consider our reservoir initially filled with oil and no active aquifer; then, we set the initial single fluid with the following parameters, = 1 cP and = 1014 Kg/m 3 (Figs. 3 and 4).

SPE9 dynamic heterogeneity and sweep efficiency
After getting the flow field we need for diagnostics, we can compute the F-ϕ quantities and plot the diagram F versus ϕ heterogeneity developed by Shook and Mitchell 2009 to get the dynamic heterogeneity expression. We note that for a homogeneous displacement, the F-ϕ gives a straight line; all flow paths breakthrough at the same time and gives a Lorenz coefficient of value zero. On the other hand, for a heterogeneous flow motion, F-ϕ plot gives a concave curve with the steep initial slope that denotes the high flow regions giving early breakthrough and the flat regions that correspond to the low flow or stagnant regions where there is a breakthrough delay. The Lorenz coefficient takes the values nonzero to one. In the case of SPE9, the Lorenz coefficient is equal to 0.5340. Stepping further, we can define the sweep efficiency diagram which is the measure of the injected fluid volumetric effectiveness. It is the ratio of oil contracted by water at a time (t) and an infinite time (Fig. 5).

SPE9 flow model and volumetric connections
The sum of the backward and forward time of flight gives the residence time, which is the total travel time of particle from an injector well to a producer. This quantity is used to determine high flow zones with small residence time and low flow or stagnant zones with high residence time. The next step of this part after computing the time of flight is the study of the reservoir volumetric sweep. The addition of  injector tracers gives the flooded regions, also amassing the majority of the vote overall tracers, we determine that I2 is the injector that contributes the most (Figs. 6 and 7).
By associating all producer tracers, we can determine the drainage regions, and using the majority of the vote, we determine that P6 is the producer that contributes the most. In addition by combining the flood and the drainage regions, we get the well pair regions. Finally, after studying the drainage and tracer partitions, we can get well pairs and compute the pore volume of regions that are linked with them (Fig. 8).

ZHNBA reservoir model setup
Here, the reservoir setup consists of eleven wells, eight producers and three injectors forming a repetitive well five spots scheme. All wells are controlled by flow rate, for the producers we set 100 m 3 ∕day and injectors 175 m 3 ∕day . We consider the reservoir initially filled with oil and no active aquifer; then, we set the initial single fluid with the following parameters, = 1 cP and = 1014 Kg/m 3 . Before going straight to the computation of rock properties, we first inspect the whole model. We plot the model outlines and identify all faults in red color. Then, we distinguish the active and the inactive reservoir zones. We can see the active part appearing in yellow; we create a new unstructured grid with a non-active part removed. This active model is the basic grid structure we will use for our study (Fig. 9).
Then, we can compute and plot the rock properties, porosity, permeability and wells (Fig. 10).

ZHNBA measure of dynamic heterogeneity and sweep efficiency
We start our flow diagnostics study by computing the F-ϕ diagram for the measure of dynamic heterogeneity, and  we get the Lorenz coefficient evaluation. From the Lorenz coefficient computation, we get Lc = 0.6623 ; this shows that ZHNBA is a highly heterogeneous reservoir (Fig. 11).

ZHNBA flow model and volumetric connections
We compute the sweep region to distinguish the zones in the reservoir which are going to be affected by the injector wells. After identifying the sweep regions, we can compute the drainage to see the zones which are going to be affected by each producer well. We refine partitions by using majority vote to get the sweep and drainage region interactions, and we blend in gray, the zones influenced by a sweep and drainage combined. Also, we compute the well pairs which are the connection between each injection, tracer partitions and pore volumes of regions associated (Fig. 12).

Waterflooding optimization models
Generally, the way to optimize waterflooding requires mathematical methods with several forward simulations; these kind of full and rigorous simulations are costly. Here in this work, we will first present a framework that can be used as a waterflooding optimization proxy. Then secondly, we will optimize the wells flow rate by using flow diagnostics interactive tools. The idea behind the first method is to use flow diagnostics to design a simple and lightweight algorithm easy to implement for waterflooding optimization. On the other hand, the second method is based on the instant TOF snapshot studying. Here, we will see how flow diagnostics with low-cost tools can iteratively improve water injection configurations.

Optimization by using the objective function
In this part inspired by Møyner et al. 2015, our reservoir model is based on the flow field and time of flight (TOF); in the aim to derive a full discretization, we write all the model equations in the residual form: In the above equation, we have the time of flight ( ), Darcy's velocity ( ⃗ u ) and the pressure ( p ) as unknown. The full flow diagnostics is practically one part based on forward and backward time of flight equations; the other part on tracer equations, but in this work, we will precisely examine the forward time of flight equation and simply consider the treatment of the tracer equations as it is analogous.
To create fluid motion within the reservoir, we set n bh wells to operate with bottom hole pressure (bhp) and n r wells operate with rate; from there, we can assume n w = n bh + n r been the total wells having n pf perforations. Our well model is constructed following the Peaceman standard well model per perforated cell; we will have to define the well-perforated fluxes q pf and the bottom hole pressure p bh , and we let n w (j) being the well index of perforation j . n pf (k) represents the indices of perforations of well k , and w j pf represents the Peaceman well index belongings to perforation j ; finally, we can define the well model through the following equation for each perforated cell We need to give precision for each pressure and rate controlled wells Furthermore, to continue our waterflooding optimization, we need to discretize the Lorenz coefficient equation using the standard two-point flux approximation with the aim to get a reservoir simulation model. Herein, the continuous expression of pressure p and the time of flight are changed by vectors p and τ which with their values determined in the center of each cell, and ⃗ u is replaced by the vector fluxes across all cell faces u. Then, p(q, ⃗ u) = 0 becomes P h ( , ) = 0 , with P h the expression of the discretized system. Since in this work, we will only work with the discretized fluid flow equations; we can drop h on the P h , and then we can introduce the discrete model equations as representation of an implicit system which has his linearization coming from the Newton linearized system. Our system will not be constructed explicitly because the pressure and time of flight equations are linked and will be solved sequentially; the system is written as follows: Our approach to obtain an optimal water injection is based on the conception of getting an ideal piston-like displacement that carries us to deeply focus or observe the Lorenz heterogeneity coefficient. The Lorenz coefficient highlights the fluid flow capacity and storage capacity based on piston-like fluid motion; this makes us considering the optimized Lorenz coefficient as our objective function to improve waterflooding performance. Herein, the optimized Lorenz coefficient will seek to equalize the total travel time through the standard two-point flux approximation in all cells and then improve the sweep.

Solution strategy and rate optimization
To optimize waterflooding by using the equalized Lorenz coefficient, we need to focus down on the time of flight and tracers; in particular, we need to solve Eq. (35) and compute the objective function gradient. To this hand, we can start by assuming that we have zero pressure and fluxes because all equations are linear and no iterations are required. We design our solution to aim with the minimum computational workload by solving Eq. (35) in sub-steps. We first assemble the pressure equation and solve it, and then set that the total mobility in Eq. (32) is constant; then, the following linear system after using the two-point flux approximation for the spatial discretization is obtained as follows: The next step after all fluxes computed is the time of flight equations become linear, and then we set 0 = 0 but knowing from Eq. (35) that = , we can get the following simplification Following the same steps, expression (37) can be extended to tracers if required. When we finish with the computation of all forward equations, we then start looking for the objective function gradient, and to carry this investigation, we can use two ways, applying a numerical differentiation or adjoint method.
In this work, we use the adjoint method although its application requires alteration of the model equations; it is more efficient than a numerical differentiation which is simple to implement and non-intrusive but can be computationally costly because all forward equations need to be solved on each perturbation.
We begin the adjoint model by introducing the equations derivation tool the stacked vector of solution quantities × T = (p T , u T , T , q T , p T bh ) ; here, each value depends on the control w which is the target value of BHP or target rate of each well as expressed by Eqs. (34). Our objective function is formulated as follows: G[×(w)] and stands as a function for solution equation quantities. The target here is to find a set of constraints g [×(w), w] = 0 that represent the reservoir Eqs. (32) and (34). To this hand, we define the Lagrange function to describe our problem function that has the similitude values with our objective function and written as: Since we aim to find our objective function derivative, we can differentiate to get: The simplification of the last equation can be done by elimination of all terms which contain d× dw and defining (solving) by the adjoint equation: We can furthermore simplify if we consider a state x that makes g[×(w), w] = 0 , and then we get: Going deeply, our objective function will not depend on the control w exclusively but on state variable x ; in particular, the term G w can disappear. We can conclude this solution strategy study by highlighting the most important steps which are first, solving the forward system Eq. (35) using the block-wise approach explained previously to get state variables and upwind directions, and secondly solving the backward problem Eq. (40) with the same approach to get the objective function sensitivities.
After designing our equations solutions, we set a pathway to optimize well rates by using the steepest descent algorithm to get an optimal control; herein, our rate optimization experiment is assimilated to the control variable w which precisely represents well rates. In our approach production wells are not designed to become injectors and vice versa, also the injection rate is kept constant. All these requirements are sustained by the following equation: Here, represent the step size. Our target in the last equation is to describe a stationary point in the aim to get the gradient equal zero. To switch our rate optimization, we proceed iteration in the steepest descent algorithm to evaluate the objective function for the update controls; in the case we have no improvement, we reduce the step size until we get an improvement. When the projected gradient small enough and the objective function improvement between two iterations is smaller than the prescribed tolerance, the algorithm stops.

Application on SPE9
We first define the objective function which is computed based on the grid cells, the rock properties (permeability and porosity) and the equalized Lorenz coefficient. Then, we deduct the minimum rate for each cell and launch the algorithm that targets to improve the time of flight and tracers, and the waterflooding well rate (Fig. 13). We compute the optimized Lorenz coefficient and plot the new F-diagram and the reservoir sweep efficiency (Fig. 14).
In Tables 1 and 2, we observe the initial and the optimized well configuration. Here, we can clearly notice the change in injector rate which in the optimized well configuration appeared with their optimal values. We convert these values from the flow field unity to cubic meter per day and put them in a table with the initial and optimized Lorenz coefficient (Table 3).

Application on ZHNBA China oilfield
Before starting the waterflooding optimization process, we first launch the calling algorithm of our flow diagnostic interactive tools named GUI, and it appears with its input control screen and plot compartment showing our reservoir model (Fig. 15). We start the waterflooding optimization experiment by watching every displacement front snapshot of flooding volumes based on residence time flow from injection wells. We set as base case our ZHNBA initial setup; all eleven wells are controlled by rate as follows: the three injectors 175 m 3 ∕day and producers 100 m 3 ∕day . We run the base case injectors displacement front simulation to notify early breakthroughs in producer wells; then, we rank and optimize the injection well rates in the aim to equalize all producer breakthroughs. In our flow experiment study at the base case, the producers ranking breakthrough is identified as follows (Fig. 16): • The producer P6 records the earliest water breakthrough from the injector I1 • Following by the producers, P8 and P7 get their breakthrough from I3 • The remaining producer water breakthroughs are considered as late; then, they are not recorded including those coming from the injector I2.
From this flow experiment, we can elaborate our waterflooding optimization strategy to harmonize the producer breakthroughs, and then decreasing the Lorenz coefficient and optimizing the sweep. To this hand, we first increase the rate of I2 from 175 to 250 m 3 ∕day keeping the other well rates at the base case setup. We run the flow diagnostics simulation and compute the Lorenz coefficient which is decreasing; we call this step case one. Case 2: we increase again all the injector rates; we assign I1 to 200 m 3 ∕day , I2 to 300 m 3 ∕day and I3 to 200 m 3 ∕day ; we run the simulation and compute the Lorenz coefficient. Then, we restart the same process for the last time in case 3; we increase the rate of injector I1 to 210 m 3 ∕day and injector I3 to 250 m 3 ∕day ; we run the simulation and compute the Lorenz coefficient. We end the optimization process by plotting all flow capacity and storage capacity diagrams and arranging in a table the injector wells optimal rates (Table 4 and Fig. 17).

Conclusion and summary
Waterflooding is one of the best methods used to improve oil recovery and has sown its efficiency in many applications. However, herein, we have presented a pathway to improve its performance using finite volume discretizationbased flow diagnostics simulation. By drawing its sweep efficiency improvement scheme based on the optimization of the dynamic heterogeneity Lorenz coefficient, we harmonized the flooding breakthroughs. Designing a waterflooding optimization performance framework using flow diagnostics-based finite volume method is having the reservoir flow quantities fluxes computed on volumetric average portioning. In this work, the finite volume discretization showed its perfect ability to handle both Cartesian and unstructured grids, when flow diagnostic tools deepen our understanding of reservoir flow characterization summarized in two parts: firstly, the solution of the incompressible steady-state flow equation through the standard two-point flux approximation to get the fluid bulk motion. Secondly, the solution of the mixture time of flight and tracer equations to split the reservoir into volumetric flow and describing the flow path regions.