Daniel W. Schmid

MILAMIN: MATLAB-based finite element method solver for large problems

The finite element method (FEM) combined with unstructured meshes forms an elegant and versatile approach capable of dealing with the complexities of problems in Earth science. Practical applications often require high-resolution models that necessitate advanced computational strategies. We therefore developed “Million a Minute” (MILAMIN), an efficient MATLAB implementation of FEM that is capable of setting up, solving, and postprocessing two-dimensional problems with one million unknowns in one minute on a modern desktop computer. MILAMIN allows the user to achieve numerical resolutions that are necessary to resolve the heterogeneous nature of geological materials. In this paper we provide the technical knowledge required to develop such models without the need to buy a commercial FEM package, programming compiler-language code, or hiring a computer specialist. It has been our special aim that all the components of MILAMIN perform efficiently, individually and as a package. While some of the components rely on readily available routines, we develop others from scratch and make sure that all of them work together efficiently. One of the main technical focuses of this paper is the optimization of the global matrix computations. The performance bottlenecks of the standard FEM algorithm are analyzed. An alternative approach is developed that sustains high performance for any system size. Applied optimizations eliminate Basic Linear Algebra Subprograms (BLAS) drawbacks when multiplying small matrices, reduce operation count and memory requirements when dealing with symmetric matrices, and increase data transfer efficiency by maximizing cache reuse. Applying loop interchange allows us to use BLAS on large matrices. In order to avoid unnecessary data transfers between RAM and CPU cache we introduce loop blocking. The optimization techniques are useful in many areas as demonstrated with our MILAMIN applications for thermal and incompressible flow (Stokes) problems. We use these to provide performance comparisons to other open source as well as commercial packages and find that MILAMIN is among the best performing solutions, in terms of both speed and memory usage. The corresponding MATLAB source code for the entire MILAMIN, including input generation, FEM solver, and postprocessing, is available from the authors (http://www.milamin.org) and can be downloaded as auxiliary material.

Enhanced mass transfer through short-circuit diffusion: Growth of garnet reaction rims at eclogite facies conditions

Abstract In the Monte Rosa area, Northern Italy, the assemblage garnet + phengite + quartz was formed in polymetamorphic metapelites from pre-existing biotite and plagioclase during an eclogite-facies Alpine metamorphic overprint. While phengite nucleated within plagioclase, garnet formed 10 to 20 μm wide continuous rims along the original biotite-plagioclase interphase boundaries. Garnet formation involved diffusion of Ca-, Fe-, and Mg-bearing species across the growing rims. The garnet rims show an overall asymmetric compositional zoning and distinct nanometer scale chemical patterns across grain boundaries within the garnet polycrystal. The compositional patterns resulted from the interplay of diffusion along grain and interphase boundaries and volume diffusion. Individual garnet crystals are separated by low-angle grain boundaries, which are predominantly oriented perpendicular to the reaction rim giving rise to an overall palisade microstructure. Grain boundaries contain arrays of closely spaced channels about 2 nm wide, which are filled with an amorphous phase. These channels represent direct links between the garnet-biotite and garnet-plagioclase reaction fronts and provide pathways for fast diffusion. From numerical simulation of the observed chemical patterns, diffusion of Ca and Fe is inferred to have been 5.6 orders of magnitude faster within grain boundaries than through the interior of garnet grains. Although the amorphous grain boundary phase only takes up a small fraction (2.5 per mil) of the total volume of the garnet polycrystal, short-circuit diffusion along grain boundaries contributed substantially to material transfer across the growing garnet rim. Despite the presence of fast pathways, bulk diffusion was too slow to allow for chemical equilibration of the phases involved in garnet rim formation, even on a micrometer scale. Relying on published garnet volume diffusion data, Dgb (where gb = grain boundary), is estimated to be on the order of 10-(19.20) m2/s at the estimated reaction conditions of 650 °C and 12.5 kbar, where DCagb is slower than DFegb by a factor of ten.

Folding in power-law viscous multi-layers

We study high-amplitude folding in layered rocks with two-dimensional numerical simulations. We employ the finite-element method to model shortening of an incompressible multi-layer with power-law viscous rheology. The Lagrangian numerical mesh is deformed and re-meshed to accurately follow the layer interfaces. Three settings are considered: (i) pure shearing of a confined multi-layer, (ii) simple shearing of a multi-layer above a detachment, and (iii) slump folding owing to gravity sliding. In our pure shear simulations, finite-amplitude folds always develop despite confinement and thin weak interlayers. The fold shapes can be significantly irregular, resulting from initial geometrical heterogeneities that are perturbations of the layer interfaces and differences in layer thickness. The bulk normal viscosity of the multi-layer decreases significantly with progressive folding. This structural softening decreases the bulk normal viscosities by a factor of 2–20. For simple shear, the multi-layer does not develop asymmetric fold shapes significantly. Fold axial planes in the multi-layer are mostly curved and not parallel. For slump folding, fold shapes can be significantly asymmetric exhibiting strongly curved fold axial planes and overturned fold limbs. The rheology of the competent layers has a major impact on the fold shapes for gravity-driven multi-layer folding.

Fold amplification rates and dominant wavelength selection in multilayer stacks

A combination of thin- and thick-plate theories, and finite element models is used to systematically analyze folding in multilayer stacks. We show that if the interlayer spacing is large, individual layers fold as single layers, if the spacing is small the entire stack folds as one effective single layer. In between, a third folding mode exists that is characterised by a dominant wavelength that scales with n 1/3, irrespective of total number of layers, n. The maximum growth rates in the true multilayer-folding mode are higher than the corresponding single layer growth rates, increase with n and are bounded by a saturation value that is directly proportional to the viscosity contrast. This growth rate saturation as well as the applicability of the true multilayer-folding mode with respect to interlayer spacing can be explained by the normal and inverse contact strain theory. The true multilayer-folding mode is expected to be the most frequent mode in nature, because it exhibits the highest growth rates and has a relatively large applicability range with respect to interlayer spacing. The increased growth rates in multilayer folding are especially important for systems where the corresponding single layer values are not sufficient to drive the folding instability, such as folding in low-viscosity contrast layers and detachment folding.

Automated thermotectonostratigraphic basin reconstruction: Viking Graben case study

We present a generic algorithm for automating sedimentary basin reconstruction. Automation is achieved through the coupling of a two-dimensional thermotectonostratigraphic forward model to an inverse scheme that updates the model parameters until the input stratigraphy is fitted to a desired accuracy. The forward model solves for lithospheric thinning, flexural isostasy, sediment deposition, and transient heat flow. The inverse model updates the crustal- and mantle-thinning factors and paleowater depth. Both models combined allow for automated forward modeling of the structural and thermal evolution of extensional sedimentary basins. The potential and robustness of this method is demonstrated through a reconstruction case study of the northern Viking Graben in the North Sea. This reconstruction fits present stratigraphy, borehole temperatures, vitrinite reflectance data, and paleowater depth. The predictive power of the model is illustrated through the successful identification of possible targets along the transect, where the principal source rocks are in the oil and gas windows. These locations coincide well with known oil and gas occurrences. The key benefits of the presented algorithm are as follows: (1) only standard input data are required, (2) crustal- and mantle-thinning factors and paleowater depth are automatically computed, and (3) sedimentary basin reconstruction is greatly facilitated and can thus be more easily integrated into basin analysis and exploration risk assessment.

Brittle fracture during folding of rocks: A finite element study

The goal of the present work is the development of a novel computational analysis tool to elaborate folding-induced fracture of geological structures. Discrete failure of brittle rocks is characterised by three sets of governing equations: the bulk problem, the interface problem and the crack problem. The former two sets which define the deformation field are highly nonlinear and strongly coupled. They are solved iteratively within a Hansbo-type finite element setting. The latter set defines the crack kinematics. It is linear and solved in a single post-processing step. To elaborate the features of the computational algorithm, we define a unique benchmark problem of a single, geometrically nonlinear plate, which is subjected to layer-parallel in-plane compression combined with different levels of superposed in-plane shear. The resulting folding, or buckling, induces brittle failure in the tensile regime. By systematically increasing the shear strain at constant compression, we develop crack deviation angle versus shear-to-compression ratio tables. We determine the corresponding damage zones, analyse the folding modes and elaborate the force versus amplification diagrams. The proposed two-field folding-induced fracture algorithm can ultimately be applied to interpret natural folded rocks and understand their evolution, structural development and histology.

Basin modelling of a transform margin setting: structural, thermal and hydrocarbon evolution of the Tano Basin, Ghana

ABSTRACT This study explores the structural and thermal evolution of the Ghana transform margin. The main objective is to explore how the opening of the Atlantic Ocean and subsequent interaction with the Mid-Atlantic Ridge (MAR) has affected the margins structural and thermal evolution. Two representative evolution scenarios are described: a reference case that neglects the influence of continental breakup and a second scenario that accounts for a possible heat influx during the passage of the MAR as well as magmatic underplating. These two scenarios have further been analysed for the implications for the hydrocarbon potential of the region. The scenario analysis builds on a suite of 2D realizations performed with TECMOD2D, modelling software for automated basin reconstructions. As the observed stratigraphy is input, the structural and thermal evolution of the basin is automatically reconstructed. This is achieved through the coupling of a lithosphere scale forward model with an inverse algorithm for model parameter optimization. We find that lateral heat transport from the passing MAR in combination with flexure of the lithosphere can explain the observed uplift of the margin. These results were obtained for a broken plate elasticity solution with a relative large value for the effective elastic thickness (Te=15) and necking level (15 km). Lateral heat flow from oceanic lithosphere is clearly visible in elevated basement heat flow values up to 50 km away from the ocean–continent transition (OCT). This influx of heat does not seem to have affected the maturation history along the margin significantly. Only the deepest sediments close to the OCT show slightly elevated vitrinite reflectance in simulations that account for the passage of the MAR. In conclusion, it appears that that lateral heat transport from the oceanic lithosphere is instrumental in shaping the Ghana transform margin but seems to have only limited control on the maturation history.

Ray-based seismic modeling of geologic models: Understanding and analyzing seismic images efficiently

AbstractOften, interpreters only have access to seismic sections and, at times, well data, when making an interpretation of structures and depositional features in the subsurface. The validity of the final interpretation is based on how well the seismic data are able to reproduce the actual geology, and seismic modeling can help constrain that. Ideally, modeling should create complete seismograms, which is often best achieved by finite-difference modeling with postprocessing to produce synthetic seismic sections for comparison purposes. Such extensive modeling is, however, not routinely affordable. A far more efficient option, using the simpler 1D convolution model with reflectivity logs extracted along verticals in velocity models, generates poor modeling results when lateral velocity variations are expected. A third and intermediate option is to use the various ray-based approaches available, which are efficient and flexible. However, standard ray methods, such as the normal-incidence point for unmigrat...

On the Stokes-Brinkman Equations for Modeling Flow in Carbonate Reservoirs

Cavities and fractures can significantly affect the flow paths of carbonate reservoirs and should be accurately accounted for during flow simulation. Herein, our goal is to compute the effective permeability of rock samples based on high-resolution 3D CT-scans containing millions of voxels. Hence, we need a flow model that properly accounts for the effects of Darcy flow in the porous material and Stokes flow in the void volumes on relevant scales. The presence of different length scales and large contrasts in the petrophysical parameters lead to highly ill-conditioned linear systems that make such a flow model very difficult to solve, even on large-scale parallel computers. To identify simplifications that render the problem computationally tractable, we analyze the relative importance of the Stokes and Darcy terms for a wide variety of parameter ranges on an idealized 2D model. We find that a system with a through-going free flow region surrounded by a low permeable matrix can be accurately modeled by ignoring the Darcy matrix and simulating only the Stokes flow. Using this insight, we are able to compute the effective permeability of a specific model from a CT-scan that contains more than eight million voxels.

Rigid polygons in shear

Abstract Clasts, inclusions and intrusions in shear are potential recorders of strain, stress, rheology and metamorphism. In order to extract the recorded information, it is essential to have analytical and numerical theories that describe the deformation mechanics of such bodies. To overcome the simplifications of the commonly employed ellipsoid-based shape approximation, a combination of Muskhelishvili-type analytical solutions and finite-element method calculations is used to study the behaviour of (quasi) rigid polygons in shear. The results confirm that the polygon rotation and the pressure perturbation outside rigid polygonal clasts are well approximated by ellipse-based theories. However, this observation does not hold for the inside of these polygons, which show strongly varying values of pressure perturbation and maximum shear stress. For example, pressure perturbations inside the polygons are usually the opposite of the neighbouring matrix values across the polygon-matrix interface. This complex behaviour is summarized in the ellipse decomposition rule that allows for the qualitative understanding of the pressure perturbation in and around a wide range of polygons in shear. Other quantities studied include maximum values of over-pressure relative to the shortening stress, and the area that undergoes overpressure with respect to the clast size. The results demonstrate that overpressure can be twice as large as the rock strength.