Carbon capture and sequestration initiatives make new demands on modern reservoir simulators. To find optimal locations and volumes of CO2 to inject into a subsurface to maximize CO2 storage, we must simulate a large ensemble of injection cases. One possible solution to the computational complexity of this task is to employ machine-learning models which, after a one-off overhead cost of training, can infer and predict future states of a reservoir several orders of magnitude faster than traditional methods. Most previous work in the literature has primarily focused on either convolution-based methods or, more recently, neural operator-based methods, to predict the evolution of state variables. These architectures have shown promise in predicting on structured reservoir grids but lack the capability to extend the same level of accuracy to unstructured grids. Graph neural networks (GNNs) overcome this bottleneck by incorporating inductive biases arising from local message-passing mechanisms, facilitating convolution operations over complex graphs and meshes. In this work, we present a novel autoregressive GNN autoencoder to predict time-varying state variables for an ensemble of CO2 injection cases. We implement a graph convolution network for the message-passing protocol and incorporate physics-informed edge weights between cell connections to guide flow. An exhaustive set of node features are used to train the model on the hyperbolic evolution of phase saturations while preserving the ellipticity in the pressure. We test the performance of the GNN model for (1) its ability to predict state variables for varying injection rates of CO2, (2) for the post-injection phase, and (3) under different unseen geological configurations. Training and testing are performed by constructing ensembles of 2D, 3D, and real field cases that best represent these scenarios. For the 2D regular grid case, we observe that the model can capture pressure and saturation values accurately, even for highly varying injection rates and with only a limited amount of data. This performance is maintained in the post-injection phase. A key advantage of GNNs is that they show a distinct ability for transfer learning on ensembles of unseen geological configurations. We observe that the model can predict the shape and intensity of wavefronts of certain cases with no prior exposure to the specific static properties during training. Similar results are produced for 3D grids and real field cases.

Commercial reservoir simulators must be very robust and fast. Moreover, current hardware requires the simulators to scale over multiple number of computing nodes and for a fixed (‘strong scalability’) as well as an increasing problem size per computing node (‘weak scalability’). In most current commercial reservoir simulators, due to the different geological structures and properties of hydrocarbon reservoirs and the use of enhanced oil recovery (EOR) techniques, the governing equations are strongly nonlinear and hard to solve. The Jacobian system is solved by FGMRES preconditioned by the two-level constrained pressure residual (CPR) preconditioner. The driving force of the CPR preconditioner is the solution of the pressure equation. The industry standard for solving the pressure equation is the algebraic multigrid (AMG) solver. AMG is well known for its ‘weak scalability’. However, in these applications, AMG has unfavorable ‘strong’ scalability properties. This degradation in scalability is due to the increased level of inter-processor communication in the algorithm. In this paper, a monotone non-Galerkin AMG (MNG-AMG) method is presented. The aim of the method is to reduce the overall communication in MNG-AMG by enforcing a predefined nonzero pattern and monotonicity property (i.e., M-matrices) on each multigrid level. This paper describes the application of the MNG-AMG method in the context of reservoir simulations. We will compare the parallel scalability of the default solver with the MNG-AMG solver and discuss the optimal values for the MNG-AMG solver for a variety of test cases based on full field reservoir simulations.

Extreme and isolated eigenvalues are known to be harmful to the convergence of an iterative solver. These eigenvalues can be produced by strong heterogeneity in the underlying physics. We can improve the quality of the spectrum by ‘deflating’ the harmful eigenvalues. In this work, deflation is applied to linear systems in reservoir simulation. In particular, large, sudden differences in the permeability produce extreme eigenvalues. The number and magnitude of these eigenvalues is linked to the number and magnitude of the permeability jumps. Two deflation methods are discussed. Firstly, we state that harmonic Ritz eigenvector deflation, which computes the deflation vectors from the information produced by the linear solver, is unfeasible in modern reservoir simulation due to high costs and lack of parallelism. Secondly, we test a physics-based subdomain-levelset deflation algorithm that constructs the deflation vectors a priori. Numerical experiments show that both methods can improve the performance of the linear solver. We highlight the fact that subdomain-levelset deflation is particularly suitable for a parallel implementation. For cases with well-defined permeability jumps of a factor 104or higher, parallel physics-based deflation has potential in commercial applications. In particular, the good scalability of parallel subdomain-levelset deflation combined with the robust parallel preconditioner for deflated system suggests the use of this method as an alternative for AMG.