A conservative sequential fully implicit method is derived for compositional reservoir simulation. Multi-phase flow in porous media comprises coupled complex processes: i.e. elliptic flow equation, hyperbolic transport equation and highly nonlinear phase equilibrium equation. These processes contain very different mathematical characteristics that cannot be efficiently solved by one numerical method. As a result, the fully implicit method may become numerically complex and inefficient because the Jacobian includes the derivatives w.r.t. the variables from all of the different processes involved. Jenny et al. (2004) [12] and Lee et al. (2015) [20] demonstrated that flow (pressure) and transport (saturation) for multi-phase flow without compositional effect can be efficiently solved by a sequential fully implicit method. However, the characteristics of the phase equilibrium equations are very different from those of the transport equations. This paper proposes an iterative method that solves the flow, transport and phase equilibrium equations in a sequential manner. The transport of hydrocarbons through porous media is governed by the multi-phase Darcy's equation, which is used to compute the phase velocities. The hydrocarbon components belonging to the same phase are transported with the same phase velocity. Upon arrival in the destination grid cell, these components are redistributed via a phase equilibrium calculation. This observation leads to simplification of the governing equations by reducing primary variables to four (i.e., pressure and three phase saturations). The nonlinear solution scheme composed of the stages outlined above is proven to preserve mass conservation, while a new degree of freedom, “thermodynamic flux”, is introduced to ensure volume conservation. The sequential algorithm is solved iteratively until pressure, saturation, and phase composition are fully converged. It is well-known that sequential solution schemes may require many iterations or fail to converge if the phase equilibrium calculation involves phase transition with a large volume change. This indicates that the current governing equations may not adequately describe fluid flux during rapid phase transition. With numerical examples we demonstrate that such numerical difficulties are successfully resolved via the thermodynamic flux term.

@en

Despite welcome increases in the adoption of renewable energy sources, oil and natural gas are likely to remain the main ingredient in the global energy diet for the decades to come. Therefore, the efficient exploitation of existing suburface reserves is essential for the well-being of society. This has stimulated recent developments in computer models able to provide critical insight into the evolution of the flow of water, gas and hydrocarbons through rock pores. Any such endeavour, however, has to tackle a number of challenges, including the considerable size of the domain, the highly heterogeneous spatial distribution of geological properties, as well as the intrinsic uncertainty and limitations associated with field data acquisition. In addition, the naturally-formed or artificially induced networks of fractures, present in the rock, require special treatment, due to their complex geometry and crucial impact on fluid flow patterns. From a numerical point of view, a reservoir simulator’s operation entails the solution of a series linear systems, as dictated by the spatial and temporal discretization of the governing equations. The difficulty lies in the properties of these systems, which are large, ill-conditioned and often have an irregular sparsity pattern. Therefore, a brute-force approach, where the solutions are directly computed at the original fine-scale resolution, is often an impractically expensive venture, despite recent advances in parallel computing hardware. On the other hand, switching to a coarser resolution to obtain faster results, runs the risk of omitting important features of the flow, which is especially true in the case of fractured porous media. This thesis describes an algebraic multiscale approach for fractured reservoir simulation. Its purpose is to offer a middle-ground, by delivering results at the original resolution, while solving the equations on the coarse-scale. This is made possible by the so-called basis functions – a set of locally-supported cross-scale interpolators, conforming to the heterogeneities in the domain. The novelty of the work lies in the extension of these methods to capture the effect of fractures. Importantly, this is done in fully algebraic fashion, i.e. without making any assumptions regarding geometry or conductivity properties. In order to elicit the generality of the proposed approach, a series of sensitivity studies are conducted on a proof-of-concept implementation. The results, which include both CPU times and convergence behaviour, are discussed and compared to those obtained using an industrial-grade AMG package. They serve as benchmarks, recommending the inclusion of multiscale methods in next-generation commercial reservoir simulators.@en

This paper describes a novel sensitivity analysis method, able to handle dependency relationships between model parameters. The starting point is the popular Morris (1991) algorithm, which was initially devised under the assumption of parameter independence. This important limitation is tackled by allowing the user to incorporate dependency information through a copula. The set of model runs obtained using latin hypercube sampling, are then used for deriving appropriate sensitivity measures. Delft3D-WAQ (Deltares, 2010) is a sediment transport model with strong correlations between input parameters. Despite this, the parameter ranking obtained with the newly proposed method is in accordance with the knowledge obtained from expert judgment. However, under the same conditions, the classic Morris method elicits its results from model runs which break the assumptions of the underlying physical processes. This leads to the conclusion that the proposed extension is superior to the classic Morris algorithm and can accommodate a wide range of use cases.@en

This work presents a new discrete fracture model, namely the Projection-based Embedded Discrete Fracture Model (pEDFM). Similar to the existing EDFM approach, pEDFM constructs independent grids for the matrix and fracture domains, and delivers strictly conservative velocity fields. However, as a significant step forward, it is able to accurately model the effect of fractures with general conductivity contrasts relative to the matrix, including impermeable flow barriers. This is achieved by automatically adjusting the matrix transmissibility field, in accordance to the conductivity of neighboring fracture networks, alongside the introduction of additional matrix-fracture connections. The performance of pEDFM is investigated extensively for two- and three-dimensional scenarios involving single-phase as well as multiphase flows. These numerical experiments are targeted at determining the sensitivity of the model towards the fracture position within the matrix control volume, grid resolution and the conductivity contrast towards the matrix. The pEDFM significantly outperforms the original EDFM and produces results comparable to those obtained when using DFM on unstructured grids, therefore proving to be a flexible model for field-scale simulation of flow in naturally fractured reservoirs.@en

A novel multiscale method for discrete fracture modeling on unstructured grids (MS-DFM) is developed. To this end, the DFM fine-scale discrete system is constructed using unstructured conforming cells for the matrix with lower-dimensional fracture elements placed at their interfaces. On this unstructured fine grid, MS-DFM imposes independent unstructured coarse grids for the fracture and matrix domains. While the conservative coarse-scale system is solved over these coarse-grid cells, overlapping dual-coarse blocks are also formed in order to provide local supports for the multiscale basis functions. To increase the accuracy, but maintaining the computational efficiency, fracture-matrix coupling is considered only for the basis functions inside the matrix domain. This results in additional (enriching) fracture basis functions in the matrix. By construction, basis functions form the partition of unity for both fracture and matrix sub-domains. Furthermore, to enable error reduction to any desired level, a convergent iterative strategy is developed, where MS-DFM is employed along with a fine-scale smoother in order to resolve low- and high-frequency modes in the error. The performance of MS-DFM is assessed for several 2D and 3D test cases. The proposed method achieves accurate results for several test cases even without iterations, and for challenging ones with only a few iterations. MS-DFM is the first of its kind, and thus extends the application of multiscale methods to unstructured discrete fracture models. As such, it provides a promising framework for real-field application of unstructured DFM.@en

A multiscale method for Discrete Fracture Modeling (DFM) using unstructured grids is developed. The fine-scale discrete system is obtained by imposing tetrahedron (triangular for 2D domains) shaped grid cells, while lower-dimensional fractures are imposed at the grid interfaces. The DFM approach is then used to describe the transmissibility coefficients for all the interfaces, including those with the lower-dimensional fractures. On this fine-scale discrete system, a new algebraic multiscale formulation is developed, which first imposes two sets of coarseand dual-coarse grids. The former grid is essential for conservative multiscale formulations, while the second one is used for the calculation of local multiscale basis functions. The coarse-scale partitioning of the fracture and matrix domains is flexible and totally independent. Moreover, the multiscale basis functions are constructed for both the matrix and fracture domains. By construction, the basis functions are a partition of unity. For 2D and 3D test cases, the performance of the multiscale method is systematically assessed. It is shown that the method (with no multiscale iterations) provides accurate results, even for complex fractured systems. The presented multiscale method is a promising framework for real-field application of DFM models.

@en

A multiscale method for Discrete Fracture Modeling (DFM) using unstructured grids is developed. The fine-scale discrete system is obtained by imposing tetrahedron (triangular for 2D domains) shaped grid cells, while lower-dimensional fractures are imposed at the grid interfaces. The DFM approach is then used to describe the transmissibility coefficients for all the interfaces, including those with the lower-dimensional fractures. On this fine-scale discrete system, a new algebraic multiscale formulation is developed, which first imposes two sets of coarseand dual-coarse grids. The former grid is essential for conservative multiscale formulations, while the second one is used for the calculation of local multiscale basis functions. The coarse-scale partitioning of the fracture and matrix domains is flexible and totally independent. Moreover, the multiscale basis functions are constructed for both the matrix and fracture domains. By construction, the basis functions are a partition of unity. For 2D and 3D test cases, the performance of the multiscale method is systematically assessed. It is shown that the method (with no multiscale iterations) provides accurate results, even for complex fractured systems. The presented multiscale method is a promising framework for real-field application of DFM models.@en

Any endeavour to accurately model flow through fractured porous media at the field-scale must overcome two important challenges. First, the discretized representation of the medium needs to accommodate the complex geometry of intersecting small-scale fractures with various lengths, apertures and orientations. Second, the model formulation must ensure that the conductivity of these fractures, which can be orders of magnitude higher or lower than that of the host rock, is properly taken into account when computing the pressure map. The Embedded Discrete Fracture Model (EDFM) is well known in the literature for its flexibility in representing fractures. More specifically, in EDFM, fractures are lower-dimensional features, discretized independently from the matrix. Their effect on the flow is captured by defining fluxes between the fracture control volumes and the matrix grid cells they intersect. EDFM was proven effective in capturing the flow behaviour through porous media containing highly conductive fractures. However, its formulation fails to represent the effect of low-permeable features, such as embedded flow barriers. In this work, a novel projection-based Embedded Discrete Fracture Model (pEDFM) is introduced for flow simulation in fractured porous media with general conductivity contrasts. Similar to EDFM, pEDFM constructs independent grids for the fracture and matrix domains. As an additional step, the transmissibilities at matrix interfaces are automatically adjusted to account for the conductivity of neighbouring fracture networks, via a scaling factor proportional to their geometric projections. The performance of pEDFM is investigated extensively for two- and three-dimensional scenarios involving single- as well as multiphase flows. These numerical experiments are targeted at determining the sensitivity of the model towards the grid resolution, fracture position and orientation, as well as the conductivity contrast towards the matrix. The results of these studies support the conclusion that pEDFM significantly outperforms the original EDFM model and is a viable method for field-scale simulation of flow in naturally fractured reservoirs.@en

This paper introduces an Algebraic MultiScale method for simulation of flow in heteroge-neous porous media with embedded discrete Fractures (F-AMS). First, multiscale coarse grids are independently constructed for both porous matrix and fracture networks. Then, amap between coarse-and fine-scale is obtained by algebraically computing basis functions with local support. In order to extend the localization assumption to the fractured media, four types of basis functions are investigated: (1)Decoupled-AMS, in which the two media are completely decoupled, (2)Frac-AMS and (3)Rock-AMS, which take into account only one-way transmissibilities, and (4)Coupled-AMS, in which the matrix and fracture interpolators are fully coupled. In order to ensure scalability, the F-AMS framework permits full flexibility in terms of the resolution of the fracture coarse grids. Numerical results are presented for two-and three-dimensional heterogeneous test cases. During these experiments, the performance of F-AMS, paired with ILU(0) as second-stage smoother in a convergent iterative procedure, is studied by monitoring CPU times and convergence rates. Finally, in order to investigate the scalability of the method, an extensive benchmark study is conducted, where a commercial algebraic multigrid solver is used as reference. The results show that, given an appropriate coarsening strategy, F-AMS is insensitive to severe fracture and matrix conductivity contrasts, as well as the length of the fracture networks. Its unique feature is that a fine-scale mass conservative flux field can be reconstructed after any iteration, providing efficient approximate solutions in time-dependent simulations.@en

Mathematical formulations describing flow in porous media typically entail highly heterogeneous coefficients, changing over several orders of magnitude through the entirety of the domain. In addition, many of the target geological formations are fractured. Fractures are lower dimensional manifolds with properties that differ greatly from those of the surrounding porous rock. Therefore, given their significant role in establishing the patterns of the flow regime, accurate representation of fractures within flow and transport models is crucial for many geoscientific applications, including groundwater flow and geothermal energy exploitations. Embedded Discrete Fracture Model (EDFM) employs independent grids for matrix and fractures. This results in efficient computations, specially for complex fracture geometries, and cases with dynamic fracture creations (and closures). Even though small-scale fractures are homogenized within the matrix rock, the remaining explicit fractures (bigger than fine-scale grid resolution) along with heterogeneous matrix, for realistic cases, lead to linear systems which are beyond the scope of classical simulation methods. Multiscale finite element and volume (MSFE and MSFV, respectively) methods have been developed mainly for heterogeneous, but non-fractured, porous media. In order to extend them to account for flow in heterogeneous fractured formations, here, F-AMS is introduced as a novel Algebraic Multiscale Solver. It operates by defining coarse grids for both the porous matrix and the embedded discrete fractures. Then, by computing local basis functions, a general map between the fine- and coarse-scale systems, i.e. the prolongation operator, is obtained. These basis functions form a partition of unity and, in their present formulation, they allow for four degrees of fracture-matrix coupling: (1) Decoupled-AMS, in which the two media are completely decoupled, (2) Frac-AMS allows one-way coupling, where the fracture coarse solutions also affect the matrix fine-scale pressure, (3) Rock-AMS is the counterpart of Frac-AMS, where the matrix coarse solution is also employed to find the fracture fine-scale pressure, and (4) Coupled-AMS, in which matrix and fracture interpolators are fully coupled. If only one coarse degree of freedom (DOF) is considered for each fracture network, the Frac-AMS strategy becomes equivalent to the earlier method proposed by Hajibeygi et al. However, in order to maintain efficiency for general cases, the F-AMS framework permits full flexibility in terms of the definition of the fracture coarse grids and the level of matrix-fracture prolongation coupling. Moreover, by using the Finite Volume restriction operator after any iteration, a mass conservative velocity can be reconstructed and be used to solve the transport equations. Systematic numerical experiments for 3D heterogeneous fractured domains (from $10^5$ to $10^7$ grid cells, and fracture-matrix transmissibility contrasts of $10^1$ to $10^8$) are presented and discussed. In addition, the F-AMS is benchmarked agains SAMG, a commercial Algebraic Multigrid solver. These results illustrate that F-AMS is an efficient multiscale procedure for large-scale fractured reservoirs. It is important to note that for multiphase flow scenarios, only a few F-AMS iterations are sufficient to obtain good quality pressure solutions. These lead to the conclusion that F-AMS is an important multiscale development for the efficient simulation of flow in naturally fractured porous media. @en

A novel multiscale method for multiphase flow in heterogeneous fractured porous media is devised. The discrete fine-scale system is described using an embedded fracture modeling approach, in which the heterogeneous rock (matrix) and highly-conductive fractures are represented on independent grids. Given this fine-scale discrete system, the method first partitions the fine-scale volumetric grid representing the matrix and the lower-dimensional grids representing fractures into independent coarse grids. Then, basis functions for matrix and fractures are constructed by restricted smoothing, which gives a flexible and robust treatment of complex geometrical features and heterogeneous coefficients. From the basis functions one constructs a prolongation operator that maps between the coarse- and fine-scale systems. The resulting method allows for general coupling of matrix and fracture basis functions, giving efficient treatment of a large variety of fracture conductivities. In addition, basis functions can be adaptively updated using efficient global smoothing strategies to account for multiphase flow effects. The method is conservative and because it is described and implemented in algebraic form, it is straightforward to employ it to both rectilinear and unstructured grids. Through a series of challenging test cases for single and multiphase flow, in which synthetic and realistic fracture maps are combined with heterogeneous petrophysical matrix properties, we validate the method and conclude that it is an efficient and accurate approach for simulating flow in complex, large-scale, fractured media.@en

This work presents the formulation of a novel Projection-based Embedded Discrete Fracture Model (pEDFM), and its integration into an algebraic multiscale procedure. Similar to EDFM, pEDFM constructs independent grids for the matrix and fracture domains. However, as a significant step forward, it is able to accurately model the effect of fractures with general conductivity contrasts relative to the matrix, including impermeable flow barriers. This is achieved by automatically adjusting the matrix transmissibility field, in accordance to the conductivity of neighboring fracture networks. Then, in order to extend the pEDFM to real-field applications, F-AMS-pEDFM is introduced, which is an extension of the recently developed algebraic multiscale solver, F-AMS [Ţene et al., 2016], to include pEDFM. The performance (efficiency and scalability) of F-AMS-pEDFM is investigated extensively for challenging two- and three-dimensional scenarios with complex fracture geometries and a wide range of conductivity contrasts. Moreover, FAMS- pEDFM is benchmarked against the commercial SAMG solver, where CPU time is monitored during both the setup and solution phases. The results support the conclusions that (1) pEDFM significantly outperforms the original EDFM model, and (2) the F-AMS-pEDFM approach proposed in this work is an accurate and efficient method for field-scale simulation of flow in fractured reservoirs. @en

This work investigates the possibility of using an incompletely converged pressure solution in the sequentially implicit loop, traditionally used to simulate multiphase flow through porous media. This is made possible through the Multiscale Finite Volume (MSFV) formulation, which was recently extended to include fully integrated basis functions which account for contributions from wells and fractures [Tene et al, SPE RSS 2015]. Here, we take this method one step further towards multiphase flows, by adding a conservative flux reconstruction stage and pairing it with a fine-scale transport solver. In our experiments we monitor the convergence behaviour while adaptively varying the tolerance level for the pressure solution. Finally, we conclude that only few iterations of the multiscale pressure solver are sufficient to obtain results which are meaningful for decision-making, especially given the fact that, in practice, the reservoir specification is subject to geological uncertainties.@en