High Performance Computing in Basin Modeling

Abstract

The evolution of sedimentary basins is a complex set of processes developing over millions of years. Often, tectonic activity will create accommodation space that is filled in with sediments. These sediments undergo changes as pressure and geothermal energy are applied to them. Typically, sedimentary basins will be in a marine or fluvial environment, thus the sediments are saturated with a pore fluid, water or brine. Such a depositional environment, treated as a saturated porous medium, will undergo loading and compaction, as more sediments eject pore fluid out of the pore spaces around the solid matrix. The geochemistry of these sediment layers, strata, is of interest for possible hydrocarbon resource production. <br /> Using mathematical models to approximate the evolution of physical states in sedimentary basins, called basin modeling, has progressed in the last decades due to the increasing computational resources available. The computational approach to approximately solving the model equations, formulated as partial differential equations, relies on reducing the physical domain to a "cloud of points", known as a mesh or grid. Physical states are then only defined at these points (or objects derived from them, e.g., faces). The quality of the grid's representation of physical geometry is essential to the quality of the approximate solutions obtained computationally. <br /> In geological basins, the layering of the strata can be very complex. Portions of the domain might receive more or less sediments of varying physical characteristics, such as permeability (the ability to allow fluid to pass through) or compressibility (the change in porosity with respect to a change in stress). The layered sediments will change their geometry as they compact during deposition. Other examples of complexity include sharply discontinuous material properties, e.g., one layer might have drastic differences in permeability compared to its neighbors. Resolving where this change in permeability occurs would then be required to produce accurate simulations. Even more pronounced deformation comes from various tectonic processes that might have caused folding or even fracturing of the rock layers, as is commonly revealed through seismic imaging. In this thesis, we will not consider tectonic processes or fracturing, though these are important features in many geological basins and prime examples of geometrical complexity. <br /> One of the main challenges in accurately simulating physical processes in geological contexts is the ability of the discrete mesh to conform to the physical states and material properties, given that the geometry of the layers can be deformed in ways that make it challenging to line up on a grid. As basin modeling is concerned with the evolution of sediment layers, a time-dependent context, one approach to this problem is to frequently recreate or modify the grid. The physical changes in the domain are adapted to a correspondingly updated grid. This has the advantage that the phyical quantities and material properties will be conforming to the new grid for the solution at the next time step. Recreating the mesh can be an expensive operation (in computational work), and can increase the amount of communication time between processors, when dealing with distributed parallel computing. <br /> As one of the goals of this project has been to consider model applications in the context of High-Performance Computing (HPC), the approach of frequently modifying the grid presents a possible impediment to scaling to higher processor counts. The purpose of this doctoral thesis is to consider techniques to represent the complex geometries found in geological basins in a way that does not require re-meshing. The idea is to implicitly represent the layered structure of a basin, beginning with a fixed global domain, through which sediments are deposited at a certain speed. The interfaces between layers are implicitly given, allowing the identification at a point in the domain of what rock type is there, and so what material properties are appropriate. A well-known technique to evolve an interface, encoded as a fixed value level set of an auxiliary function, is the level set method. The interface moves according to how a given "speed function'' evolves the auxiliary potential function. Whether the potential function at a point is above or below the value (often taken to be 0), determines whether the point is "in or out'', and thus we can reconstruct the identity of the material, what layer it is in, at that point. <br /> The novel contribution of this work is to apply the numerical technique of the level set method for interfaces coupled to physical processes in the modeling of sedimentary basins. First, using the principle of vertical effective stress, the work considers in 1-D the fundamental coupling between (1) the layer interfaces equipped with a speed function and a (2) pore pressure, combining to influence the movement of a rock layer based upon the sedimentation rate and the rate of compaction. As there are no publicly available benchmarks or standard solutions by which to verify a basin simulator, the developed model has been verified against one of the only available analytical solutions relevant to sedimentary basins, formulated as porous media, the Gibson consolidation problem. Further examples of multi-layer setups are presented. <br /> In the second part, we consider a finite element discretization of a pore fluid flow solver coupled to a representation of geologic interfaces through a level set method. A description of the model equations, their discretization and their implementation in the finite element library, deal.II, is given. Subsequently, the verification example is considered, as well as some elaborations to more complex demonstration cases in 2-D, including a basin with a temperature field. A scaling study has been performed on Forschungszentrum Juelich's supercomputer JSC/Juwels to profile the code's parallel performance with distributed memory. <br /> In the third part, three elaborations are briefly considered: (1) a small 3-D example derived from the Gibson consolidation problem, (2) the possibility of incorporating an effect similar to buoyancy in the level set speed function, and (3) the foundations required for adaptive mesh refinement. <br /> A conclusion and some possibilities for extension complete the thesis.