CSMP++ simulation approach
... and its benefits for reservoir simulation and CO2 geostorage performance prediction.
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 pinchouts, 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 wellknown problem of nonneighbour 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. Crossflow or flow in the fault plane can only be approximated with adhoc cell connections. When inclined faults are represented using stairstep arrays of gridblocks, this fails to create oblique flow continuity because the finitedifference 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 flowfocussing potential.

Implicit wells
The lack of spatial adaptivity and refinement also precludes a discrete representation of wells in conventional simulation models. In stead, semianalytical 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 nearwellbore 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 gridblocks, internal processes that have a decisive impact on production, like gravitytonguing 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 2point flux approximations
Although higherorder FD stencils and/or socalled 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 2point flux approximations, accepting severe gridorientation 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 (pointbased) 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 timestepping
In conventional reservoir simulation, the socalled CourantFredrichLevy 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, filebased couplings between different software packages are commonly used to implement multiphysics simulations. However, such loose sequential coupling strategies are illsuited to capture the nonlinear dynamics that ensue from the THMC couplings. Multiapplication approaches also suffer from inconsistencies, limited predictive capabilities, and a lack of robustness.

Goalbased 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 goalbased simulation is highly desirable and already the standard in Computational Fluid Dynamics where unstructured grids are adapted dynamically so that emergent subgridscale features are captured. Goalbased 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 billioncell 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 spatiotemporal resolution makes it questionable whether they are suited for the longterm 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, pinchout and largeaspect ratio features Unstructured grids directly from structural framework model, gridtogrid conservative property mapping Implicit wells Discrete lowerdimensional well representations and 1Dwell flow equations 2point flux approximations Multidimensional fluxes with compact support from finiteelement framework Slow serial computation GPU parallelised simulation kernel Slow global timestepping 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 Hessianbased posterior error estimates supporting dynamic mesh adaptation Slow linear algebraic solvers Parallelised Algebraic Multigrid for Systems (SAMG) solves couple geomechanicspressure equations.
CUDAparallelised linear solvers
No reactive transport calculations CSMPGEMs runtime integration facilitating computation of ratelimited 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 standardcompliant 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 blackoil (Geiger et al, 2009) and CO2 injection functionality (Annewandtner, Main and Geiger, 2013). There are >50 ISI peerreviewed 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 finiteelement (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 ellipticparabolic 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 multigrid 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 timestepping methods, goalbased simulation and dynamic mesh adaptivity which has already been applied to model fracture propagation (Paluszny and Matthai, 2010). Multidomain 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 magmatichydrothermal systems (Weiss et al, 2014).