If secondary hydrocarbon recovery methods, like water flooding, fail because of the occurrence of viscous fingering one can turn to an enhanced oil recovery method (EOR) like the injection of foam. The generation of foam in a porousmedium can be described by a set of partial differential equations with strongly non-linear functions, which impose challenges for the numerical modeling. Former studies [1–3] show the occurrence of strongly temporally oscillating solutions when using forward simulation models, that are entirely due to discretization artifacts. We describe the foam process by an immiscible two-phase flow model where gas is injected in a porousmedium filled with a mixture of water and surfactants. The change from pure gas into foam is incorporated in the model through a reduction in the gas mobility. Hence, the two-phase description of the flow stays intact. Since the total pressure drop in the reservoir is small, both fluids can be considered incompressible [3]. However, whereas the fractional flow function for a gas-flooding process is a smooth function of water saturation, the generation of foam will cause a rapid increase of the flux function over a very small saturation scale. Consequently, the derivatives of the flux function can become extremely large and impose a severe constraint on the time step. We address the stability issues of the foam model, by numerous numerical approaches that improve the accuracy of the solutions. First, we study several averaging schemes and introduce a novel way of approximating the foam mobility functions on the grid interfaces in a finite volume framework. This will lead to solutions that are significantly smoother than can be achieved with standard averaging schemes. Next, we discuss several novel discretization schemes where the discontinuity is incorporated in the numerical fluxes for a simplified compressible flow model. These include the indirect addition of an extra grid interface at the location of the discontinuity, to preserve monotonicity of the solutions in time. Variations on this method, are the addition of an extra grid cell around the highly non-linear phase transition and the adaption of the flux terms based on the location of the discontinuity or non-linearity in the grid. As a practical example to demonstrate these techniques we study a simplified model for foam flow in porous media. The model is then extended to a two-dimensional reservoir, where the accuracy of the solutions is a main concern. The two-dimensional simulator that is used for this, was build and tested for the foam model. It includes higher-order hyperbolic Riemann solvers, and flux correction schemes to compute the saturation of the different fluid phases in the model. The elliptic solver for the pressure equation is also adapted to the stiffness of the problem. With this simulator we perform a quantitative study of the stability characteristics of the flow, to gain more insight in the important wave-lengths and scales of the foam model. This insight forms an essential step towards the design of a suitable computational solver that captures all the appropriate scales, while retaining computational efficiency. In addition, we present a qualitative analysis of the effect of different reservoir and fluid properties on the foam fingering behavior. In particular, we consider the effect of heterogeneity of the reservoir, injection rates, and foam quality. This leads to interesting observations about the influence of the different foam parameters on the stability of the solutions, and we are able to predict the flow stability for different foam qualities. Finally, we discuss several other approaches that were addressed during this PhD-project to increase the understanding of solving highly non-linear flow problems in a porous medium. @en

Accurate field-scale simulations of foam enhanced oil recovery are challenging, due to the sharp transition between gas and foam. Hence, unpredictable numerical and physical behavior is often observed, casting doubt on the validity of the simulation results. In this paper, a thorough stability analysis of the foam model is presented to validate the simulation results. We study the effect of a strongly non-monotonous total mobility function arising from foam models on the stability characteristics of the flow. To this end, we apply the linear stability analysis to nearly discontinuous relative permeability functions and compare the results with those of highly accurate numerical simulations. In addition, we present a qualitative analysis of the effect of different reservoir and fluid properties on the foam fingering behavior. In particular, we consider the effect of heterogeneity of the reservoir, injection rates, and foam quality. Relative permeability functions play an important role in the onset of fingering behavior of the injected fluid. Hence, we can deduce that stability properties are highly dependent on the non-linearity of the foam transition. The foam-water interface is governed by a very small total mobility ratio, implying a stable front. The transition between gas and foam, however, exhibits a huge total mobility ratio, leading to instabilities in the form of viscous fingering. This implies that there is an unstable pattern behind the front. We deduce that instabilities are able to grow behind the front but are later absorbed by the expanding wave. Moreover, the stability analysis, validated by numerical simulations, provides valuable insights about the important scales and wavelengths of the foam model. In this way, we remove the ambiguity regarding the effect of grid resolution on the convergence of the solutions. This insight forms an essential step toward the design of a suitable computational solver that captures all the appropriate scales, while retaining computational efficiency.@en

Many enhanced oil recovery (EOR) processes can be described using partial differential equations with parameters that are strongly non-linear functions of one or more of the state variables. Typically these nonlinearities result in solution components changing several orders of magnitude over small spatial or temporal distances. The numerical simulation of such processes with the aid of finite volume or finite element techniques poses challenges. In particular, temporally oscillating state variable values are observed for realistic grid sizes when conventional discretization schemes are used. These oscillations, which do not represent a physical process but are discretization artifacts, hamper the use of the forward simulation model for optimization purposes. To analyze these problems, we study the dynamics of a simple foam model describing the interaction of water, gas and surfactants in a porous medium. It contains sharp gradients due to the formation of foam. The simplicity of the model allows us to gain a better understanding of the underlying processes and difficulties of the problem. The foam equations are discretized by a first-order finite volume method. Instead of using a finite volume method with a standard interpolation procedure, we opt for an integral average, which smooths out the discontinuity caused by foam generation. We introduce this method by applying it to the heat equation with discontinuous thermal conductivity. A similar technique is then applied to the foam model, reducing the oscillations drastically, but not removing them. @en