Simulating Multi-phase Flow in Fractured Reservoirs

Abstract

A method for performing fluid extraction from a fractured subsurface formation includes receiving a discrete fracture model representing the fractured subsurface formation and receiving pressure values and saturation values for multiple fluid phases across the discrete fracture model. Based on the pressure values for the multiple fluid phases across the discrete fracture model, face-centroid velocities are generated for the cells and the pressure values for the multiple fluid phases are updated by performing operations including a mimetic finite difference analysis. Based on the generated face-centroid velocities, an exit face and time-of-flight is determined for each cell and the saturation values are updated for the multiple fluid phases across the discrete fracture model based on the exit face and time-of-flight for each cell.

Description

TECHNICAL FIELD

[0001] This specification relates to simulating multiphase flow in subsurface formations, particularly in fractured reservoirs.BACKGROUND

[0002] The modeling of multi-phase flow in fractured formations is used in various energy and environmental applications such as geological carbon dioxide (CO2) sequestration, geothermal energy extraction, radioactive waste management in the subsurface, oil and gas recovery process in fractured reservoirs. Detailed characterization of fractured formations is commonly described by a discrete-fracture model (DFM), in which the fractures and rock-matrix are explicitly represented by unstructured grid elements.

[0003] Various numerical methods based on DFM have been used to simulate multi-phase flow in fractured formations. Examples include Galerkin finite element (FE), cell centered finite volume (FV), control volume FE, mixed FE, boundary element, and mimetic finite difference (MFD) methods.SUMMARY

[0004] This specification describes an approach to simulating multiphase flow in subsurface formations, particularly in fractured reservoirs. The systems and methods of this approach use a streamline-based simulator based on DFM, for multi-phase flow in fractured reservoirs. These systems and methods provides a hybridize mimetic finite difference and streamline (MFD-SL) approach, in which mimetic finite difference (MFD) is employed to discretize the pressure equation, and the streamline (SL) method is adopted to solve the saturation equation along 1D streamlines. The hybrid formulation is implemented on a discrete fracture model (DFM) and implicitly addresses pressure and explicitly addresses saturation. This approach provides a simple and practical streamline tracing method applicable to the triangular and tetrahedral grids that are commonly implemented in DFMs. This approach has been demonstrated to provide good results with lower computation costs than earlier formulations. These improvements are provided in part by using streamline tracing, which reduces the numerical diffusion error and improves the computation efficiency.

[0005] The approach disclosed in this specification can be used to simulate multiphase flow in fractured reservoirs with improved calculation accuracy, computation efficiency, and general applicability, especially for reservoir-scale simulations. In particular, this approach can be used to help plan, manage, and control field operations (e.g., planning and controlling operation of injection and production pumping during hydrocarbon recovery).

[0006] The use of MFD provides the general applicability by enabling polygonal grids of any shape with full-tensor permeability. Additionally, DFMs are applicable to more general flow scenarios than embedded discrete fracture models (EDFM) in which fracture and matrix grids are constructed independently and then coupled to each other via source / sink relations.

[0007] The MFD-SL approach also provides improved computational efficiency. The implementation of MFD only requires one cell for the calculation of numerical flux. It just requires coordinates information rather than needing extra basis functions and integration. In addition, the use of SL tracking significantly improves computational efficiency compared to conventional approaches.

[0008] The MFD-SL approach also provides improved accuracy. The DFM enables the accurate description of matrix-fracture interaction involved more complex physics. In addition, the streamline tracking method allows accurate description of flow within fracture grids.

[0009] The details of one or more embodiments of the invention are set forth in the accompanying drawings and the description below. Other features, objects, and advantages of the invention will be apparent from the description and drawings, and from the claims.

[0010] The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee.DESCRIPTION OF DRAWINGS

[0011] FIG. 1 is a schematic view of field operations being performed to map subsurface features in and produce hydrocarbons from a subsurface formation.

[0012] FIG. 2 is a flowchart illustrating a method for simulating multiphase flow in subsurface formations.

[0013] FIG. 3A is a schematic illustrating two fracture cells and FIG. 3B is a schematic illustrating multiple neighboring fracture cells.

[0014] FIG. 4 is a schematic illustrating streamline tracing in a two-dimensional triangular matrix cell i.

[0015] FIG. 5 is a flowchart of streamline tracing in 2D triangular matrix cell.

[0016] FIG. 6A is a schematic illustrating streamline tracing in multiple fracture cells. FIG. 6B is a schematic illustrating streamlines flowing out from upstream fracture cells fi to downstream ones. FIG. 6C is a schematic illustrating streamlines flowing into downstream fracture cells fi from upstream ones.

[0017] FIG. 7 is a schematic illustrating streamline tracing in triangular prism.

[0018] FIGS. 8A-8D illustrate a fracture model with a single diagonal fracture and four meshes in comparing the current approach with other approaches.

[0019] FIGS. 9A and 9B illustrate the streamlines generated based on the fracture model with the single diagonal fracture and the coarse mesh. FIGS. 9C-9H present numerical results for comparison.

[0020] FIGS. 10A and 9B illustrate the streamlines generated based on the fracture model with the single diagonal fracture and the fine mesh. FIGS. 9C-9H present numerical results for comparison.

[0021] FIG. 11A illustrates a reservoir model with two intersecting highly conductive features and a flow barrier. FIGS. 11B and 11C illustrate two meshes used in comparing the current approach with other approaches.

[0022] FIGS. 12A and 12B illustrate the streamlines generated based on the reservoir model with two intersecting highly conductive features and a flow barrier. FIGS. 12C-12H present numerical results for comparison.

[0023] FIG. 13A illustrates a reservoir model with a complicated fracture network. FIG. 13B presents an unstructured mesh for the reservoir model. FIGS. 14A-14C present the properties, respectively, fracture permeability, fracture porosity, and fracture aperture, assigned to the fractures.

[0024] FIGS. 15A and 15B illustrate the streamlines generated based on the reservoir model of FIG. 13A. FIGS. 15C-15F present numerical results for comparison.

[0025] FIG. 16 illustrates hydrocarbon production operations.

[0026] Like reference symbols in the various drawings indicate like elements.DETAILED DESCRIPTION

[0027] This specification describes an approach to simulating multiphase flow in subsurface formations, particularly in fractured reservoirs. The systems and methods of this approach use a streamline-based simulator based on DFM, for multi-phase flow in fractured reservoirs. This approach hybridizes mimetic finite difference and streamline (MFD-SL) methods, in which mimetic finite difference (MFD) is employed to discretize the pressure equation, and the streamline (SL) method is adopted to solve the saturation equation along 1D streamlines. The hybrid formulation is implemented on discrete fracture model (DFM) and implicitly addresses pressure and explicitly addresses saturation. This approach provides a simple and practical streamline tracing method applicable to the triangular and tetrahedral grids that are commonly implemented in DFMs. This approach has been demonstrated to provide good results with lower computation costs than earlier formulations. These improvements are provided in part by using streamline tracing, which reduces the numerical diffusion error and improves the computation efficiency.

[0028] FIG. 1 is a schematic view of field operations being performed to map subsurface features in and produce hydrocarbons from a subsurface formation 100. These field activities provide the underlying basis for implementation of the systems and methods described with reference to FIG. 2.

[0029] The subsurface formation 100 includes a layer of impermeable cap rock 102 at the surface. Facies underlying the impermeable cap rocks 102 include three other layers 104, 106, and 108. A fault line 110 extends across the layer 104 and the layer 106.

[0030] Oil and gas tend to rise through permeable reservoir rock until further upward migration is blocked, for example, by the layer of impermeable cap rock 102. Seismic surveys attempt to identify locations where interaction between layers of the subsurface formation 100 are likely to trap oil and gas by limiting this upward migration. For example, FIG. 1 shows an anticline trap 107, where the layer of impermeable cap rock 102 has an upward convex configuration, and a fault trap 109, where the fault line 110 might allow oil and gas to flow in with clay material between the walls traps the petroleum. Other traps include salt domes and stratigraphic traps.

[0031] A seismic survey is being performed using a seismic source 112 (for example, a seismic vibrator or an explosion) that generates seismic waves that propagate in the earth. Although illustrated as a single component in FIG. 1, the source or sources 112 are typically a line or an array of sources 112. The generated seismic waves include seismic body waves 114 that travel into the ground and seismic surface waves 115 travel along the ground surface and diminish as they get further from the surface.

[0032] The seismic body waves 114 are received by a sensor or sensors 116. Although illustrated as a single component in FIG. 1, the sensor or sensors 116 are typically a line or an array of sensors 116 that generate an output signal in response to received seismic waves including waves reflected by the horizons in the subsurface formation 100. The sensors 116 can be geophone-receivers that produce electrical output signals transmitted as input data, for example, to a computer 118 on a seismic control truck 120. Based on the input data, the computer 118 may generate a seismic data output, for example, a seismic two-way response time plot.

[0033] The seismic surface waves 115 travel more slowly than seismic body waves 114. Analysis of the time it takes seismic surface waves 115 to travel from source to sensor can provide information about near surface features.

[0034] A control center 122 can be operatively coupled to the seismic control truck 120 and other data acquisition and wellsite systems. The control center 122 may have computer facilities for receiving, storing, processing, and analyzing data from the seismic control truck 120 and other data acquisition and wellsite systems that provide additional information about the subsurface formation. For example, the control center 122 can receive data from a computer 119 associated with a well logging unit 121.

[0035] The computer systems 124 can be located in a different location than the control center 122. Some computer systems are provided with functionality for manipulating and analyzing the data, such as performing seismic interpretation or borehole resistivity image log interpretation to identify geological surfaces in the subsurface formation or performing simulation, planning, and optimization of production operations of the wellsite systems.

[0036] A wellbore 130 that has been drilled in the subsurface formation 100 is being logged in a well logging operation 128. The wellbore 130 extends downhole from a wellhead 132. The wellbore 130 is a vertical wellbore but well logging can also be performed in other wellbores, for example, slanted or horizontal wellbores. In the well logging operation 128, the wellbore 130 penetrates through three layers 102, 104, and 106 of a subsurface formation 100. A control truck 121 lowers a logging tool 134 down the wellbore 130 on a wireline 136. The control truck 121 can provide data generated by well logging to the control center 122.

[0037] The computer systems 124 in the control center 122 can be configured to analyze, model, control, optimize, or perform management tasks of field operations associated with development and production of resources such as oil and gas from the subsurface formation 100. For example, an injection well 123 and a production well 125 extend into layer 104 of the subsurface formation 100. Based on data gathered by the exploratory field operations, the computer systems 124 can generate models such as a reservoir model for portions of the subsurface formation 100. These models can simulate the effects of production field operations (e.g., injecting water or carbon dioxide through the injection well 123 to increase the production of hydrocarbons through the production well 125). The simulations can be used to plan and, in some instances, control field operations (e.g., the operation of pumps associated with the injection well 123 and the production well 125).

[0038] FIG. 2 is a flowchart illustrating a method 200 for simulating multiphase flow in subsurface formations. As previously noted, these simulations can be used to plan and / or control field operations. Although developed for use to describe flow in fractured subsurface formations and described below with respect to a fractured subsurface formation, the method 200 can also be applied to subsurface formations that are not significantly impacted by fractures.

[0039] The method 200 requires data from a fractured subsurface formation. It also requires a discrete fracture model representing the fractured subsurface formation and having matrix cells and fracture cells. In some instances, the method 200 includes generating data from the fractured subsurface formation (step 210) and generating the discrete fracture model (step 212). The data can be gathered, for example, the technologies described with respect to FIG. 1. Alternatively, stored data representing the results of earlier exploratory field operations can be used. Similarly, a previously generated discrete fracture model can be used rather than generating a new discrete fracture model. The method 200 includes receiving the discrete fracture model whether it was generated as part of the method 200 or

[0040] The flow of incompressible and immiscible phases in porous media can be described by Darcy's law and the saturation equation. For a two-phase system with oil and water, the governing equations of flow are described by the saturation equation and the generalized Darcy law of each phase. By ignoring capillarity and saturation constraints leads to the formulation:ϕ⁢∂Sα∂t+∇ (vα⇀)=qα⁢ α=o⁢ or⁢ w(1)vα⇀=-kr⁢αμα⁢K↔⁢∇pα⁢ α=o⁢ or⁢ w(2)So+Sw=1(3)where ϕ is the porosity, t is the time, Sα, , qα, krα, μα and pα are the saturation, velocity, sink / source term, relative permeability, viscosity, and pressure of phase α, respectively, is the absolute permeability tensor. Eq. (1) and Eq. (2) can be rewritten in equivalent forms as:-∇ (λ⁢K↔⁢∇p)=q(4)ϕ⁢∂sw∂t-∇ (λw⁢K↔⁢∇p)=qw(5)whereλα=kr⁢αμαis the mobility of phase α, λ=λo+λw is the total mobility, and q=qo+qw is the total sink / source term.Eq. (4) and Eq. (5), together with suitable initial and boundary conditions, form a global equation system, which is solved by a hybrid MFD-SL formulation. In the method 200, pressure values and saturation values received for multiple fluid phases across the discrete fracture model (step 216) provide the initial conditions.The global equation system is decoupled and solved sequentially in an approach, in which velocities are generated from the pressure equation as discretized by MFD (step 218), and the saturation equation is solved by a streamline-based method, which transforms the advection-dominated transport problem into 1D flow problem to be solved along streamlines (step 220). The hybrid MFD-streamline is implemented based on the DFM, where fractures are geometrically simplified by using (n−1)-dimensional grid cells in an n-dimensional matrix domain. The calculated pressure and saturation values are stored in the time counter is incremented (step 222). If the time counter has not reached its end value, steps 220-222 are repeated based on the stored pressure and saturation values (step 224). If the time counter has reached its end value, the simulation run is complete. In some cases, the method 200 includes using results of the simulation to plan and / or control field operations (step 226).The MFD approximation of volumetric flux is based on integrating Eq. (4) on control volume of cell i (denoted by Ωi), which yields:-∫ Ωi∇ (λ⁢K↔⁢∇p)⁢ d⁢Ω=Q(6)in which, Q=∫Ω<sub2>i < / sub2>qdΩ. Applying the divergence theorem on the left side of Eq. (6) provides:-∑ β=1N⁢∫∂Ωiβλ⁢K↔⁢∇p⁢n⇀⁢dA=∑ β=1N⁢Fiβ=Q(7)in which, N is the number of faces of cell i, and Fi<sub2>β< / sub2> is the volumetric flux across one face of cell i (referred to by iβ). From this, the volumetric flux of Fi<sub2>β < / sub2>based on MFD is given by:Fiβ=λiβ⁢∑ γ=1N⁢Tiβ,iγ(pi-piγ)(8)in which, λi<sub2>β< / sub2> is the total mobility at the centroid of face iβ, calculated by single-point upstream scheme, represents the transmissibility matrix of cell i, pi is the pressure at the centroid of cell i, and pi<sub2>γ< / sub2> is the pressure at the centroid of face iγ. The transmissibility tensor can be obtained by:T↔i=1<semantics definitionURL="">❘<annotation encoding="Mathematica">"\[LeftBracketingBar]"< / annotation>< / semantics>Ωi<semantics definitionURL="">❘<annotation encoding="Mathematica">"\[RightBracketingBar]"< / annotation>< / semantics>⁢N⇀i⁢Ki↔⁢N⇀iT+6d⁢t⁢r⁡(Ki↔)⁢Ai↔(Ii-Qi⇀⁢Qi⇀T)⁢Ai↔(9)in which, =(|σΩi<sub2>1< / sub2>|, . . . |σΩi<sub2>β< / sub2>|, . . . , |σΩi<sub2>N< / sub2>|)T, |σΩi<sub2>β< / sub2>| is the area of face iβ, is the outward unit normal on face iβ, is the absolute permeability tensor of cell i, =diag (|σΩi<sub2>1< / sub2>|, . . . |σΩi<sub2>N< / sub2>|), =orth(), is the vector from cell centroid to face centroid.The matrix cells and fracture cells require different discretization approaches.For matrix cells, Eq. (8) is introduced into Eq. (7) to yield the discretized form of Eq. (6) for matrix cell i (referred to by superscript m) as:∑ β=1Nm⁢λiβm⁢∑ γ=1Nm⁢Tiβ,iγm(pim-piγm)=Qm(10)The following flux continuity conditions are imposed on the interface (i.e., iβ=jη) between two neighboring cells i and j to link the matrix cell i with the neighboring cell j. If j is a matrix cell:Fiβm+Fjηm=0(11)If j is a fracture cell (referred to by superscript f):Fiβm+Fjηf=0(12)For fracture cells, the discretized form of Eq. (6) for fracture cell i is:∑ β=1Nf⁢λiβf⁢Tiβf(pif-piβf)=Qf(13)The transmissibility of fracture cell i is calculated based on two-point flux approximation (TPFA) scheme due to the K-orthogonal feature of fracture cells. The flux continuity conditions are different depending on whether there are two neighboring fracture cells or there are more than two fracture cells having an intersection point.FIG. 3A is a schematic illustrating two fracture cells and FIG. 3B is a schematic illustrating multiple neighboring fracture cells. For two neighboring fracture cells (denoted by i and j) (see FIG. 3A), the flux continuity conditions are:Fiβf+Fjηf=0(14)For multiple (more than two) neighboring fracture cells (denoted by i, j, k, l) with an intersection at point o (see FIG. 3B), the sets of fracture cells are determined, such that the flux qx,of=Tx<sub2>62 < / sub2>f(pxf−po)x ∈{i, j, k, l}:qxupstream,of≥0⁢ xupstream∈x(15)qxdownstream,of<0⁢ xdownstream∈xin which, xupstream and xdownstream represent the upstream and downstream fracture cells, respectively. The flux continuity condition proposed by Hoteit and Firoozabadi (2005) is adopted as follows:∑I∈xupstreamλI⁢qI,of+∑I∈xdownstreamλI⁢qI,of=0(16)λI={λI,oI∈xupstreamλoI∈xd⁢o⁢wnstreamCombining Eq. (10) and Eq. (13) for all matrix and fracture cells, together with corresponding flux continuity equations in Eq. (11), Eq. (12), Eq. (14), and Eq. (16) completes the pressure equation system, which is implicitly solved by the MFD scheme which provides pressures at centroids of cells and faces are obtained. Further, velocities at centroids of faces can be calculated as:Viβ⇀=Fiβ<semantics definitionURL="">❘<annotation encoding="Mathematica">"\[LeftBracketingBar]"< / annotation>< / semantics>∂Ωiβ<semantics definitionURL="">❘<annotation encoding="Mathematica">"\[RightBracketingBar]"< / annotation>< / semantics>(17)in which, Fi<sub2>β< / sub2> is calculated by using Eq. (8).Based on the generated face-centroid velocities, an exit face and time-of-flight for each cell is determined and the saturation values for the multiple fluid phases are updated across the discrete fracture model based on the exit face and time-of-flight for each cell.FIG. 4 is a schematic illustrating streamline tracing in 2D triangular matrix cell i. The illustrated matrix cell i has the node coordinates Nβ, β∈{1, 2, 3} and the outward unit normal vectors ={, , }T. Given face-centroid velocities =(, , )T obtained by Eq. (17) and assume constant velocity distributed on the whole triangular domain denoted by =(Vi-x, Vi-y)T, the following overdetermined linear system can be formulated:Gi⇀⁢Vi⇀=Viβ⇀(18)with corresponding least-square solutions as:Vi⇀=(Gi⇀T⁢Gi⇀)-1⁢(Gi⇀T⁢Viβ⇀)(19)The assumption of constant velocity distributed on the triangular domain makes straight-line streamlines, which simplifies the complexity of streamlines tracing. The validity of this assumption depends on the residual error of Eq. (19). It is noted that for an overdetermined equation system, the greater number of equations, the bigger expected residual errors and the number of equations correspond to the number of faces of polygonal cells. The application of quadrilateral cells is supposed to honor bigger residual error than triangular cells. This is the case for 3D cases. For the examples discussed in this specification, triangular and tetrahedral cells are implemented for 2D and 3D cases, respectively, instead of other type of cells.Based on the node coordinate of the streamline inlet Ninlet=(xinlet, yinlet), the streamline equation can be parameterized as:A=Ninlet+Δ⁢t⁢V⇀i ⁢ Δ⁢t>0(20)Which face the streamlines exit from and the corresponding time of flight is then determined. The parametric equation of face β is given:B=Nβ+Δ⁢r⁡(Nγ-Nβ)⁢ Δ⁢r>0(21)Assuming streamlines exit from the face β: WP=U(22)in which, W=(, Nβ-Nγ), P=(Δt, Δr)T, U=(Nβ-Ninlet).FIG. 5 is a flowchart of streamline tracing in 2D triangular matrix cell. If the matrix W in Eq. (22) is singular, it indicates streamlines are parallel to the face β (coincidence occurs in the case of streamlines entering from vertices of the face β). A new face is then evaluated. If W is nonsingular, unique solutions of (Δt, Δr) exist for Eq. (22). If Δt<0, it means the streamlines enter the matrix cell from the face β or intersect with the face β in the opposite direction of velocity . Face β, thus, is not the face where streamlines exit from. Followed by a new face is evaluated. If Δt>0, in the cases of Δr<0 or Δr>1, it indicates streamlines intersect with the extension of the face β. A new face is then evaluated. In the cases of 0≤Δr≤1, streamlines intersect with the face β, from which the streamlines exit. Time of flight (TOF) is given by: TOF=ϕ⁢ Δ⁢t(23)Accordingly, the node coordinates of the streamline outlet can be determined by:Noutlet=Nβ+Δ⁢r⁡(Nγ-Nβ)(24)Difficult streamline tracing occurs in cases where the distance between the streamline outlet and vertices of the intersecting face is too small. To avoid this, a threshold parameter denoted by rd (e.g., 0.1) is introduced based on which Δr is determined:0≤Δ⁢r≤rd⁢ Δ⁢r=rd(25)rd<Δ⁢r<1-rd⁢ Δ⁢r=Δ⁢r1-rd≤Δ⁢r≤1⁢ Δ⁢r=1-rdStreamlines tracing on 3D tetrahedral matrix cells can follow the same procedure of 2D triangular case.FIG. 6A is a schematic illustrating streamline tracing in multiple fracture cells. FIG. 6B is a schematic illustrating streamlines flowing out from upstream fracture cells fi to downstream ones. FIG. 6C is a schematic illustrating streamlines flowing into downstream fracture cells f1 from upstream ones.For the 2D discrete fracture model, fracture cells in the gird domain are geometrically simplified by 1D lines, but they are considered as rectangular cells with corresponding lengths and widths (i.e., apertures) in streamlining tracing. The classical Pollock method can be directly implemented with multiple neighboring fracture cells (denoted by i, j, k, l) with corresponding apertures referred to by wif, wjf, wkf, wlf, in which i and k are assumed to be upstream fracture cells while j and l are downstream ones (see FIG. 5A). The face of the fracture cell at the intersection o is denoted by Γx,of with x ∈{i, j, k, l}.Taking upstream fi (see FIG. 6B), streamlines from upstream fi to all downstream cells are determined. Accordingly, Γi,of is divided into two sub-faces (referred to by df<sub2>i→l < / sub2>and df<sub2>i→j< / sub2>) corresponding to areas from which streamlines flow out of upstream fi. The sub-faces are distributed from the leftmost of the face Γi,of following the order of downstream fracture cells in a clockwise manner.dfi→l=wif⁢qflqfl+qfj⁢dfi→j=wif⁢qfjqfl+qfj(26)In terms of downstream fl (see FIG. 6C), streamlines flow into downstream fl from all upstream cells are determined. Accordingly, Γl,of is divided into two sub-faces (referred to by df<sub2>k→l < / sub2>and df<sub2>i→l< / sub2>) corresponding to areas from which streamlines flow into downstream fl. The sub-faces are distributed from the rightmost of the face Γl,of following the order of upstream fracture cells in a clockwise manner.dfk→l=wlf⁢qfkqfk+qfi⁢dfi→l=wlf⁢qfiqfk+qfi(27)The classical Pollock method is then applied to these pre-processed upstream and downstream fracture cells for streamline tracing.FIG. 7 is a schematic illustrating streamline tracing in triangular prism. For the 3D discrete fracture model, fracture cells are treated as triangular prism with height equal to its apertures wf. In this case, the classical Pollock method is inadequate for streamlining tracing. A modified Pollock method is then developed with the following assumptions: linearly distributed velocity in the direction perpendicular to the fracture plane (i.e., z direction) and constant velocity on the triangular fracture planes (x-y planes). The specific procedure is detailed as follows.According to the first assumption, travel time in the z-direction can be calculated by following similar steps in the classical Pollock method:Δ⁢tz=1mz⁢ln⁢vz,outletvz,inlet(28)The streamlines traveling in the x-y planes is same as that in 2D triangular cells. Thus, based on Eq. (22), travel time in the x-y plane (denoted by Δtx-y) and uniform velocity on the x-y plane (denoted by =(Vx, Vy)) can be obtained.Node coordinates of streamline outlet can also be obtained. If Δtz<Δtx-y, streamline exits from the triangular face and node coordinates of streamline outlet can be obtained by:xoutlet=xinlet+Δ⁢tz⁢Vx(29)youtlet=yinlet+Δ⁢tz⁢Vyzoutlet={wfvz,o>00vz,o=0If Δtx-y<Δtz, the streamline exits from the rectangular face, in which the classical Pollock method can be employed for calculate zoutlet and outlet coordinates in x and y directions can be obtained by Eq. (24).The water saturation is calculated by streamline-based method, which transforms the advection-dominated transport problem into 1D transport equation to be solved along each streamline. The 1D transport equation along the streamline is given by:ϕ⁢∂Sw∂t+v⇀⁢∂fw∂l=0(30)in which fw is the water fractional function, is the total velocity along the streamline, and l is the distance along the streamline.The TOF required to travel the distance l along the streamline can be The TOF required to travel the distance l along the streamline can be computed as:τ=∫0 lϕv⇀⁢ ds(31)Differentiating Eq. (31) and Eq. (30) can be rewritten in equivalent form as:∂Sw∂t+∂fw∂τ=0(32)Based on Eq. (32), water saturation can be calculated as long as the time of flight along each streamline is determined. In a prototype whose results are described later, the upwind finite difference scheme was adopted to solve the Eq. (32) with the discretized form as:Sw,in+1=Sw,in-Δ⁢tsΔ⁢τ⁢(fw,in-fw,i-1n)(33)in which Δts is the local time step, Δτis the time-of-flight step, which honors the CFL condition to ensure regular spaced τgirds. It is noted that the limit of CFL condition on Δτ only depends on the selected Δts regardless of cell size, which differs from the traditional IMPES approach.The cell-average water saturation of the cell i can be calculated as:Sw,i_=∑ η∈ΛiSw,iη⁢Δ⁢τiη∑ η∈ΛiΔ⁢τiη(34)in which, Λi represents the set of streamlines traveling through cell i and iη denotes the η th streamline within cell i with the corresponding time-of-flight interval Δτi<sub2>η< / sub2>.A prototype system implementing this approach was developed was developed and used to compare the accuracy and computational time of this approach to earlier methods, including a full MFD scheme on DFM and EDFM. Three numerical examples with varying complexities are provided. Example 1 shows the accuracy of the hybrid MFD-SL on a simple reservoir model. Example 2 illustrates its flexibility in handling flow barriers (i.e., permeability is set to be 0). Example 3 illustrates the effectiveness of the hybrid MFD-SL in complex fracture networks with heterogeneous fracture properties. In the three examples, linear relative permeability with the total mobility equal to one was applied. Under this condition, the water fractional function (fw) is a linear function of water saturation (Sw). Accordingly, the solutions exhibited a discontinuous water-flooding front consistent with the Buckley-Leverett theory. Therefore, the bandwidth of the water-flooding front reflects the magnitude of numerical diffusion error, i.e., the narrower the bandwidth, the smaller the numerical diffusion error. CPU time was provided to evaluate computation efficiency for all tested methods.Example 1FIGS. 8A-8D illustrate a fracture model with a single diagonal fracture and four meshes in comparing the current approach with other approaches. This example shows the accuracy of the hybrid MFD-SL formulation on a simple reservoir model with one single diagonal fracture. The corresponding unstructured (see FIGS. 8A and 8B) and Cartesian meshes (see FIGS. 8C and 8D) were generated based on DFM and EDFM, respectively, with almost the same number of cells for coarse and fine schemes. The domain was initially saturated with oil, with an injector well at the lower left corner and a producer well at the opposite corner. Relevant data of Example 1 is given in Table 1.TABLE 1Relevant data for Example 1Domain dimensions:100 m × 100 mMatrix properties:Km = 10 mD  ϕm = 0.2Fracture properties:Kf = 105 mD   af = 1 cm  ϕf = 0.4Fluid properties:μo = μw = 1 cpRelative permeabilities:linear relationResidual saturations:Srw = Sro = 0.15Well settings:Q inj=Q pro=2⁢0⁢m3dayFIGS. 9A and 9B illustrate the streamlines generated based on the fracture model with the single diagonal fracture and the coarse mesh. FIGS. 9C-9H present numerical results for comparison. The results are at 600 and 800 days for the hybrid MFD-SL (see FIGS. 9C and 9D), full MFD based on DFM (see FIGS. 9E and 9F), and full MFD based on EDFM (see FIGS. 9G and 9H) on the coarse meshes. The streamline distribution obtained using the hybrid MFD-SL method (see FIGS. 9A and 9B) corresponds to the oil saturation in FIGS. 9C and 9D.FIGS. 10A and 10B illustrate the streamlines generated based on the fracture model with the single diagonal fracture and the fine mesh. FIGS. 10C-10H present numerical results for comparison. The results are at 600 and 800 days for the hybrid MFD-SL (see FIGS. 10C and 10D), full MFD based on DFM (see FIGS. 10E and 10F), and full MFD based on EDFM (see FIGS. 10G and 10H) on the fine meshes. The streamline distribution obtained using the hybrid MFD-SL method (see FIGS. 10A and 10B) corresponds to the oil saturation in FIGS. 10C and 10D.It can be observed that the hybrid MFD-SL method significantly reduced the numerical diffusion errors with the smallest bandwidth of the water-flooding front compared to full MFD schemes on DFM and EDFM. As the number of cells increases, the numerical diffusion errors of full MFD schemes on DFM and EDFM show a decreasing trend but still remain larger than those of the proposed method.In addition, the proposed method provides higher computation efficiency. Table 1, as shown in Table 2. As a result, the hybrid MFD-SL method provides a twofold advantage over full MFD schemes based on DFM and EDFM, excelling in both accuracy and computation efficiency.TABLE 2CPU time in seconds for all tested methods in Examples 1 and 2.Example 1Example 1(Coarse)(Fine)Example 2Hybrid MFD-SL298059Full MFD on DFM94724538Full MFD on EDFM624074195Example 2: Two Intersecting Fractures and One Flow BarrierFIG. 11A illustrates a reservoir model with two intersecting highly conductive features and a flow barrier. FIGS. 11B and 11C illustrate two meshes used in comparing the current approach with other approaches.This example illustrates the flexibility of the hybrid MFD-SL in dealing with complex geological settings with a reservoir model with two intersecting highly conductive fractures (represented by blue lines) and one flow barrier (red line), as illustrated in FIGS. 11A-11C. The injector is positioned at the centroid of the lower side, with two producers situated at the two corners of the upper side. All other relevant data was adopted from Example 1, except the production rates for both producers are set to be10⁢ m3day.FIGS. 12A and 12B illustrate streamline distribution. FIGS. 12C-12H show oil saturation profiles for the hybrid MFD-SL method (see FIGS. 12C and 12D), full MFD based on DFM (FIGS. 12E and 12F), and full MFD based on EDFM (FIGS. 12G and 12H).It can be observed that EDFM fails to properly capture the water-flooding front compared to DFM, especially at 200 and 300 days when the water-flooding front encounters the flow barrier. This highlights the inherent limitations of EDFM in simulating two-phase flow across low-conductive fractures, especially when dealing with flow barriers. Hence, streamline simulator based on EDFM limits its applicability in specific flow scenarios.The proposed method also provides more accurate results with smaller numerical diffusion errors than the full MFD based on DFM. Furthermore, the proposed method enables an accurate description of the shape and distribution of the water-flooding front with more details included. The CPU time for all tested methods was previously summarized in Table 2. The hybrid MFD-SL formulation surpasses other classical methods in applicability, accuracy, and efficiency.Example 3: Complex Fracture Networks

[0089] FIG. 12A illustrates a reservoir model with a complicated fracture network. FIG. 12B presents an unstructured mesh for the reservoir model. FIGS. 13A-13C present the properties, respectively, fracture permeability, fracture porosity, and fracture aperture, assigned to the fractures.

[0090] This example demonstrates the effectiveness of the hybrid MFD-SL on reservoir models featuring intricate fracture networks. The setup includes two injectors and four producers positioned at the middle and four corners of the domain, respectively (see FIG. 13A). The reservoir model adopts an unstructured mesh, as depicted in FIG. 13B. Detailed data related to the simulation can be found in Table 1, except the production rates for all producers are set to be10⁢ m3day.Furthermore, heterogeneous properties in the fractures, such as permeability, porosity, and aperture, are added with the values shown in FIGS. 14A-14C, respectively. This adds complexity and realism to the model, enabling a more accurate representation of real-world reservoir conditions.FIGS. 15A and 15B illustrate the streamlines generated based on the reservoir model of FIG. 13A. FIGS. 15C-15F present numerical results for comparison. FIGS. 15 and 15B illustrates the streamline distribution and FIGS. 15C-15F compare oil saturation profiles between the hybrid MFD-SL (see FIGS. 15E and 15F) and full MFD based on DFM (see FIGS. 15C and 15D). Results show that the hybrid MFD-SL provides similar solutions as the full MFD on DFM in terms of the shape of the water-flooding area with lower numerical diffusion errors (i.e., narrow bandwidth of flooding front). Furthermore, the implementation of streamline techniques enables larger time steps compared to the full MFD scheme while maintaining accuracy, thus improving the overall efficiency of the simulation process. This advantage makes the hybrid MFD-SL method particularly valuable in effectively analyzing reservoir behaviors and performing quick flow diagnostics.Hydrocarbon Operations

[0092] FIG. 16 illustrates hydrocarbon production operations 1600 that include both one or more field operations 1610 and one or more computational operations 1612, which exchange information and control exploration for the production of hydrocarbons. In some implementations, outputs of techniques of the present disclosure can be performed before, during, or in combination with the hydrocarbon production operations 1600, specifically, for example, either as field operations 1610 or computational operations 1612, or both.

[0093] Examples of field operations 1610 include forming / drilling a wellbore, hydraulic fracturing, producing through the wellbore, injecting fluids (such as water) through the wellbore, to name a few. In some implementations, methods of the present disclosure can trigger or control the field operations 1610. For example, the methods of the present disclosure can generate data from hardware / software including sensors and physical data gathering equipment (e.g., seismic sensors, well logging tools, flow meters, and temperature and pressure sensors). The methods of the present disclosure can include transmitting the data from the hardware / software to the field operations 1610 and responsively triggering the field operations 1610 including, for example, generating plans and signals that provide feedback to and control physical components of the field operations 1610. Alternatively or in addition, the field operations 1610 can trigger the methods of the present disclosure. For example, implementing physical components (including, for example, hardware, such as sensors) deployed in the field operations 1610 can generate plans and signals that can be provided as input or feedback (or both) to the methods of the present disclosure.

[0094] Examples of computational operations 1612 include one or more computer systems 1620 that include one or more processors and computer-readable media (e.g., non-transitory computer-readable media) operatively coupled to the one or more processors to execute computer operations to perform the methods of the present disclosure. The computational operations 1612 can be implemented using one or more databases 1618, which store data received from the field operations 1610 and / or generated internally within the computational operations 1612 (e.g., by implementing the methods of the present disclosure) or both. For example, the one or more computer systems 1620 process inputs from the field operations 1610 to assess conditions in the physical world, the outputs of which are stored in the databases 1618. For example, seismic sensors of the field operations 1610 can be used to perform a seismic survey to map subterranean features, such as facies and faults. In performing a seismic survey, seismic sources (e.g., seismic vibrators or explosions) generate seismic waves that propagate in the earth and seismic receivers (e.g., geophones) measure reflections generated as the seismic waves interact with boundaries between layers of a subsurface formation. The source and received signals are provided to the computational operations 1612 where they are stored in the databases 1618 and analyzed by the one or more computer systems 1620.

[0095] In some implementations, one or more outputs1622 generated by the one or more computer systems 1620 can be provided as feedback / input to the field operations 1610 (either as direct input or stored in the databases 1618). The field operations 1610 can use the feedback / input to control physical components used to perform the field operations 1610 in the real world.

[0096] For example, the computational operations 1612 can process the seismic data to generate three-dimensional (3D) maps of the subsurface formation. The computational operations 1612 can use these 3D maps to provide plans for locating and drilling exploratory wells. In some operations, the exploratory wells are drilled using logging-while-drilling (LWD) techniques which incorporate logging tools into the drill string. LWD techniques can enable the computational operations 1612 to process new information about the formation and control the drilling to adjust to the observed conditions in real-time.

[0097] The one or more computer systems 1620 can update the 3D maps of the subsurface formation as information from one exploration well is received and the computational operations 1612 can adjust the location of the next exploration well based on the updated 3D maps. Similarly, the data received from production operations can be used by the computational operations 1612 to control components of the production operations. For example, production well and pipeline data can be analyzed to predict slugging in pipelines leading to a refinery and the computational operations 1612 can control machine operated valves upstream of the refinery to reduce the likelihood of plant disruptions that run the risk of taking the plant offline.

[0098] In some implementations of the computational operations 1612, customized user interfaces can present intermediate or final results of the above-described processes to a user. Information can be presented in one or more textual, tabular, or graphical formats, such as through a dashboard. The information can be presented at one or more on-site locations (such as at an oil well or other facility), on the Internet (such as on a webpage), on a mobile application (or app), or at a central processing facility.

[0099] The presented information can include feedback, such as changes in parameters or processing inputs, that the user can select to improve a production environment, such as in the exploration, production, and / or testing of petrochemical processes or facilities. For example, the feedback can include parameters that, when selected by the user, can cause a change to, or an improvement in, drilling parameters (including drill bit speed and direction) or overall production of a gas or oil well. The feedback, when implemented by the user, can improve the speed and accuracy of calculations, streamline processes, improve models, and solve problems related to efficiency, performance, safety, reliability, costs, downtime, and the need for human interaction.

[0100] In some implementations, the feedback can be implemented in real-time, such as to provide an immediate or near-immediate change in operations or in a model. The term real-time (or similar terms as understood by one of ordinary skill in the art) means that an action and a response are temporally proximate such that an individual perceives the action and the response occurring substantially simultaneously. For example, the time difference for a response to display (or for an initiation of a display) of data following the individual's action to access the data can be less than 1 millisecond (ms), less than 1 second(s), or less than 5 s. While the requested data need not be displayed (or initiated for display) instantaneously, it is displayed (or initiated for display) without any intentional delay, taking into account processing limitations of a described computing system and time required to, for example, gather, accurately measure, analyze, process, store, or transmit the data.

[0101] Events can include readings or measurements captured by downhole equipment such as sensors, pumps, bottom hole assemblies, or other equipment. The readings or measurements can be analyzed at the surface, such as by using applications that can include modeling applications and machine learning. The analysis can be used to generate changes to settings of downhole equipment, such as drilling equipment. In some implementations, values of parameters or other variables that are determined can be used automatically (such as through using rules) to implement changes in oil or gas well exploration, production / drilling, or testing. For example, outputs of the present disclosure can be used as inputs to other equipment and / or systems at a facility. This can be especially useful for systems or various pieces of equipment that are located several meters or several miles apart, or are located in different countries or other jurisdictions.In Conclusion

[0102] The hybrid MFD-SL is applicable to efficient modeling of two-phase flow in fractured reservoirs. The hybrid approach is implemented within the discrete fracture model (DFM) framework.

[0103] The hybrid MFD-SL approach can achieve more accurate results compared to full MFD schemes when applied to DFM and EDFM with lower computation costs. These advancements are contributed by the use of streamline, which significantly alleviates the numerical diffusion error and improves the computation efficiency. The proposed approach provides wider applicability than approaches based on EDFM due to the inherent limitation of EDFM in handling complex cases, such as multi-phase flow across flow barriers and highly conductive fractures. The proposed approach provides a practical and easy-to-implement streamline tracing method developed on triangular and tetrahedral grids commonly implemented in DFM.EXAMPLES

[0104] In some implementations, methods include: receiving a discrete fracture model representing the fractured subsurface formation, the discrete fracture model comprising matrix cells and fracture cells; receiving pressure values and saturation values for multiple fluid phases across the discrete fracture model; based on the pressure values for the multiple fluid phases across the discrete fracture model, generating face-centroid velocities and for the matrix cells and the fracture cells and updating the pressure values for the multiple fluid phases across the discrete fracture model by performing operations including a mimetic finite difference analysis; based on the generated face-centroid velocities, determining an exit face and time-of-flight for each cell and updating the saturation values for the multiple fluid phases across the discrete fracture model based on the exit face and time-of-flight for each cell; storing the pressure values and the saturation values for discrete fracture model and incrementing a time counter; iteratively repeating the generating, determining, and storing steps until the time counter reaches a set value; and based on stored the pressure values and the saturation values, controlling pumping from a production well extending into the fractured subsurface formation.

[0105] In some implementations, methods include: receiving, by a processor, a discrete fracture model comprising matrix cells and fracture cells representing the fractured subsurface formation; receiving, by the processor, pressure values and saturation values for multiple fluid phases across the discrete fracture model; based on the pressure values for the multiple fluid phases across the discrete fracture model, generating, by the processor, face-centroid velocities and for the matrix cells and the fracture cells and updating the pressure values for the multiple fluid phases across the discrete fracture model by performing operations including a mimetic finite difference analysis; based on the generated face-centroid velocities, determining, by the processor, an exit face and time-of-flight for each cell and updating the saturation values for the multiple fluid phases across the discrete fracture model based on the exit face and time-of-flight for each cell; storing the pressure values and the saturation values for discrete fracture model and incrementing a time counter; and iteratively repeating the generating, determining, and storing steps until the time counter reaches a set value.

[0106] In an example implementation combinable with any other example implementation, methods also include measuring the pressure values and the saturation values.

[0107] In an example implementation combinable with any other example implementation, methods also include generating data from the fractured subsurface formation and generating the discrete fracture model.

[0108] In an example implementation combinable with any other example implementation, the multiple fluid phases comprise oil and water. In some cases, the multiple fluid phases consist of oil and water.

[0109] In an example implementation combinable with any other example implementation, receiving pressure values and saturation values comprises

[0110] In an example implementation combinable with any other example implementation, receiving the discrete fracture model comprises receiving a previously generated discrete fracture model.

[0111] In an example implementation combinable with any other example implementation, the saturation values are calculated using a streamline-based method. In some cases, the streamline-based method comprises evaluating, for each face of a matrix cell, whether streamlines interest the face being evaluated. In some cases, the methods also include identifying which face of the matrix cell being evaluated the streamlines exit through.

[0112] A number of embodiments of the systems and methods have been described. Nevertheless, it will be understood that various modifications may be made without departing from the spirit and scope of this specification. Accordingly, other embodiments are within the scope of the following claims.

Claims

1. A method for performing hydrocarbon extraction from a fractured subsurface formation, the method comprising:receiving a discrete fracture model representing the fractured subsurface formation, the discrete fracture model comprising matrix cells and fracture cells;receiving pressure values and saturation values for multiple fluid phases across the discrete fracture model;based on the pressure values for the multiple fluid phases across the discrete fracture model, generating face-centroid velocities and for the matrix cells and the fracture cells and updating the pressure values for the multiple fluid phases across the discrete fracture model by performing operations including a mimetic finite difference analysis;based on the generated face-centroid velocities, determining an exit face and time-of-flight for each cell and updating the saturation values for the multiple fluid phases across the discrete fracture model based on the exit face and time-of-flight for each cell;storing the pressure values and the saturation values for discrete fracture model and incrementing a time counter;iteratively repeating the generating, determining, and storing steps until the time counter reaches a set value; andbased on stored the pressure values and the saturation values, controlling pumping from a production well extending into the fractured subsurface formation.

2. The method of claim 1, further comprising measuring the pressure values and the saturation values.

3. The method of claim 1, further comprising generating data from the fractured subsurface formation and generating the discrete fracture model.

4. The method of claim 1, wherein the multiple fluid phases comprise oil and water.

5. The method of claim 4, wherein the multiple fluid phases consist of oil and water.

6. The method of claim 1, wherein receiving pressure values and saturation values comprises receiving previously measured pressure values from a database.

7. The method of claim 1, wherein receiving the discrete fracture model comprises receiving a previously generated discrete fracture model.

8. The method of claim 1, wherein the saturation values are calculated using a streamline-based method.

9. The method of claim 8, wherein the streamline-based method comprises evaluating, for each face of a matrix cell, whether streamlines interest the face being evaluated.

10. The method of claim 9, further comprising identifying which face of the matrix cell being evaluated the streamlines exit through.

11. A method for performing hydrocarbon extraction from a fractured subsurface formation, the method comprising:receiving, by a processor, a discrete fracture model comprising matrix cells and fracture cells representing the fractured subsurface formation;receiving, by the processor, pressure values and saturation values for multiple fluid phases across the discrete fracture model;based on the pressure values for the multiple fluid phases across the discrete fracture model, generating, by the processor, face-centroid velocities and for the matrix cells and the fracture cells and updating the pressure values for the multiple fluid phases across the discrete fracture model by performing operations including a mimetic finite difference analysis;based on the generated face-centroid velocities, determining, by the processor, an exit face and time-of-flight for each cell and updating the saturation values for the multiple fluid phases across the discrete fracture model based on the exit face and time-of-flight for each cell;storing the pressure values and the saturation values for discrete fracture model and incrementing a time counter; anditeratively repeating the generating, determining, and storing steps until the time counter reaches a set value.

12. The method of claim 11, wherein the multiple fluid phases comprise oil and water.

13. The method of claim 12, wherein the multiple fluid phases consist of oil and water.

14. The method of claim 12, wherein receiving pressure values and saturation values comprises receiving previously measured pressure values from a database.

15. The method of claim 14, wherein receiving the discrete fracture model comprises receiving a previously generated discrete fracture model.

16. The method of claim 15, wherein the saturation values are calculated using a streamline-based method.

17. The method of claim 16, wherein the streamline-based method comprises evaluating, for each face of a matrix cell, whether streamlines interest the face being evaluated.

18. The method of claim 17, further comprising identifying which face of the matrix cell being evaluated the streamlines exit through.

Images