CSMP++ simulation approach

The main technical deficiencies of commercial reservoir simulation tools with respect to the realistic simulation of multiphase flow and geomechanic processes in the subsurface are:

• Lack of geometrical flexibility to capture geological detail

  For current engineering practice this implies that many important geologic features and structural details of the storage complex must be ignored to obtain a workable simulation model. Key obstacle: the use of structured grids. These dictate a uniform resolution and regularized geometry that makes it impossible to represent oblique large aspect ratio features such as faults, layer pinch-outs, etc (Matringe et al, 2008). Furthermore the grid cannot be refined to resolve local details of interest and / or achieve a uniform distribution of discretization errors.

• Oblique faults cannot be represented

  This is the well-known problem of non-neighbour connections. Displaced or truncated layers usually are represented by offsets in the structured grid, and the grid is aligned with such boundaries. Thus, faults as such have do not have a discrete grid representation. Cross-flow or flow in the fault plane can only be approximated with ad-hoc cell connections. When inclined faults are represented using stair-step arrays of grid-blocks, this fails to create oblique flow continuity because the finite-difference stencils only consider transmissibilities along the principal axes of the grid. Finally, unless a prohibitively fine resolution is used, fault thickness tends to be over represented.  If so, this distorts flow velocities and leads to a different flow-focussing potential.

• Implicit wells

  The lack of spatial adaptivity and refinement also precludes a discrete representation of wells in conventional simulation models. In stead, semi-analytical well representations are used to relate source / sink terms to corner point pressures. For wells in structurally complex reservoirs, this treatment is inadequate (Durlofsky et al, 2000). In addition, it precludes model inclusion of near-wellbore flow physics.

• Key phenomena / features of the flow are not resolved

  The inability to refine regions of interest implies that most reservoir simulation models locally are insufficiently refined.  If highly permeable strata are represented by a single layer of grid-blocks, internal processes that have a decisive impact on production, like gravity-tonguing or override, will not show. In this case, also the unique properties of the layer will be smeared along the vertical grid axis because of harmonic averaging with adjacent layers in the transmissibility calculations. Furthermore, standard hydrocarbon reservoir simulators will produce a cylindrical plume irrespective of injection rate and viscosity contrast with the formation water. In reality, viscous fingering would occur at high flow rates. By contrast with laboratory observations, constitutive saturation functions are treated as rate independent and will therefore only capture a subset of the CO2 dynamics. If these deficiencies are not addressed, the benefits of a better geologic characterisation of the Australian subsurface will not be realised.

• Limitations of 2-point flux approximations

  Although higher-order FD stencils and/or so-called multipoint flux approximations are available in many commercial reservoir simulators these are difficult to construct near material interfaces and their application slows down simulation. Engineers therefore use 2-point flux approximations, accepting severe grid-orientation effects (Chen et al., 1993), suppressed flow focussing, an unrealistically high sweep and artificially stable displacement fronts. For the accurate prediction of the unstable spreading of a CO2 plume, this lack of fidelity in the representation of the flow processes is intolerable.

• Smearing of material interfaces

  Classic (point-based) finite difference or finite volume methods require averaging of flow properties across material interfaces. Not only does this blur the representation of material interfaces in the model. If different flow physics apply on opposite sides of such boundaries, these averages are physically meaningless.

• Globally driven time-stepping

  In conventional reservoir simulation, the so-called Courant-Fredrich-Levy condition dictates the size of time steps (the discretization of time continuous time). However, the size of optimal time steps decreases with flow rate and the latter varies over many orders of magnitude within a single reservoir: Saturation fronts in the far field move at a speed of tens of centimetres per year, while velocities near the wellbore approach meters per second. Therefore any global time stepping scheme cannot be optimal. In addition, only a few tens of years of behaviour need to be simulated / forecast for hydrocarbon reservoirs, while thousands of years are mandatory for CO2 storage complexes.

• Inability to do multiphysics

  Since there is no single software that solves multiphase advection, diffusion, and mechanics problems equally well, file-based couplings between different software packages are commonly used to implement multiphysics simulations. However, such loose sequential coupling strategies are ill-suited to capture the nonlinear dynamics that ensue from the THMC couplings. Multi-application approaches also suffer from inconsistencies, limited predictive capabilities, and a lack of robustness.

• Goal-based simulation impossible

  Any numerical estimate of a physical quantity should come with an error bar and an uncertainty quantification accounting for the natural variability and limited knowledge of the input parameters. However, numerical error estimates for water breakthrough or recovery rate are rarely presented in reservoir simulation studies because errors cannot be quantified for complex simulations where gradients vary over time and the grid might be optimal for one of the equations that are solved, but not for the others. It follows, that the simulation engineer cannot prescribe target accuracy, i.e. set prediction goals determining the computational effort. Such goal-based simulation is highly desirable and already the standard in Computational Fluid Dynamics where unstructured grids are adapted dynamically so that emergent subgrid-scale features are captured. Goal-based simulation carries a high potential for reservoir simulation but has yet to be introduced in this domain.

• Limited parallelism

  Moore’s law no longer applies to CPU clock speed. Thus, dramatic speed gains are possible only by making use of improved hardware parallelism.  For billion-cell geologically realistic models, short runtimes can only be achieved if they are implemented in a highly granular way on dedicated massively parallel hardware like GPUs. This has not happened yet.

  While conventional simulators are highly optimized to minimize computational cost, and although they can match reservoir pressure and production responses after calibration with historic data, their lack of physical realism and spatio-temporal resolution makes it questionable whether they are suited for the long-term prediction of the performance of CO2 storage complexes.

  The alternative approach pursued with CSMP++ to overcome the enlisted shortcomings of conventional reservoir simulation and the underpinning improved simulation technology are summarised in Table 1.

• Table 1

  Limitations of commercial reservoir simulation software and solutions pursued with CSMP++

  Simulation limitations	How we can overcome them
  Lack of geometric flexibility to represent faults, pinch-out and large-aspect ratio features	Unstructured grids directly from structural framework model, grid-to-grid conservative property mapping
  Implicit wells	Discrete lower-dimensional well representations and 1D-well flow equations
  2-point flux approximations	Multi-dimensional fluxes with compact support from finite-element framework
  Slow serial computation	GPU parallelised simulation kernel
  Slow global time-stepping	Asynchronous spatially adaptive time stepping
  Inability to do multiphysics in a single software	Single THMC framework that combines geometric flexibility of FEM with local conservativeness of FVM methods
  No estimation of discretisation errors	Hessian-based posterior error estimates supporting dynamic mesh adaptation
  Slow linear algebraic solvers	Parallelised Algebraic Multigrid for Systems (SAMG) solves couple geomechanics-pressure equations. CUDA-parallelised linear solvers
  No reactive transport calculations	CSMP-GEMs runtime integration facilitating computation of rate-limited solid fluid chemical equilibration

  Much of the functionality enlisted in Table 1, already exists in Complex Systems Modelling Platform (CSMP++). This is a C++14 ANSI / ISO standard-compliant library of numerical methods for the combined simulation of thermal – hydrological – mechanical and chemical (THMC) processes designed by the lead investigator (e.g., Matthai et al, 2007), including black-oil (Geiger et al, 2009) and CO2 injection functionality (Annewandtner, Main and Geiger, 2013). There are >50 ISI peer-reviewed applications that illustrate the functionality of this code, some accompanied by journal covers and publication highlights. Two simulation studies were published in Science (Coumou et al, 2008; Weiss et al, 2012). CSMP++ contains implementations of the finite-element- (FEM) and finite volume (FVM) methods and combinations thereof (Geiger et al, 2004, 2009; Paluszny et al, 2007, Bazrafkan et al, 2014). In addition to solution methods for elliptic-parabolic and hyberbolic partial differential equations, CSMP++ provides implementations of constitutive relationships for THMC processes, equations of state (Driesner and Heinrich, 2007; Matthai and Nick, 2009), and property modelling and simulation analysis functionality (e.g., Mosser 2013). It also integrates Fraunhofer’s SAMG algebraic multi-grid solver for systems of coupled equations (Stüben 2001) and it is based on a unique data model: There are neither arrays nor cell numbers. Material properties and dependent variables are stored directly on a hierarchical tree structure representing the mesh. This implementation supports the proposed local time-stepping methods, goal-based simulation and dynamic mesh adaptivity which has already been applied to model fracture propagation (Paluszny and Matthai, 2010). Multi-domain simulations where different physics are modelled in different regions have been implemented (Milliotte et al, 2013; Matthai and Milliotte, 2014). CSMP++ also forms the basis of a wide variety of simulation programs, from discrete fracture and matrix modelling (Matthai et al, 2012), upscaling and reservoir simulation (Matthai and Nick, 2009) to flow processes with phase separation at extreme conditions such as magmatic-hydrothermal systems (Weiss et al, 2014).