How the mesh controls accuracy in geomechanical-numerical models

Given that in-situ data of the stress field are limited, geomechanical models are used to receive a continuous prediction of the 3-D stress in the subsurface. For the numerical solution of these models the Finite Element (FE) method is often used as it allows to discretize complex structures using unstructured meshes. Therefore, the results strongly depend on the FE-mesh resolution, the element types, and the element order used. In a comprehensive approach to optimizing a mesh, it is necessary to find a balance between a detailed high-resolution mesh and the resulting computing time or available computing capacity.

To investigate such a question an extensive model series that change the model geometry and its resolution is required. However, such test series can only be conducted using simplified models, as the effort involved in producing the FE-mesh would otherwise be too substantial. For this reason, the complexity of a real 3-D structure was reduced to a 2-D profile section. As a template for this approach, we use a geomechanical-numerical model of the siting region of Nördlich Lägern for a deep geological repository for radioactive waste in northern Switzerland. Geologically, it consists of the crystalline basement, south-dipping Mesozoic units, and a cover of Cenozoic deposits. The key purpose of our study is to investigate the impact of the FE-mesh on the predicted 3‑D stress field within the thin stratigraphic units of the Mesozoic. The mesh resolution, element type, element order, and solver-specific elements with reduced integration points are tested. All models are calibrated separately with the same set of in-situ stress magnitude data from a borehole to find the best-fit displacement boundary conditions. All results are also displayed along the same borehole trajectory, which is located exactly on the profile section, in comparison to the available in-situ stress magnitude data. The result shows that horizontally elongated hexahedrons are more suitable for thin layers in comparison to tetrahedron elements; higher order elements also offer little added value in such static case.

For our study of Nördlich Lägern we use three different model geometry realizations that resulted from the different stages of the exploration process during the past 15 years. To achieve better comparability, the mechanical properties were also harmonized as far as possible. If almost identical rock properties are used, only small differences in the predicted stress field are visible, mainly in the areas where the stratigraphic boundaries differ between the models. Differences become more significant when the original and deviating rock properties are used. This indicates that the rock properties have a large influence on the model estimates. However, in comparison to predicted bandwidth of the predicted stress field that is controlled by the probability distribution of the rock stiffness in each lithology, the changes to the different model realization are small.