Numerical Methods

Key Takeaways

• Finite difference methods (FDM) discretize the spatial and temporal derivatives in the governing equations using Taylor series approximations on a regular grid, replacing the continuous differential equation with a system of algebraic equations at each grid node: the most common FDM scheme for the reservoir flow equation uses a first-order backward difference in time (implicit scheme, unconditionally stable) and central differences in space, yielding the standard five-point (2D) or seven-point (3D) stencil used in all major commercial reservoir simulators (Eclipse, CMG IMEX, Nexus); the IMPES (implicit pressure, explicit saturation) formulation solves the pressure equation implicitly at each time step and updates saturations explicitly, which is computationally efficient but requires small time steps for stability; the fully implicit (FI) formulation solves all variables simultaneously, allowing larger time steps at the cost of solving a larger, coupled linear system at each time step; adaptive implicit methods (AIM) apply the fully implicit scheme only in grid cells where the flow is rapid (near wells, faults, high-permeability streaks) and the IMPES scheme elsewhere, balancing efficiency and stability; grid orientation effects (numerical artifacts that cause flow to align preferentially along grid axes rather than following the true pressure gradient direction) are a fundamental limitation of FDM on Cartesian grids and are partially mitigated by using nine-point or MPFA (multi-point flux approximation) spatial differencing schemes.
• Finite element methods (FEM) use piecewise polynomial basis functions defined on unstructured triangular or tetrahedral meshes to approximate the solution of the governing equations, providing greater geometric flexibility than FDM (arbitrary domain shapes, local mesh refinement in zones of high gradient) and superior accuracy for problems with complex boundaries such as wellbore stress analysis, hydraulic fracture propagation, and subsidence modeling: the FEM solution is obtained by minimizing the weighted residual of the governing equation over each element (the Galerkin method) and assembling the element contributions into a global stiffness matrix that is solved by direct or iterative linear algebra; geomechanical finite element models (using software such as ABAQUS, DIANA, or ELFEN) are used to analyze wellbore stability (computing the stress concentration around the borehole in anisotropic stress fields), casing integrity (computing the bending and buckling response of casing in subsiding formations), hydraulic fracture geometry (computing the coupled fluid flow and rock deformation during fracture propagation), and reservoir compaction (computing the vertical and lateral strain in the reservoir and overburden in response to production-induced pore pressure decline); coupled FEM-FDM approaches (using FEM for the geomechanics and FDM for the fluid flow) are increasingly used for problems where pore pressure and stress are tightly coupled, such as chalk reservoir compaction (Ekofisk), unconsolidated sand production, and CO2 injection-induced seismicity assessment.
• Finite volume methods (FVM) are the dominant discretization approach for computational fluid dynamics (CFD) problems in petroleum engineering, including wellbore multiphase flow (slug flow, annular flow, churn flow regimes in tubing and pipelines), separator design (liquid holdup, gas-liquid separation efficiency), and subsea pipeline flow assurance (hydrate formation, wax deposition, severe slugging prediction): FVM divides the computational domain into control volumes (cells) and applies the conservation laws (mass, momentum, energy) in integral form over each cell, with fluxes computed at cell faces and accumulated to update cell-averaged quantities; the FVM approach naturally conserves mass and momentum (unlike FDM which may not be strictly conservative on non-uniform grids) and can handle complex geometries using body-fitted curvilinear grids or unstructured meshes; commercial flow assurance tools (OLGA, LedaFlow, Pipesim) use FVM or variant schemes to simulate transient multiphase pipeline flow, with the governing equations including gas-liquid momentum exchange, slip velocity correlations, and wall friction models calibrated against field data; reservoir simulation on corner-point geometry (CPG) grids -- the industry standard grid format for representing faulted, folded reservoir geometries in commercial simulators -- uses a variant of FVM with transmissibility-weighted inter-cell fluxes to handle non-orthogonal cell connections at faults and unconformities.
• Linear and nonlinear systems arising from the discretization of the governing equations must be solved efficiently to make simulation of large reservoir and wellbore models computationally tractable: a modern full-field black-oil reservoir simulation model with 5 to 50 million grid cells generates a linear system of equations with 5 to 150 million unknowns at each time step, which must be solved in seconds to minutes for the simulation to complete in practical run times; direct solvers (Gaussian elimination and its variants, LU factorization) are exact but scale poorly with problem size (O(n^3) operations for a dense matrix, reduced but still unfavorable for sparse matrices by fill-in during factorization); iterative solvers (Krylov subspace methods including conjugate gradient, GMRES, and BiCGSTAB preconditioned by incomplete LU factorization, constrained pressure residual (CPR) preconditioning for pressure-saturation systems) scale approximately linearly with the number of unknowns for well-structured problems and are the dominant approach in modern reservoir simulators; parallel computing (distributing the grid cells across multiple CPU cores or GPU processors using domain decomposition) reduces wall-clock time by factors of 10 to 1,000 for large models, with commercial simulators (Eclipse, CMG, Intersect) offering parallel scaling up to several thousand CPU cores on high-performance computing clusters; GPU acceleration for the linear algebra step (using the CUDA or OpenCL frameworks) provides additional speedup factors of 5 to 20 for the matrix-vector products that dominate Krylov solver cost.
• Uncertainty quantification in numerical models requires running large ensembles of simulations with different parameter realizations to characterize the range of possible subsurface outcomes: a reservoir model constrained by seismic, well log, and production data has uncertain parameters including permeability distribution, porosity, fault transmissibility, relative permeability endpoints, and fluid contact depths, each of which can take a range of plausible values consistent with the available data; Monte Carlo simulation (random sampling from the parameter distributions and running a forward simulation for each sample) provides the most straightforward uncertainty quantification but requires thousands to tens of thousands of simulation runs to adequately sample high-dimensional parameter spaces; proxy models or emulators (fast surrogate models trained on a small number of simulation runs that can predict the output for new parameter combinations without running the full simulator) are used to reduce the computational cost of ensemble-based uncertainty quantification, with machine learning approaches (Gaussian processes, neural networks trained on simulation output) increasingly replacing the polynomial response surfaces used in earlier generations of uncertainty quantification workflows; ensemble Kalman filter (EnKF) and related data assimilation methods update the ensemble of models sequentially as new production data becomes available, providing a computationally tractable approach to continuous history matching and uncertainty tracking in producing fields.