Numerical methods with reduced grid dependency for enhanced oil recovery

In recent years there has been a resurgence of interest in enhanced oil recovery techniques, such as gas injection and in-situ combustion, as the need to access a wider array of hydrocarbon resources increases. Along with this, CO2 sequestration as a carbon storage mechanism is increasingly important. These subsurface flow problems can be physically unstable at the scale modeled. This means that the displacement front will generally exhibit an unstable growth pattern known as viscous fingering. When simulating these processes, numerical errors which can be strongly correlated with the underlying computational grid can trigger, or at least bias, the formation of these fingers in ways not dictated by the underlying rock properties. This can give rise to the so-called grid orientation effect, where varying the orientations of the computational grid results in convergence to fundamentally different solutions.;The interest of this thesis is the development of numerical techniques that reduce grid dependency for these processes. We do this in a manner that draws upon fundamental numerical analysis techniques in conjunction with a physical intuition about the underlying governing equations. The governing equations for subsurface flow are of a mixed elliptic-hyperbolic character; here we focus on the hyperbolic portion but seek methods compatible with general elliptic discretizations and computational grid topologies.;Some aspects of the grid dependency can be traced back directly to the handling of the injection wells. Perturbation created early in time and near the injection wells propagate into the interior of the domain setting up an initial biasing for the simulation. We demonstrate that the handling of injection wells drastically affects the final result. We propose two methods, a well-sponge method and a local embedding technique, which extends the well region based on near-well models. Both of these methods drastically increases the similarity between the solutions on different computational grids.;The numerical method used in the interior region can also have a significant impact on the computed solution. For multidimensional (multi-D) problems the structure of the numerical diffusion tensor creates preferential flow directions. These preferential flow directions can be strongly coupled to the grid, especially when a 1-D method is applied dimension by dimension. This motivates the development of truly multi-D transport discretizations that incorporate local flow information into the flux calculation. We develop a family of compact, upstream biased multi-D finite volume methods for 2-D simulation. The methods control the structure of the diffusion tensor and incorporate characteristic flow information through the use of a local coupling methodology based on interaction regions. The family of schemes is positive for linear equations and monotone for a class of scalar nonlinear equations.;Two multi-D schemes, tight multi-D upstream weighting and smooth multi-D upstream weighting, are proposed and analyzed. Since transverse diffusion largely controls finger formation in linear problems, a third scheme is proposed called the Flat Scheme, which has constant transverse diffusion. All the multi-D methods are appropriate for both explicit and implicit time stepping. The methods are applied to both miscible and immiscible displacement problems and significantly reduce biasing due to the grid as compared with the commonly used single point upstream weighting method.