Search

The governing equations for multiphase flow in porous media have a mixed character, with both nearly elliptic and nearly hyperbolic variables. The flux for each phase can be decomposed into two parts: (1) a geometry- and rock-dependent term that resembles a single-phase flux; and (2) a mobility term representing fluid properties and rock–fluid interactions. The first term is commonly discretized by two- or multipoint flux approximations (TPFA and MPFA, respectively). The mobility is usually treated with single-point upstream weighting (SPU), also known as dimensional or donor cell upstream weighting. It is well known that when simulating processes with adverse mobility ratios, SPU suffers from grid orientation effects. An important example of this, which will be considered in this work, is the displacement of a heavy oil by water. For these adverse mobility ratio flows, the governing equations are unstable at the modeling scale, rendering a challenging numerical problem. These challenges must be addressed in order to avoid systematic biasing of simulation results. In this work, we present a framework for multidimensional upstream weighting for multiphase flow with buoyancy on general two-dimensional grids. The methodology is based on a dual grid, and the resulting transport methods are provably monotone. The multidimensional transport methods are coupled with MPFA methods to solve the pressure equation. Both explicit and fully implicit approaches are considered for time integration of the transport equations. The results show considerable reduction of grid orientation effects compared to SPU, and the explicit multidimensional approach allows larger time steps. For the implicit method, the total number of non-linear iterations is also reduced when multidimensional upstream weighting is used.

Flow in a porous medium can be described by a set of non-linear partial differential equations. The pressure variable satisfies a maximum principle, which guarantees that the solution will have no oscillations. A discretisation of the pressure equation should preserve this monotonicity property. Whether a numerical method is monotone will depend both on the medium and on the grid. We study monotonicity of Multi-point Flux Approximation methods on triangular grids. We derive necessary conditions for monotonicity on uniform grids. Further, we study the robustness of the methods on rough grids, and quantify the violations of the maximum principle. These investigations are done for single phase flow, however, they are supported by two-phase simulations.

Control-volume multipoint flux approximations (MPFA) are discussed for the simulation of complex near-well flow using geometrically flexible grids. Due to the strong non-linearity of the near-well flow, a linear model will, in general, be inefficient. Instead, a model accounting for the logarithmic pressure behavior in the well vicinity is advocated. This involves a non-uniform refinement of the grid in the radial direction. The model accounts for both near-well anisotropies and heterogeneities. For a full simulation involving multiple wells, this single-well approach can easily be coupled with the reservoir model. Numerical simulations demonstrate the convergence behavior of this model using various MPFA schemes under different near-well conditions for single-phase flow regimes. Two-phase simulations support the results of the single-phase simulations.

Numerical ocean modelling is computationally very demanding. Traditionally, the hydrostatic approximation has been applied to reduce the computational burden. This approximation is valid in large scale studies with coarse grid resolution. With faster computers and gradually smaller grid sizes, we may expect that more studies will be performed with non-hydrostatic ocean models. In recent papers several methods for including non-hydrostatic pressure in σ -coordinate models have been suggested. In this paper the sensitivity of the non-hydrostatic pressure field, the velocity fields, and the density fields to changes in the method for computing non-hydrostatic pressure in σ -coordinate ocean models is addressed. The first test case used involves the propagation and breaking of an internal wave at an incline in a tank. The other test case concerns tidally driven flow over a sill in a stratified fjord. The results from our numerical exercises suggest that the velocity and density fields are very robust to the model choices investigated here. The differences between the model results are of the same order as the uncertainty due to the internal pressure gradient error, and they are smaller than an estimate of the uncertainty due to subgrid scale closure.

To address the engineering questions on security issues of geological carbon sequestration (GCS), a broad range of computational models with different levels of complexity have been developed - from multiphase multicomponent fully coupled three-dimensional models to simplified analytical solutions. Within this wide range of models, a family of so-called vertical equilibrium (VE) models has been developed. These models assume that CO2 and brine have segregated due to buoyancy and reached a hydrostatic pressure distribution in the vertical direction. Such VE models are computationally efficient due to the dimensional reduction in the vertical, and accurate as long as the VE assumption is satisfied. However, a study comparing results from a VE model to results from a full three- dimensional model found that there are realistic conditions for which the VE assumption is not justified, especially for geological formations that have low vertical permeability, on the order of 10 milliDarcy or lower [1] . In an attempt to overcome the VE limitation of the vertically-integrated models, a new type of vertically-integrated model which relaxes the VE assumption while still using a vertically-integrated framework has been developed [2] . The new vertically-integrated model is cast into a multiscale framework, where the coarse scale is the horizontal domain and the fine scale is the vertical domain corresponding to the thickness of a geologic formation. This type of model maintains much of the computational advantages of the VE models, while allowing a much wider range of problem to be modelled. In this paper, we extend the model in [2] to include horizontally layered geologic heterogeneities and develop a new dynamic reconstruction model, which we refer to as a “multi-layer dynamic reconstruction” model. The model in [2] is called “single-layer dynamic reconstruction” model to distinguish the modelling approaches. We apply both dynamic reconstruction models to field injection sites, including a hypothetical injection scenario into the Mount Simon formation in the Illinois Basin, USA and the well-known industrial-scale injection into the Utsira formation at the Sleipner site in Norway. The modelling results show that the multi-layer dynamic reconstruction model is capable of dealing with horizontally layered heterogeneities and gives results that agree reasonably well with results from the full multi-dimensional model, although in geologic layers with high permeability the reconstruction algorithm is not able to fully capture the horizontal driving forces due to buoyancy. This could be important over long time scales for highly permeable geologic layers. The single- layer dynamic reconstruction model was shown be to the right model choice for homogeneous formations with relatively low permeability where it takes a long time for CO2 and brine to reach vertical equilibrium [2] . In this study, we found that for homogeneous formations with high permeability and steep capillary pressure curves, the single-layer dynamic reconstruction model gives results that are analogous to vertical equilibrium models, with the fast segregation dynamics requiring small time steps in the dynamic reconstruction algorithm. As such vertical equilibrium models are the appropriate choice for systems with high permeability, although we also note that the behaviour of the brine relative permeability curve at low brine saturations is an additional important consideration and must be included in the overall analysis.

We consider domain decomposition techniques for a non-linear elasticity problem. Our main focus is on non-linear preconditioning, realized in the framework of additive Schwarz preconditioned inexact Newton (ASPIN) methods. The standard 1-level ASPIN method is extended to a 2-level method by adding a non-linear coarse solver. Numerical experiments show that the coarse component is necessary for scalability in terms of linear iterations inside the Newton loop. Moreover, for problems that are dominated by nonlinearities that are not localized in space the non-linear coarse iterations are crucial for achieving computational efficiency. publishedVersion

The governing equations for flow and transport in porous media are derived assuming conservation of mass. To ensure stability of the simulations significant attention is given to ensure that the discrete system retains the conservation property. Due to discretization errors and parameter uncertainty it is natural to consider an inexact solution strategy for the resulting system of equations. However most linear solvers are not designed by the same principles as the underlying discretization and will thus not produce inexact solutions that preserve the conservation property. In this work we illustrate how inexact yet conservative linear solvers can be realized for porous media applications. The linear solver is formulated as a multi-level control volume methods and produces a conservative flux field for all approximated solutions.

The standard approximation for the flow-pressure relationship in porous media is Darcy's law that was originally derived for infiltration of water in fine homogeneous sands. Ever since there have been numerous attempts to generalize it for handling more complex flows. Those include upscaling of standard continuum mechanics flow equations from the fine scale. In this work we present a heterogeneous multiscale method that utilizes fine scale information directly to solve problems for general single phase flow on the Darcy scale. On the coarse scale it only assumes mathematically justified conservation of mass on control volumes, that is, no phenomenological Darcy-type relationship for velocity is presumed. The fluid fluxes are instead provided by a fine scale Navier--Stokes mixed finite element solver. This work also considers several choices of quadrature for data estimation in the multiscale method and compares them. We prove that for an essentially linear regime, when the fine scale is governed by Stokes flow, our method converges to a rigorously derived homogenization solution---Darcy's law. Moreover the method gives the flexibility to solve problems with faster nonlinear flow regimes that is important in a number of applications, such as flows that may occur near wells and in fractured regions in subsurface. Those flows are also common for industrial and near surface porous media. The numerical examples presented in the paper verify the estimate and emphasize the importance of good data estimation. publishedVersion

Control volume methods are prevailing for solving the potential equation arising in porous media flow. The continuous form of this equation is known to satisfy a maximum principle, and it is desirable that the numerical approximation shares this quality. Recently, sufficient criteria were derived guaranteeing a discrete maximum principle for a class of control volume methods. We show that most of these criteria are also necessary. An implication of our work is that no linear nine-point control volume method can be constructed for quadrilateral grids in 2D that is exact for linear solutions while remaining monotone for general problems.

Understanding flow in subsurface porous media is of great importance for society due to applications such as energy extraction and waste disposal. The governing equations for subsurface flow are a set of non-linear partial differential equations of mixed elliptichyperbolic type, and the parameter fields are highly heterogeneous with characteristic features on a continuum of length scales. This calls for robust discretization methods that balance the challenges in designing efficient and accurate methods. In this chapter we focus on a class of linear solvers for elliptic systems that aims at providing fast approximate solutions, preferably in one iteration, but fall back to being iterative methods with good convergence properties if higher accuracy is needed.

When injected in deep saline aquifers, CO2 moves radially away from the injection well and progressively higher in the formation because of buoyancy forces. Analyzes have shown that after the injection period, CO2 will potentially migrate over several kilometers in the horizontal direction but only tens of meters in the vertical direction, limited by the aquifer caprock [1, 2]. Because of the large horizontal plume dimensions, three-dimensional numerical simulations of the plume migration over long periods of time are computationally intensive. Thus, to get results within a reasonable time frame, one is typically forced to use coarse meshes and long time steps which result in inaccurate results because of numerical errors in resolving the plume tip. Given the large aspect ratio between the vertical and horizontal plume dimensions, it is reasonable to approximate the CO2 migration using vertically averaged models. Such models can, in many cases, be more accurate than coarse three-dimensional computations. In particular, models based on vertical equilibrium (VE) [3] are attractive to simulate the long-term fate of CO2 sequestered into deep saline aquifers. The reduced spatial dimensionality resulting from the vertical integration ensures that the computational performance of VE models exceeds the performance of standard three-dimensional models. Thus, VE models are suitable to study the long-time and large-scale behavior of plumes in real large-scale CO2-injection projects [4, 1, 2, 5]. We investigate the use of VE models to simulate CO2 migration in a real large-scale eld case based on data from the Sleipner site in the North Sea. We discuss the potential and limitations of VE models and show how VE models can be used to give reliable estimates of long-term CO2 migration. In particular, we focus on a VE formulation that incorporates the aquifer geometry and heterogeneity, and that considers the eects of hydrodynamic and residual trapping. We compare the results of VE simulations with standard reservoir simulation tools on test cases and discuss their advantages and limitations and show how, provided that certain conditions are met, they can be used to give reliable estimates of long-term CO2 migration.

The geological structure of porous rocks include irregular geometries and heterogeneities on a multiple of scales. Reservoir simulations (e.g. oil recovery, CO2 storage, ground water flow) often involve large spatial scales, where local fine-scale heterogeneities might have an important impact on the global flow. Discretisations of the governing equations of multi-phase flow in general render large systems of non-linear equations. The multiscale nature of the geology is inherited by the discrete mathematical problem, which can be ill-conditioned and thereby hard to solve. In this work, we apply domain decomposition (DD) as a solution strategy for the equations. The splitting of the global system into local subproblems can be done either prior to, or after linearisation of the non-linear system. Comparison between DD as a non-linear and linear solver strategy indicates that splitting into subproblems should be considered for use as preconditioners, also for non-linear systems of equations. When applied after linearisation, DD is best suited as a preconditioner in an iterative solution procedure. In this work, we consider a two-level DD framework for linear systems. The framework is flexible, and can be used both as an upscaling procedure, and as a preconditioner.

The governing equations for multi-phase flow in porous media often have a mixed elliptic and (nearly) hyperbolic character. The total flux for each phase consists of two parts: a geometry and rock dependent term that resembles a single-phase flux and a mobility term representing fluid properties and rock-fluid interactions. The geometric term is commonly discretized by two or multi point flux approximations (TPFA and MPFA, respectively). The mobility is usually treated with single point upstream weighting (SPU), also known as dimensional or donor cell upstream weighting. It is well known that when simulating processes with adverse mobility ratios, e.g. gas injection, SPU yields grid orientation effects. For these physical processes, the governing equations are unstable on the scale at which they are discretized, rendering a challenging numerical problem. These challenges must be addressed in order to avoid systematic biasing of simulation results and to improve the overall performance prediction of enhanced oil recovery processes. In this work, we present a framework for multi-dimensional upstream weighting for multi-phase flow with gravity on general two-dimensional grids. The methodology is based on a dual grid, and the resulting transport methods are provably monotone. The multi-dimensional transport methods are coupled with MPFA methods to solve the pressure equation. Both explicit and fully implicit approaches are considered for treatment of the transport equations. The results show considerable reduction of grid orientation effects compared to SPU, and the explicit multi-dimensional approach allows larger time steps. For the implicit method, the total number of non-linear iterations is also reduced when multi-dimensional upstream weighting is used.

In reservoir simulation, a proper treatment of wells is crucial for the performance of fluid flow models. Important well parameters such as the well flow rate and the wellbore pressure are highly sensitive to the computational accuracy of near-well flow models. Near-well regions are characterised as high flow density regions, and the dominating flow pattern exhibits a radial-like nature with large pressure gradients and, for multiphase flow, large saturation gradients. Skew or horizontal wells will also in general imply a strong effect of anisotropy and heterogeneities due to geological layering of the reservoirs. Due to difference in the reservoir scale and the wellbore radius, reservoir flow models do not capture the true flow behavior in the well vicinity. Thus, more flexible models with local grid refinement in near-well regions are needed. The models should also be adapted to handle complex geological near-well structures and multiphase flow simulations. Existing near-well models are in general based on homogeneous media in the well vicinity. Anisotropies and strong heterogeneities are less accounted for in near-well flow simulations. In this work, we construct analytical solutions for near-well flow which is not aligned with a radial inflow pattern. These solutions resemble strongly heterogeneous, possible anisotropic media. We compare different control volume discretisation schemes and radial-type grids for such cases and give their convergence behavior. The objective is to obtain a clearer view on the accuracy of the near-well grids and discretisation schemes for large contrast in permeability, and hence, to determine which is the preferable grid and a suitable numerical scheme given certain near-well conditions. The simulations are performed for singlephase flow, however, the results will also have relevance for multiphase flow simulations.

In multi-phase flow simulation, the prevailing approach to discretizing flux terms treats the elliptic and the hyperbolic terms in the equations separately. With this concept, the flux is calculated analogously as for single phase flow, and later multiplied by the upstream mobility. This approach is valid when the mobility is a scalar quantity, which is the case for most traditional models. However, tensorial relative permeability (and thus tensorial mobilities) in general arise on all scales, as is seen in both laboratory and field experiments. Furthermore, upscaling methods almost invariably lead to anisotropic relative permeabilities. The saturation dependency of the fluid permeability tensor means that the upstream direction is no longer uniquely defined, which challenges common numerical schemes. In this work, we give examples of how the relative permeability has preferential flow directions on different scales. We then study how to incorporate tensorial relative permeability fields into control volume methods. In particular, we address two immediate challenges: Firstly, the non-determinacy of the upstream direction invalidates the use of upstream weighting for the saturation equation. We discuss the applicability of Godunov methods to handle flow cases where interfaces have either no or two upstream directions. Secondly, there is a marked increase in computational complexity associated with the pressure equation. We discuss the possibility of mitigating this computational complexity through mass lumping methods. The validity and computational efficiency of our approaches is discussed on theoretical terms and with the support of numerical implementations.

In fractured poroelastic media under high differential stress, the shearing of fractures and faults and the corresponding propagation of wing cracks can be induced by fluid injection. Focusing on low-pressure stimulation with fluid pressures below the minimum principal stress but above the threshold required to overcome the fracture's frictional resistance to slip, this paper presents a mathematical model and a numerical solution approach for coupling fluid flow with fracture shearing and propagation. Numerical challenges are related to the strong coupling between hydraulic and mechanical processes, the material discontinuity the fractures represent in the medium, the wide range of spatial scales involved, and the strong effect that fracture deformation and propagation have on the physical processes. The solution approach is based on a multiscale strategy. In the macroscale model, flow in and poroelastic deformation of the matrix are coupled with the flow in the fractures and fracture contact mechanics, allowing fractures to frictionally slide. Fracture propagation is handled at the microscale, where the maximum tangential stress criterion triggers the propagation of fractures, and Paris' law governs the fracture growth processes. Simulations show how the shearing of a fracture due to fluid injection is linked to fracture propagation, including cases with hydraulically and mechanically interacting fractures.

Staggered schemes for compressible flows are highly non linear and the stability analysis has historically been performed with a heuristic approach and the tuning of numerical parameters [12]. We investigate the \(L^2\)-stability of staggered schemes by analysing their numerical diffusion operator. The analysis of the numerical diffusion operator gives new insight into the scheme and is a step towards a proof of linear stability or stability for almost constant initial data. For most classical staggered schemes [9, 10, 11, 14], we are able to prove the positivity of the numerical diffusion only in specific cases (constant sign velocities). We then propose a class of linearly \(L^2\)-stable staggered schemes for the isentropic Euler equations based on a carefully chosen numerical diffusion operator. We give an example of such a scheme and present some first numerical results on a Riemann problem.

International audience; We introduce a three-dimensional Hybrid High-Order method for magnetostatic problems. The proposed method is easy to implement, supports general polyhedral meshes, and allows for arbitrary orders of approximation.