Voronoi Meshing to Accurately Capture Geological Structure in Subsurface Simulations

Mesh generation lies at the interface of geological modeling and reservoir simulation. Highly skewed or very small grid cells may be necessary to accurately capture the geometry of geological features, but the resulting poorly scaled or small grid cells can have a substantial negative impact on simulator accuracy and speed. One way to minimize numerical errors caused by gridding complex structures is to simulate on high-quality Voronoi meshes, which reduce grid orientation effects in fluid flow. This work presents a complete methodology to create Voronoi simulation grids, model fluid flow in complex geological systems, and visualize the results. A recently developed Voronoi meshing method that can automatically generate provably good unstructured meshes that conform to input surfaces creating closed volumes is used. Initially an analytical benchmark simulation is presented to validate the quality of the meshes and simulation results and demonstrate the superiority of simulation results using Voronoi meshes over flexed-hexahedral meshes on a domain with internal features. Next, meshes are created for test structures representing four of the most common geological features in the subsurface: layering, pinch-out, an interior lens that tapers to zero thickness on all sides and a fault with offset. Two benchmark flow simulations are run for each test structure. Finally, a realistic geological example for CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document} injection into an anticline is simulated. Three realizations of the Voronoi mesh at the same resolution are generated for the simulations. Each mesh is highly refined near the injection wells and coarse in areas of less interest. These three meshes are used to model the CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document} plume in the subsurface as it migrates to the top of the structure and then fills downward. Simulations on the meshes with randomly generated elements inside the input volumes each give slightly different fingering patterns for the viscous-unstable buoyant gas flow. The results presented in this work show a promising step towards utilizing fully automated Voronoi meshing for subsurface flow simulations in complex geology.

then fills downward. Simulations on the meshes with randomly generated elements inside the input volumes each give slightly different fingering patterns for the viscousunstable buoyant gas flow. The results presented in this work show a promising step towards utilizing fully automated Voronoi meshing for subsurface flow simulations in complex geology.

Introduction
Simulation mesh development is an often overlooked aspect of moving from geological modeling to reservoir simulation. Hexahedral meshes can be flexed to capture the geometry of geological features within the simulation domain, but grid cells near the features may become poorly scaled (Berge 2016). These non-orthogonal grid cells can have a substantial negative impact on simulator accuracy (Wu and Parashkevov 2009) and small grid cells can slow down simulations. Without orthogonality, the direction of the calculated flux between two grid cell centers is not perpendicular to the cell faces and there is an error in the flux calculation. These errors can accumulate, resulting in large errors in the simulation result even in relatively simple reservoir geometries (Berge 2016;Kozdon et al. 2011;Lunde 2007;Palagi and Aziz 1994).
One way to minimize numerical errors introduced by the grid is to simulate on Voronoi meshes. Voronoi meshes are extremely advantageous for meshing complex geometries for subsurface simulations as there is no constraint on cell topology, so any closed volume can be meshed to the desired degree of accuracy with sufficient grid refinement. By construction, in a Voronoi mesh the direction of the calculated flux between two adjacent cells is perpendicular to the element face connecting them. Thus, an important source of numerical error introduced by using flexed hexahedral or tetrahedral meshes is minimized. Voronoi meshes are particularly promising for simulators that use two-point flux approximations (TPFA), also known as single-point upstream weighting (SPU) as they are more susceptible to errors on non-orthogonal meshes than other flux approximation methods (Lunde 2007). An unstructured mesh also limits grid orientation effects, which can be an issue in unstable displacements such as buoyant CO 2 migration, regardless of the flux approximation used (see, e.g., Kozdon et al. 2011).
Despite the fact that errors caused by the resolution of the mesh or the way it represents the geological structure can be significant, it is often omitted from sensitivity analysis/uncertainty quantification (SA/UQ) studies because meshing error has historically been difficult and time-consuming to analyze in subsurface simulations. However, mesh uncertainty could, and likely should, be assessed as part of the workflow that is used to study uncertainty in material and fluid properties.
Voronoi meshes for subsurface flow simulation have been studied intermittently since the 1990s. Studies show that simulations on two-dimensional Voronoi meshes have high accuracy compared to other meshes (Palagi and Aziz 1994;Du and Gunzburger 2002). The challenges faced by most groups are the difficulty in creating a Voronoi mesh for a realistic level of geological complexity and the additional computational overhead required to simulate on a Voronoi mesh.
Voronoi grids are also known as perpendicular bisection (PEBI) grids, particularly in the petroleum industry literature (Heinemann et al. 1989;Palagi and Aziz 1994). Many of the PEBI grids in the literature are not genuinely Voronoi in three dimensions, but are two-dimensional Voronoi meshes extruded into the third dimension. If the extrusion is perpendicular to the plane of the two-dimensional mesh, then the resultant 2.5-dimensional mesh is Voronoi. However, if the extrusion is flexed to fit non-planar surfaces or extruded to a plane that is not parallel to the two-dimensional mesh, the final mesh is not Voronoi. The calculated fluxes between adjacent grid cells will not all be perpendicular to grid cell faces, though they may be close, and the mesh may still be far superior to a traditional hexahedral mesh for simulation. However, the accuracy may be dependent on both the geometry of the simulation domain and the direction of flow (Wu and Parashkevov 2009).
There is a substantial body of literature on extruding two-dimensional Voronoi to 2.5-dimensional or three-dimensional meshes for subsurface flow simulation (see, for example, Berge (2016), Berge et al. (2019), Branets et al. (2008), Deng et al. (2011), Fung et al. (2014, and Lunde (2007)). Notably, Branets et al. (2008) generated twodimensional Voronoi meshes for a layer of the highly heterogeneous SPE10 (Christie and Blunt 2001) benchmark case, and demonstrated that the resultant simulations have similar accuracy as simulations on the original mesh, but with far fewer grid cells. Lunde (2007) demonstrated that TPFA may not be accurate on unstructured threedimensional meshes, but with 2.5-dimensional Voronoi meshes, the error in TPFA simulations is of the same order as the simulated error for mimetic finite-differences discretization, which is suited for non-orthogonal unstructured meshes. Berge (2016) and Berge et al. (2019) created two-dimensional meshes that conform to wells and complex fault networks, though fault geometry is only approximately preserved near intersections. These meshes are extruded into 2.5-dimensional models and shown to have reduced grid orientation effects compared to traditional mesh flexing.
There are fewer studies of fully three-dimensional Voronoi meshes, due to the challenge of creating Voronoi meshes in three dimensions. For example, Merland et al. (2014) and Bonduà et al. (2017) presented three-dimensional meshes that are strictly Voronoi; however, they may both approximate interior boundary locations. Bonduà et al. (2017) accurately modeled gas flow in a complex reservoir structure in a hybrid Voronoi-hexahedral mesh with far fewer cells than a comparable hexahedral mesh. Freeman et al. (2014) presented the development of a Voronoi meshing software, MeshVoro, for the TOUGH family of codes. Their example problems include complex geological models and well configurations. They note that their algorithm may struggle to capture two mesh objects such as faults and wells that are very close to each other. In their discussion, Freeman et al. (2014) also highlight the fact that the Jacobian matrix required to simulate on Voronoi meshes is intrinsically less sparse and less well-conditioned than the Jacobian of a structured mesh. This results in simulations that will take longer for the same number of grid cells, which is also seen in the workflow presented in this study (see section 6.3.3 of LaForce et al. 2020). Freeman et al. (2014 emphasize that Voronoi meshes should only be used when they are suited to the demands of the problem and that the power of Voronoi meshes is in applications where it is not possible to represent the geometry efficiently using Cartesian meshes. The same holds true for the present work. All of the software in the present work is free and open-source except for the Voronoi meshing software, and an open-source version of it should be released this calendar year. In previous studies, Gross et al. (2019), Sevougian et al. (2019) and LaForce et al. (2020) detailed work using earlier versions of the VoroCrust software for subsurface flow calculations. Gross et al. (2019) generated VoroCrust (Abdelkader et al. 2020) meshes for four test structures representative of common geological features. At that time, VoroCrust output meshes did not contain all the information needed by flow simulators, so Voronoi meshes created in LaGriT (Los Alamos National Laboratory 2017) were used for 16 benchmark simulations in PFLOTRAN (Hammond et al. 2014;Lichtner et al. 2020) and FEHM (Los Alamos National Laboratory 2018). PFLOTRAN simulations on VoroCrust Voronoi meshes of these four test structures were later studied by Sevougian et al. (2019). LaForce et al. (2020) presented PFLO-TRAN simulations on VoroCrust meshes for two analytical benchmark cases and two large-scale models using surfaces from geological mapping software.
In the present work, meshes are generated using VoroCrust (Abdelkader et al. 2020), a fully automated Voronoi meshing software that has been adapted to create simulation grids for fluid flow in complex geological systems. This software has the ability to create a mesh that satisfies the Voronoi/Delaunay criteria and conforms exactly to difficult material geometries for non-manifold surfaces. Much like Freeman et al. (2014) and Bonduà et al. (2017), the present work uses Voronoi meshes to improve the accuracy of simulations using the TPFA for complex reservoir structures. VoroCrust generates meshes based on closed volumes with watertight input surfaces.
All simulations in the present work use PFLOTRAN (Hammond et al. 2014;Lichtner et al. 2020), which is an open-source, community-developed software with a broad range of capabilities for one-and two-phase fluid flow and chemical transport through porous media that may also include thermal impacts. In this work, thermal diffusion in a water-saturated porous medium, single-phase, single-phase variably saturated, and two-phase flow with chemical component partitioning between the fluid phases are simulated. PFLOTRAN has been shown to scale well to thousands of processors (Hammond et al. 2014). PFLOTRAN uses the finite-volume approximation for fluid flow and TPFA. Visualization is done in ParaView, which is also open-source software (Ahrens et al. 2005).
Simulations are initially benchmarked against a two-dimensional analytical solution for heat diffusion that has been modified from LaForce et al. (2020), and their accuracy is compared against simulations on a flexed hexahedral mesh. Four common reservoir structures are then meshed for further benchmarking. The four geometries are layering, pinch-out of a layer, an interior lens with pinch-out on all sides, and offset fault, and were originally presented in Gross et al. (2019). Finally, a realistic example for CO 2 injection into a doubly plunging anticline is shown. The large-scale surfaces are the same as used by Stauffer et al. (2009), Deng et al. (2012 and Harp et al. (2017), though a different part of the structure and simpler flow parameters are chosen. Fig. 1 A two-dimensional illustration of the VoroCrust algorithm: First, we cover the input geometry with a set of overlapping spheres (circles in two dimensions). We place seeds at the uncovered intersection points of each sphere triplet (pair in two dimensions). The resulting Voronoi cells split the domain around the geometry into two subregions bounded by an interface of Voronoi facets that represents the input geometry. Finally, we carefully add seeds in the interior to generate cells with good aspect ratios 2 Voronoi Meshing

VoroCrust Meshing Software
The main challenge in generating Voronoi meshes for complex geometric models is how to force the Voronoi cells to conform to curved surfaces, which may include sharp curves and/or corners. VoroCrust (Abdelkader et al. 2018;Abdelkader et al. 2020) has a novel automated solution to this problem. The way the algorithm works is to completely cover the input surface with a set of spheres. The radius of each sphere is selected automatically based on the geometric features associated with the input geometry. Then Voronoi seeds are placed at every intersection point of an overlapping sphere triplet that is not covered by a fourth sphere. This forces the triangle connecting the centers of the three spheres to form a Voronoi face in the output Voronoi tessellation. Repeating this for all the sphere triplets results in a set of Voronoi faces that are a provably good Boissonnat and Oudot (2005) representation of all input interfaces regardless of their geometric/topological complexity, even when the regions to be meshed are very thin.
The software analyzes the input to detect sharp features and the curvature of the input interfaces. During this process, a sizing function is constructed to maintain the geometric features of the input model. Handling sharp corners and curves is discussed in more detail in Abdelkader et al. (2020). Once the surface seeds are generated, additional seeds are added to the interior of the domain to ensure an upper bound on the aspect ratio of each Voronoi cell in the output. The resulting Voronoi mesh has grid cells that always have triangular faces on input surfaces, but are general polyhedra in the interior. An illustration of the VoroCrust algorithm is given in Figs 1 and 2.
The published results in Abdelkader et al. (2020) focused on the provable correctness of the method. The published algorithm was slow, requiring hours in some cases to generate a Voronoi mesh with a million cells. A variety of algorithmic changes have since been made, including multithreading, which has led to improved meshing speed. The current code can generate a one-million-cell mesh in less than 10 min on a 24-core dual-threaded Linux box in an automated process that detects all features on its own and generates a provably good mesh at the end.
Another advantage of VoroCrust is that it relies on a random sphere packing process. Hence, every time the code is executed for an input file, it generates a different  2021) The most recent development is to generate the output Voronoi mesh as an Exodus II file (Schoof and Yarberry 1994).
If the physics of the simulation require additional refinement, the user can impose a smaller sizing function requirement on top of the geometric sizing requirements. The refinement imposed by the user can be global or point-wise, as needed for the simulation.
Finally, the user can specify a set of "monitoring points" to include as seeds in the Voronoi mesh. The constraint there is that none of these seeds can interfere with the surface seeds that are generated to capture an interface. In the case that an input seed is too close to the surface (e.g., it lies inside one of the surface spheres), this seed will be discarded and not included in the final mesh. If the domain has a narrow region and the input sizing function does not allow several layers of Voronoi cells to capture that region, it can be challenging to place monitoring points in these locations. This can sometimes be alleviated by defining a small sizing function and is an area of future work.

Mesh Generation
The process for creating a Voronoi mesh and converting it for use in PFLOTRAN is discussed in detail in Section 6.3.2 and Appendix C of LaForce et al. (2020) and is only briefly summarized here. VoroCrust accepts watertight Wavefront Object (.obj) files as input surfaces for meshing. It can also read in two optional comma-separated variable (.csv) files. The first file, monitoring_points.csv, is a list of the desired monitoring points in the domain. They are used as initial seeds, but sometimes during the meshing process these points are eliminated, as discussed in the preceding section. A second optional file of points where the grid is to be refined, regions.csv, can also be read in. The maximum grid cell radius at any point in the domain is defined based on the nearest user-specified grid refinement point.
The Wavefront Object surface files are created by open-source software LaGriT available at https://github.com/lanl/LaGriT, (Los Alamos National Laboratory 2017), which was developed to create meshes for geological applications. LaGriT reads in surfaces created by geological mapping software and creates simulation meshes based on these files, or can be used independently to create meshes for geometric shapes. Both capabilities are used here. LaGriT is used to output Stereolithography (.stl) files containing only the external and internal surfaces of the model. These files are converted to .obj files using ParaView, available at https://www.paraview.org (Ahrens et al. 2005). Finally, the surface file and optional monitoring point and refinement point files are read into VoroCrust, which creates meshes in an automated process.
The mesh is output in two formats: a mesh.vcg file that contains only the information needed to simulate using the finite-volume method, and an Exodus II file that will later be used for visualization of simulation results. A Python script (voro2pflotran.py) is used to split the information in the mesh.vcg file into three types of files that PFLO-TRAN needs for simulations; a mesh file, a material ID file that will be used to define properties for each region in the domain, and boundary files (for rectangular domains only). This script is available in the PFLOTRAN source code repository at https:// www.pflotran.org in the folder pflotran/src/python/vorocrust.

Postprocessing
PFLOTRAN output files on Voronoi meshes do not contain sufficient information to create three-dimensional visualizations. A Python post-processing script (pflo-tran2exodus.py) is used to add all simulation outputs to the Exodus II mesh file. This script is also available in the PFLOTRAN source code repository. The combined file is read in ParaView for three-dimensional visualization.

Benchmark Model
One analytical benchmark simulation is considered in the present work. Simulations on a Voronoi mesh and a flexed-hexahedral mesh are compared to demonstrate the effect of mesh flexing on the local accuracy of finite-volume simulations and to show that simulations on the Voronoi meshes are uniformly accurate.
The problem domain consists of a point heat source in the center of two concentric cylinders with high contrast in thermal conductivity between them. This problem tests the ability of the mesh to accurately capture interior boundaries and temperature diffusion. A version of this benchmark was presented in LaForce et al. (2020); however, in that work, the increase in temperature in the outer cylinder was very small, which made quantitative analysis of the numerical solution difficult. The present work extends the domain radially and increases the power of the heat source significantly so that the dimensionless integral error metric can be used to study the quality of the simulation result. The dimensionless error is defined as where T a is temperature in the analytical solution at monitoring point a and T j is the simulated temperature at the monitoring point. Approximation of the integral in Eq. 1 is most accurate if monitoring point quantities are output at frequent, uniform time intervals, so the solutions are compared every 5 years for a 600-year simulation.

Mesh Generation
The simulation domain is an approximately 1/4 radial volume. The simulation domain is 10 m tall and extends 1000 m in the radial direction, as shown in Fig. 3. The interface between the high-thermal-conductivity inner cylinder and the low-thermalconductivity outer cylinder is 50 m from the origin.

Voronoi Mesh
In the Voronoi mesh, the outer and inner radial boundaries are approximated by nine line segments, each representing 10 • of the 1/4 circle. This was determined to be sufficiently accurate for material boundary and results in fewer grid cells at the material interface and outer boundaries than would be required for finer resolution. As can be seen in Fig. 3, the mesh captures the piecewise linear approximation to the interface between the two material types to a high degree of accuracy. The mesh on interior and exterior input surfaces is always populated by triangular faces, due to the way they are constructed by populating the surfaces before the interior of the volumes. The simulation mesh has 196,546 grid cells, with 4306 cells in the inner material and 192,240 in the outer material. The mesh is generated with no user-defined maximum limit on the size of grid cells. Nevertheless, a large number of cells is required because the meshing software generates isotropic grid cells and for high aspect ratio problems such as the 1000 × 1000 × 10 m 1/4 radial domain, the size of the grid cells in the horizontal directions is limited to a fraction of the vertical thickness. In this mesh the average cell volume is 39.8 m 3 and the largest is 102.0 m 3 .
Eight monitoring points at r = 25 m and eight at r = 62.5 m, evenly spaced every 10 • from 10 • to 80 • , were requested. Four were retained at r = 25 and two were retained at r = 62.5, as shown in Table 1 and as circles in Fig. 3. The purpose of multiple monitoring points at the same radius is to check for grid orientation effects in the 1/4 radial model as well as overall suitability of the mesh for simulation of this problem.

Flexed Hexahedral Mesh
The flexed hexahedral mesh used for comparison is shown in Fig. 4. It was generated using Cubit (Skroch et al. 2021) meshing software. The inner and outer radial boundaries are generated using the cylindrical shape type defined by the software and clipped into a 1/4 radial section. The grid resolution is chosen to be similar to the Voronoi mesh. It is uniformly five grid cells thick in the z-direction and has a total of 142,560 grid cells, with 4,320 cells in the inner material and 138,240 in the outer material. This mesh is considered to be of high quality, with all cells having a scaled Jacobian Red lines on a and c indicate the interface between the two materials and boundary of the model. In subfigures b and d, ID 2 (red) is the low-thermal-conductivity medium and 1 (blue) is the higher-thermal-conductivity medium. Light blue circles indicate monitoring points at radius 25 m from the heat source, and yellow circles indicate monitoring points 62.5 m from the heat source of between 0.780 and 1.0, where 1.0 is a brick with right angles at all corners and any value above 0.65 generally considered acceptable for PFLOTRAN simulations. All 2,115 of the cells with scaled Jacobian less than 0.95 are in the inner material.
As can be seen in Fig. 4, inside material 1 the grid cells are flexed in order to interpolate from the two linear boundaries to the curved boundary. In this area it is expected that simulated solutions will have lower accuracy because the cell boundaries are not aligned perpendicular to the radial flow direction. Outside the material boundary, the hexagonal grid cells are flexed to interpolate from the inner radial boundary to the outer radial boundary. The radial fluxes between these cells are aligned nearly perpendicular to the grid cell faces, and in this area the simulated solutions should have accuracy similar to the Voronoi mesh.  It is not possible to impose monitoring points in the hexahedral mesh without biasing the mesh, so grid cell centers nearest the six monitoring points given in Table 1

Analytical Benchmark
The analytical solution for heat conduction in two concentric cylinders requires that the materials have identical properties except for a large contrast in thermal conductivities. Both materials have porosity 0.28, rock density 2,650 kg/m 3 and specific heat 1,000 J/(kg K). The inner cylinder has thermal conductivity 16.7 W/(m K), representative of steel, while the outer cylinder has thermal conductivity 2.0 W/(m K), representative of a generic shale host rock. The outer material is assumed to be infinite in the radial direction. LaForce et al. (2020) present the one-dimensional radial analytical solution in detail along with a simpler single-phase benchmark.
The initial temperature is 25 • C and the thermal input is a point source at the origin, with constant energy of 4 × 10 3 J/s for 600 years. The domain is water-saturated, with liquid density fixed at 998.6 kg/m 3 , liquid viscosity fixed at 1 × 10 −3 Pa s, and liquid specific heat of 4,186.0 J/(kg K), and the no-gravity option is used in the simulation. This guarantees that fluid properties do not change with increasing temperature, which would cause flow away from the origin and disrupt temperature diffusion. The outer boundary is held at the initial condition to mimic an infinite domain. The x = 0 and y = 0 boundaries are closed, so they are reflective boundaries and the full radial solution is simulated with 1/4 of the domain.

Simulation Results and Error
The simulated temperature distribution on the Voronoi mesh after 200 and 600 years are shown in Fig. 5 and the temperature profiles at the monitoring points are in Fig.  6. As can be seen in Fig. 5, the high-thermal-conductivity material for radius r < 50 m results in rapid propagation of temperature from 88 • C near r = 0 to 65 • C at the material boundary, while in the low-thermal-conductivity region of the outer material, the temperature decreases to 45 • C by r = 100 m and to within 0.1 • C of the initial condition by r = 440 m and within 0.05 • C of the initial condition by r = 500 m . This supports the assumption that this model is large enough to represent a semiinfinite-acting domain. Figure 6 shows that on the Voronoi mesh, the simulated temperature at the monitoring points closely tracks the analytical solution. The dimensionless integral error, E, is calculated for the Voronoi and flexed-hexahedral meshes at the monitoring points using Eq. 1 and the temperature at 5-year intervals. Table 1 shows the errors from the analytical solution at each monitoring point. On the Voronoi mesh, the simulation error is on the order of 1 × 10 −3 [−], with error slightly higher at r = 62.5 m than at r = 25 m. The largest error on the Voronoi mesh is at TD2b, as shown in Table 1. This indicates that the mesh is of very high quality for flow simulations using TPFA, and the similarity of the errors for monitoring points with the same radius shows that the quality of the mesh is uniform across the domain.
On the flexed-hexahedral mesh, the error near r = 25 m is at least an order of magnitude larger than the error at these monitoring points on the Voronoi mesh, as can be seen in Table 1. The largest error is at the monitoring points nearest the domain boundaries TD1a and TD1d. This was anticipated, as the non-orthogonal fluxes inside   Table 1 the inner material bias the flow direction away from radial. For r > 50 m, the mesh is flexed such that the flow is nearly perpendicular to the grid cell faces, and the error is of the same magnitude as the Voronoi mesh, with smaller error than the Voronoi mesh at TD2a and larger error at TD2b.

Test Structures
The four test structures are shown in Fig. 7. Each of the test volume meshes in this section have a similar resolution to the tetrahedral meshes in Gross et al. (2019). In that work the layered model (test structure 1) has a relatively coarse mesh while the pinch-out model (test structure 2) and interior lens model (test structure 3) used a large number of grid cells to stair-step the angled interfaces. The Voronoi meshes of these geometries are generated with a similar spatial resolution by enforcing a maximum grid cell size, with the smallest maximum grid cell corresponding to the largest number of grid cells in the mesh. This wide range of spatial resolutions of the mesh on comparable sized domains serves to test the meshing algorithm and simulator on problems from tens of thousands to over a million grid cells.

Test Structure 1
Test structure 1 is a simple layered system as shown in Fig. 7a. The model is 20 × 20 × 12.0766 m. The layer thickness from top to bottom is 2.9661, 2.8597, 4.1742, and 2.0766 m. The Voronoi mesh is shown in Fig. 8 and has 40,439 elements.  Table 2 shows the input and meshed volumes for the full model and each of the four regions. As can be seen, the relative difference between the meshed and input volumes is small and of similar magnitude for all the regions. The relative difference of the full model is 1.84 ×10 −4 [−], while the difference in each individual material is smaller, between 4 and 5 ×10 −5 [−].

Test Structure 2
Test structure 2 contains three materials with material 2 pinching out in the interior of the model as shown in Fig. 7b. The model is 20 × 20 × 10 m with z in the range (0, −10) m. Material 2 pinches out on the line parallel to the y-axis at (9, y, −4.5081) m and intersects the x = 20 m boundary of the model at (20, y, −1.641) m and (20, y, −7.4393) m. The Voronoi mesh is shown in Fig. 9 and has 1,059,151 elements. Figure 10 shows the details of the surface mesh at the pinch-out. As can be seen, all of the cells have triangular faces on the interior and exterior surfaces defined by the user due to the way VoroCrust constructs the mesh. Small grid cells are needed for the isotropic mesh a short distance away from the intersection where two of the material boundaries are nearby each other. As can be seen in Table 2, the relative difference between the meshed and input volumes is again small and of similar magnitude in each region. The total relative difference is 1.75 ×10 −5 [−], while the difference in each region is slightly lower, with the lowest error in the wedge where the mesh is more refined. Coarser meshes of this structure were considered in LaForce et al. (2021).

Test Structure 3
Test structure 3 contains an interior lens surrounded on all sides by a different material as shown in Fig. 7c. The model is 20 × 20 × 10 m with z in the range (0, −10) m. The lens has a flat bottom at z = −4 m and a top surface with varying height. It pinches to zero thickness on all sides. The Voronoi mesh for test structure 3 is shown in Fig.  11 and has 442,994 elements. Figure 12 shows the details of the mesh on top of the Relative difference is calculated as   As can be seen on the right side of Fig. 12, many of the grid cells adjacent to the boundary are long and thin. One such cell is shown on the left side of Fig. 13. This is because distance between cells along the boundary between the materials is on the order of 0.2 m, but the height is limited by the thinness of the lens near the edge and distance to other material interfaces. Much like test structure 2, very small grid cells are required to mesh the top surface a short distance away from the interior boundary, where height can be as low as 0.01 m. The narrowness of the pointed region on the right side of the figure results in the smallest grid cells in the model, with volume on the order of 5×10 −7 m 3 , such as the one shown on the right side of Fig. 13.
Although the small and narrow grid cells on the boundary of the lens can negatively impact simulation time because of their small size, they will not impact simulation accuracy, as they are Voronoi cells, and fluxes between adjacent cells remain perpendicular to the cell faces. Table 2 shows that the relative difference between the meshed and input volumes is on the order of 1×10 −5 [−] for both regions and the total volume, with the smallest difference in the lens where the mesh is more refined.

Test Structure 4
Test structure 4 contains three material layers which are offset across a fault that passes through the model at a 60 • angle from horizontal as shown in Fig. 7d. The model is     Gross et al. (2019) 20 × 20 × 10 m. The plane dividing material 2 and 5 is parallel to the z-axis at z = −5 m, the plane dividing material 3 and 4 is parallel to the z-axis at z = −2.8096 m, and the plane dividing material 1 and 3 is parallel to the z-axis at z = −7.8686 m. The plane of the fault intersects the bottom of the model on the line (x, 10.7735, −10) m. This structure is similar to the fault example used in Wu and Parashkevov (2009), and large grid orientation effects were observed for vertical flow in that work. The Voronoi mesh for test structure 4 is shown in Fig. 14 and has 228,965 cells. Table 2 shows that the relative difference between the meshed and input volumes is on the order of 1 × 10 −5 [−] for all five regions and the total volume, with slightly lower difference in the smaller material volumes 1 and 4.

Test Case Simulations
In this section, the simulation results for the eight example problems on the four test structure meshes are shown. Simulation parameters are in Table 3. The meshes described in the previous section and shown in Figs. 8,9,11,and 14 Cases 1.1, 2.1, 3.1, 4

.1: Analytical Benchmarks
Test cases 1.1 to 4.1 are single-phase variably saturated flow through a homogeneous porous medium using PFLOTRAN's RICHARDS mode. There is an enforced pressure gradient driving flow upwards in the absence of gravity. The purpose of this benchmark is to check for grid orientation effects caused by the internal material boundaries in each domain. The simulation parameters are shown in Table 3. The simulated pressure profiles are shown on a slice through the interior of the domain in Fig. 15. The vertical flux through the simulation is calculated analytically using the equation In all the cases cross sectional area A = 400 m 2 , and permeability k = 1 × 10 −12 m 2 , while liquid viscosity and density are fixed at μ = 1.0 × 10 −3 Pa s and ρ = 998.6 kg/m 3 , respectively. For test cases 1.1 and 2.1, the enforced pressure gradient is ∇ P = 1 × 10 5 /12.0 Pa/m, which gives flux q = 3.33 kg/s, which is matched to seven significant figures in the test case 1 simulation (a relative error of 3.33 × 10 −8 [−]), and to five significant figures for the test case 2.1 simulation (4.0 × 10 −8 [−] relative error). For test cases 3.1 and 4.1, the enforced pressure gradient is ∇ P = 1 × 10 5 /10.0 Pa/m, which gives flux q = 3.9944 kg/s, which is matched to the precision of the simulation output (9-significant figures) for both test cases 3.1 and 4.1 simulations, indicating a relative error of less than 1.0 × 10 −9 [−]. This result, along with visual inspection of the pressure distribution in Fig. 15 demonstrates that there is little to no grid effect on the simulated flux for this benchmark.

Test Cases 1.2, 2.2, 3.2, 4.2: Two-Phase Flow in Heterogeneous Media
This benchmark simulation tests the capability of the simulation model to accurately capture two-phase flow with large property contrasts between regions using PFLO-TRAN's GENERAL mode. The model has fixed pressure at the bottom of the domain, and constant water saturation, S w = 50%. The infiltration rate of water into the top boundary is 4, 000 kg/year. The initial conditions are 6% water saturation and pressure of 1 × 10 5 Pa. These are modified from Gross et al. (2019) to avoid high suction when the Otowi is at the top boundary of the domain in test case 4.2. The boundary conditions are given in Table 3. A quantitative analysis of the grid-resolution error in test case 4.2 is shown in LaForce et al. (2021).
The permeability and two-phase flow properties for all the models are given by the Santa Fe or the Otowi rock properties in Table 4 and shown in Figs. 8,9,11 and 14. This simulation is isothermal with fluid properties representative of air and water at T = 20 • C.
The simulation results are shown in Figs. 16 and 17. As illustrated by the figures, there are two distinct equilibrium saturation profiles for the high-and low-permeability regions. In the low-permeability and high-capillary pressure Otowi rock type, the liquid saturation at the bottom of the region is above 50% and decreases upwards. The exception to this is the Otowi region at the lower left corner of test case 4.2. This region has lower saturation because of the imposed saturation S w = 50% and pressure P = 1.1 × 10 5 Pa at the base of the model. In the Santa Fe rock type, liquid saturation is on the order of 6% to 7% and relatively uniform. Liquid is able to flow at a low saturation through the Santa Fe as a result of the higher permeability, higher relative permeability to water and lower capillary pressure relative to the Otowi rock type. The pressures in the Santa Fe are on the order of the imposed pressure at the base, while pressures in the Otowi are about 1 × 10 4 Pa lower, Fig. 16 Steady-state liquid pressure for test case 1.2 (a), test case 2.2 (b), test case 3.2 (c) and test case 4.2 (d). Test cases 1.2, 2.2 and 3.2 are shown on the X Z plane on a slice through the plane y = 10 m, while the solution for test case 4.2 is shown on the Y Z plane on a slice through the plane x = 10 m. Red lines indicate boundaries of the input surfaces with a sharp transition between the two. This is due to the different capillary pressure curves in the two media. The low-pressure regions below the Otowi in Fig. 16b and to a lesser extent Fig. 16d are because there is an enforced infiltration rate across the entire top surface and constant saturation and pressure at the base of the model, regardless of rock type. The fluids are flowing around the low-permeability Otowi formation and may be flowing either in or out of the base of the model to maintain equilibrium.

CO 2 Storage in Rock Springs Uplift Geological Model
In this section, the full workflow for generating a mesh and simulating a complex model are demonstrated on a field-scale CO 2 storage project. Simulation of CO 2 storage in the subsurface is a challenging example because of the fingering phenomena that occur when CO 2 is injected into a water-saturated porous medium. These fingers form even when the porous medium is homogenous and the injection uniform and constant (Wang et al. 2013). In real geological porous media, viscous fingering can be exacerbated by rock heterogeneities and uneven top structure driving buoyant flow, and the fingers are typically smaller than the grid resolution of the mesh. As was seen in the Sleipner gas project, CO 2 migration plume migration has been difficult to reproduce using industry-standard reservoir simulators, including PFLOTRAN, the simulator in the present work (Cowton et al. 2018 and references therein.) The present example geologic and conceptual model are based on a potential CO 2 injection site that has been studied extensively by Stauffer et al. (2009), Deng et al. (2012 and Harp et al. (2017). A different part of the uplift is used in the present study, Fig. 17 Steady-state liquid saturation for test case 1.2 (a), test case 2.2 (b), test case 3.2 (c) and test case 4.2 (d). Test cases 1.2, 2.2 and 3.2 are shown on the X Z plane on a slice through the plane y = 10 m, while the solution for test case 4.2 is shown on the Y Z plane on a slice through the plane x = 10 m. Red lines indicate boundaries of the input surfaces along with simplifications to the smaller-scale rock heterogeneity and thermodynamics of the flow simulation to allow for studying meshing and geological structure. Three realizations of the mesh are generated and simulations are run on each.

Meshing
The surfaces for input are assumed to be available as height maps on regularly spaced grids, as these can be generated by geological mapping software. A slice of the Rock Springs uplift in Wyoming from EarthVision (Dynamic Graphics I 2021) commercial software is shown in the top of Fig. 18.
The geological surfaces are trimmed to the area of interest and encapsulated into a single volume with watertight intersections between all surfaces. Three surfaces representing the top Phosphoria (caprock), top Weber sandstone (reservoir) and top Madison limestone (reservoir base) are preserved and clipped to a box approximately 20 × 20 km in the x and y directions, as indicated by the selected area is shown on the bottom of Fig. 18. The simulation domain includes one flank of the uplifted structure, where CO 2 will be injected and the crest where the CO 2 is expected to migrate and pool. The over-and underlying formations are included to absorb the overpressure created by CO 2 injection, and also to allow for the possibility of CO 2 migration into the surrounding strata. A flat bottom surface is created in the crystalline basement for simplicity. An isometric view of the region definitions and surface mesh are shown in Fig. 19. A depth map of the top Weber Sandstone layer in the simulation domain is shown in Fig. 20a. The present simulation is focused on fluid flow in the Weber sandstone, so a very coarse mesh is desirable in the entire Phosphoria and crystalline basement regions. VoroCrust is allowed to mesh these areas as coarsely as possible, as can be seen in Fig. 19. A more refined mesh is desired throughout the Weber Sandstone formation. It is only 100 to 200 m thick so the largest isotropic grid cells will be 50 to 100 m tall, around half the layer thickness, which is sufficient away from where the CO 2 plume is expected to migrate. To determine areas for further refinement, a mesh where VoroCrust is allowed to create the coarsest mesh possible while honoring the input volumes was generated and a simulation of the proposed injection scheme was run to 100 years to see where the plume traveled (not shown). The reservoir volume nearest the wells and containing the plume is the primary area of interest as the vertical thickness of the CO 2 plume in the reservoir is likely to be on the order of a few meters during buoyant migration Dance et al. (2019). The plume was limited to the area within the L shape shown in Fig. 20b, which is defined as the location of refinement. A maximum grid cell size was enforced at the top Weber geological surface in the area where the plume traversed. The level of refinement was based on the x-distance from the injection wells. The area immediately around the wells was refined to a maximum cell radius of 20 m, with the maximum grid cell size gradually increasing to 50 m along the top structure. Maximum grid cell radius was increased in small increments as sharp changes in mesh refinement were observed to bias simulation results. The volume of the cells at the top of the Weber with refinement in the path of the plume is shown in Fig. 20b. A striping pattern can be seen where the mesh is coarsened up-and down-dip of the injection wells and there is also an area of coarser cells near the large y boundary of the model. The refined mesh along the top Weber with coarsening above and below can be seen in the y-slice shown in Fig. 21. The mesh is finest at the top Weber and coarsens both downward and upward into the other formations, with very large cells achieved in the crystalline basement rock.
Refined regions were also defined near each of the injection wells to a depth of 50 to 55 m below the well completion. Two of the injection wells are shown as red circles in Fig. 21. It can be seen that the refined region extends further down into the Weber Sandstone at and between the wells, to allow for more accurate modelling of the CO 2 plume in the near-well region and as it travels vertically up to the caprock.
Three realizations of the target mesh were generated to study the impact of the mesh on the CO 2 flow field. The details of each mesh are shown in Table 5. The injection well locations have y-coordinate 50 m greater in mesh 2 due to the meshing algorithm not retaining the desired injection points in the final mesh. Given the large scale of the simulation domain, this is not anticipated to have a significant impact on simulated results.

Simulation Implementation
The permeability and porosity of the regions are shown in Table 6 and are the mean of the high case from Deng et al. (2012). Principal permeability directions are aligned with coordinate axes. Permeability at cell faces not aligned with coordinate axes are calculated using appropriately averaged vertical and horizontal values, as full-tensor permeability is not currently available in PFLOTRAN. The vertical to horizontal permeability ratio of the Weber and crystalline basement are 0.1, while for the Phosphoria the ratio is 0.001, to prevent vertical intrusion of CO 2 into the overburden. Within each region properties are considered to be homogeneous. High permeabilities were chosen so that only four injection wells at 1 km spacing were required to inject CO 2 at the desired rate of 10 kg/s into each well, for a total of 1.26 million metric tons per year for 100 years. Water is also injected in each well at a rate of 1 × 10 −5 kg/s to prevent drying out. Explicit well models are not available in PFLOTRAN, so wells Injector 2 (x, y, z) ( 1 8 ,500, 9,500, −1,390) (18,500, 9,550, −1,390) (18,500, 9,500, −1,390) Injector 3 (x, y, z) ( 1 9 ,500, 8  are modeled as point sources in the reservoir at the center of the cells shown in Table  5. Using point sources as wells will negatively impact accuracy of the pressure and CO 2 distribution at the wells, but will not significantly impact the simulation away from the wells. The simulations are isothermal. These simplifications from the cases considered in Stauffer et al. (2009), Deng et al. (2012 and Harp et al. (2017) allow for study of the interaction of meshing and the geological structure on the simulation in isolation from other sources of complexity. The minimum and maximum x boundaries are open to flow and assigned initial reservoir properties. The minimum and maximum y boundaries are closed to represent sealing faults on either side of the structure. These fault seals are not present in the original geological model and are introduced so that CO 2 will pool at the top of the structure instead of migrating out of the model in the positive y-direction. The van Genuchten capillary pressure curve and the Mualem-van Genuchten relative permeability curves with parameters in Table 7 are used (Lichtner et al. 2020) with the residual saturations taken from Deng et al. (2012).
Because of the complexity of simulating two-phase flow on highly unstructured Voronoi meshes, new experimental solvers were also used in the simulation. As discussed in Bonduà et al. (2017) andFreeman et al. (2014) the Jacobian matrices for unstructured meshes are more difficult to solve than for structured meshes both because α is the inverse of the air entry pressures and M is the empirical van Genuchten exponential parameter of additional connections and the loss of the banded matrix structure. As shown in Table 5, the average number of connections for each cell is 11 for all three Rock Springs Uplift meshes and the maximum is 36 or 37 connections. For comparison, interior cells in a hexahedral mesh have six connections.
A new type of nonlinear solver called Newton trust-region dogleg Cauchy (NTRDC) has been implemented into PETSc, which is part of PFLOTRAN's parallel framework Park et al. (2021). The solver is an official capability for PETSc as of version 3.17.1. NTRDC is more capable of resolving complex nonlinear constitutive models than traditional Newton-Raphson solvers. The flexible generalized minimum residual method (FGMRES) is used as the linear solver, combined with a two-stage preconditioner, the constrained pressure residual (CPR) preconditioner using the alternate-block factorization (ABF) decoupler. The combination of FGMRES, CPR, and ABF seems to be effective for simulation domains with large contrast in permeability and unstructured grids. The combination of these preconditioners, linear, and nonlinear solvers has been effective in solving very complex constitutive models, heterogeneous domains, and unstructured grids (Park et al. 2021).
With the combination of the advanced linear and nonlinear solvers and capillary pressure and relative permeability smoothing it was possible to complete the 100-year simulations in 24 to 28 h using 144 cores of a parallel super-computer, as shown in Table 5.

Rock Springs Uplift Simulation Results
The liquid saturation at the top of the Weber Sandstone is shown in Figs. 22 and 23. In each of the simulations the CO 2 travels updip in the negative x-direction to the top of the anticline structure at early times (see Fig. 20a). Figure 24 shows that the CO 2 plume is one grid cell thick as it follows the caprock updip away from the injection wells. The so-called leopard spot pattern is caused by the mesh not fully resolving the CO 2 plume in the vertical direction, even with the small grid cells at the well. The plotted saturation for each grid cell is representative of the saturation averaged over the cell volume, where some fraction of each cell is in the plume and the rest is in the water beneath it. Every polyhedral volume is of a different size and each center is a different distance from the top of the Weber Sandstone, resulting in disparate saturations for neighboring cells. After 10 years, the CO 2 has begun to pool in a local high point in the structure, and by 20 years it is visibly migrating in the positive y-direction towards the highest point in the anticline, which is at the large y boundary of the model. Migration in the y-direction is slower and much more uniform because the incline is very small compared to the flank of the structure, as can be seen in Fig. 20a. By 50 years, the CO 2 is pooling against the sealing fault at the top of the structure, and the structure fills downward until 100 years, at the end of the simulation. It is unclear what causes the narrowing of the CO 2 plume just below the boundary of the structure (top subfigures g to l of Fig. 23). It may be an artifact of the grid coarsening at the mesh boundary.
Buoyancy-driven flow of CO 2 is a viscous-unstable displacement and forms fingers, even in homogeneous media (Wang et al. 2013). Using structured meshes these fingers are often suppressed because the uniformity of the grid tends to smooth these fronts out, which results in an under-prediction of plume migration distance (Cowton et al. 2018). Because of the random orientation of the Voronoi mesh in this work, the opposite is observed, and gas fingering is triggered by irregularities in the mesh from the first year of injection, as shown in the top of Fig. 22.
The individual fingers triggered by the randomness of the grids are not physically correct, as the mesh is not fine enough to fully resolve them. Comparison of simulated and experimental gas fingering is the topic of future work, but the meshes with randomly generated elements inside the input volumes provide a simple and rapid way of generating an ensemble of models that might eventually be used to study uncertainty created by viscous fingering in displacements. In particular, a suite of equally probable models such as these could be used to create confidence intervals on how likely CO 2 is to appear at a monitoring well updip of the injectors, or the range of breakthrough times and ultimate thickness of the plume at points on top of the anticline.

Conclusions and Future Work
This paper has presented a series of examples of the generation and use of highly unstructured VoroCrust Voronoi meshes for subsurface fluid flow using the simulator PFLOTRAN, which uses TPFA and is known to yield results that are sensitive to poor-quality meshes. A two-material benchmark problem demonstrates the ability of the mesh to accurately capture the material interface and the simulator to create results on the mesh that are more accurate than a similar-resolution flexed hexahedral mesh. Four single-phase and four two-phase benchmarks on meshes of four common reservoir structures (layering, pinch-out, interior lens, and offset fault) are also shown.
Finally, three realizations of the mesh for a mapped geological structure for CO 2 storage in the subsurface are presented. This example demonstrates the capability of the presented workflow to mesh and simulate a challenging two-phase flow problem on subsurface features with a realistic scale and level of complexity. Furthermore, the ability to easily create three realizations of the mesh with the same resolution allows for exploration of the impact of the mesh on the simulation result.
The long-term goal of the present study is to develop a meshing, simulation, uncertainty assessment, and visualization software package that is open-source, flexible, and well-suited for large-scale fluid flow simulations in complex geological models as part of the Sandia National Laboratories Geologic Disposal Safety Assessment (GDSA) Framework (https://pa.sandia.gov/). The present work addressed the meshing, simulation and visualization pieces of the proposed workflow.
Currently, LaGriT https://github.com/lanl/LaGriT, PFLOTRAN https://www. pflotran.org/, and ParaView https://www.paraview.org are open-source software. The scripts required to use VoroCrust meshes in PFLOTRAN simulations and visualize the results are available in the PFLOTRAN source code repository in the folder pflotran/src/python/vorocrust. It is a near-term goal of this project to create an open-source version of VoroCrust so that the software for the entire workflow will be freely available.
Other areas of active investigation are to run an SA/UQ workflow including uncertainty in geological structure as a sampled parameter. It is also necessary in the future to assess whether the fingers triggered by the mesh are a realistic representation of viscous fingers that are expected to form in displacements such as CO 2 injection for storage. Future work is also planned to make it possible to mesh faults that terminate within the model.