The ongoing digitalisation in geosciences, along with higher resolution data and more powerful computers, necessitates efficient workflows and tools for processing large amounts of data, modelling processes, and addressing new problems and challenges.

The Julia language offers an great basis for this by combining the benefits of a high-level language, such as ease of use and interactivity, with the features of a low-level language, including speed, efficiency, scalability, and native GPU support.

We will present recent efforts to accelerate geocomputing using Julia at scale. These efforts are supported by the PASC-funded GPU4GEO project, in collaboration with the Swiss National Supercomputing Centre. We will showcase the current HPC building blocks that we have developed, which enables geoscientists to write high-performance stencil codes that can scale from their laptops to the largest supercomputers. Additionally, we will demonstrate preliminary results on using automatic differentiation to perform inverse modelling. This allows us to efficiently constrain models with data.

Additionally, we discuss how these recent developments, along with future ones, will accelerate geocomputing and enable the education of the next generation of geoscientists.

How to cite: Omlin, S., Räss, L., and Utkin, I.: Accelerating Geocomputing with Julia at Scale, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-19607, https://doi.org/10.5194/egusphere-egu24-19607, 2024.

Pegmatites (and rare metal granites) are igneous rocks with a granitic composition, characterised by crystal growth dominated texture. They are often enriched in rare elements (such as Li, Cs, Be, Nb, Ta…) and offer valuable deposit of economic interest that belong to the list of critical raw material defined by the European commission. The aim is to use modeling to understand the formation of pegmatite, and especially the parameters that control the migration distance between their sources (granite, migmatites) and their level of emplacement.

We use a two phase flow finite difference code, in Julia, based on the porosity waves with compressible fluid in compressible medium, where the porosity is interpreted as melt [1,2], to model the magma migration inside migmatitics domes. To improve the yet existing codes, we implement temperature in our two phase flow formulation, it will be calculated from energy conservation to take into account the latent heat of the partial melting and the crystallisation. It will allow the to stop the experiments upon magma crystallisation. It will also be use to increase the viscosity of the melt while it cools down.

We here present the results of a preliminary 1D version of our model (without temperature). It shows that an increase in the ratio of matrix permeability over the fluid viscosity results in a greater distance travelled by the melt, for a constant number of time step. For a fluid viscosity of 10⁴ Pa.s the increase of matrix (dynamic) permeability from 10^-13 to 10^-11 m^-2 fasten the migration of the melt of a factor 2 and the travelled distance by a factor 1. When the energy will be implemented, we will compare the impact of matrix permeability, the fluid viscosity and the geothermal gradient on the height of the magmatic migration.

References:

[1] L. Räss, T. Duretz & Y.Y. Podladchikov (2020). https://doi.org/10.1093/gji/ggz239

[2] A. Plunder & al. (2022). https://doi.org/10.1016/j.lithos.2022.106652

How to cite: Bergogne, M., Le Pourhiet, L., Räss, L., Podladchikov, Y., Utkin, I., and Plunder, A.: 3D Modelling of pegmatites migration at the onset of partial melting, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-19212, https://doi.org/10.5194/egusphere-egu24-19212, 2024.

In this contribution we introduce ECOMAN 2.0, an open-source software package for (1) modelling the strain-/stress-induced rock fabrics and related mechanical anisotropy, and (2) performing isotropic and anisotropic inversions using real/synthetic P- and S-wave travel-times and S-wave splitting parameters.  

The strain-induced intrinsic mantle fabrics are modelled inputting the velocity, pressure, temperature and dominant creep mechanism fields from large-scale mantle flow simulations into D-Rex (Kaminski et al., 2004). This open-source software has been parallelized using a hybrid MPI and OpenMP scheme and modified to account for combined diffusion-dislocation creep mechanisms, LPO of transition zone and lower mantle polycrystalline aggregates (Wadsleyite, Bridgmanite, post-Perovskite), P-T dependence of single crystal elastic tensors, advection and non-steady-state deformation of crystal aggregates in 2D/3D cartesian/spherical grids (Faccenda, 2014; Faccenda and Capitanio, 2013). The new version of D-Rex can solve for the LPO evolution of 100.000s polycrystalline aggregates of the whole mantle in a few hours, outputting the full elastic tensor of poly-crystalline aggregates as a function of each single crystal orientation, volume fraction and elastic moduli scaled by the local P-T conditions.

Extrinsic elastic anisotropy due to grain- or rock-scale fabrics or fluid-filled cracks can also be estimated with the Differential Effective Medium (DEM) (Faccenda et al., 2019). Similarly, extrinsic viscous anisotropy can be modelled yielding viscous tensors to be used in large-scale mantle flow simulations (de Montserrat et al., 2021). 

The elastic tensors can then be interpolated in a tomographic grid for (i) visual inspection of the mantle elastic properties (such as Vp and Vs isotropic anomalies; radial, azimuthal, Vp and Vs anisotropies), (ii) generating input files for large-scale synthetic waveform modelling (e.g., SPECFEM3D format), or (iii) P- and S-wave isotropic and anisotropic inversions (e.g., Faccenda and VanderBeek, 2023). The latter can be performed with the new PSI (Platform for Seismic Imaging) module, which includes recently developed techniques for seismic anisotropic inversions of body waves (VanderBeek and Faccenda, 2021; VanderBeek et al., 2023).

How to cite: Faccenda, M., VanderBeek, B., de Montserrat, A., and Yang, J.: ECOMAN 2.0: an open-source software package for coupling geodynamic and seismological modelling., EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-3547, https://doi.org/10.5194/egusphere-egu24-3547, 2024.

Simulating the Earth’s mantle convection at full convective vigor on planetary scales is a fundamental challenge in Geodynamics even for state of the art high-performance computing (HPC) systems. Realistic Earth mantle convection simulations can contribute a decisive link between uncertain input parameters, such as the mantle viscosity structure, and testable preconditions, such as dynamic topography. The vertical deflections predicted by such models may then be tested against history of dynamic topography from stratigraphic observations. Considering realistic Earth like Rayleigh numbers (∼ 10⁸ ) a resolution of the thermal boundary layer of 10 − 50 km is necessary considering the volume of the Earth’s mantle.

Simulating Earth’s mantle convection at this level of resolution requires solving sparse indefinite systems with more than 10¹² degrees of freedom, computationally feasible only on exascale HPC systems. This is achievable only by mantle convection codes providing high degrees of parallelism and scalability. Earlier approaches with prototype frameworks using hierarchical hybrid grids (HHG) as solvers for such systems demonstrated the scalability of the underlying concept for future generations of exascale computing systems.

Building up on the TerraNeo project here we report on the progress of utilizing the improved framework HyTeG (Hybrid Tetrahedral Grids) based on matrix-free multigrid solvers in combination with highly efficient parallelisation and scalability. This will allow to solve systems with more than a trillion degrees of freedom on present and future generations of exascale computing systems. We also report on the advances in developing the scalable mantle convection code TerraNeo using the HyTeG framework to realise extreme-scale mantle convection simulations with realistic, Earth like parametrisation and a resolution in the order of ∼ 1km.

How to cite: D'Ascoli, E., Burkhart, A., Kohl, N., Bunge, H.-P., and Mohr, M.: TerraNeo: Development of a scalable mantle convection code for exascale computing, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-6508, https://doi.org/10.5194/egusphere-egu24-6508, 2024.

How to cite: Duretz, T. and Le Pourhiet, L.: One dimensional models of frictional plastic strain localisation, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-8080, https://doi.org/10.5194/egusphere-egu24-8080, 2024.

Geodynamic forward models are typically computationally expensive, they generally draw on a large set of input and testing data and will normally have quite large model domains in order to ensure that the results are practical and applicable. To reduce simulation time many programs will use some form of parallel computing in their calculations, e.g., programs may use the message passing interface (MPI) communication protocol or OpenMP API to partition the domain into reasonable chunks divided between multiple processes, this can improve performance by having multiple processes solve smaller parts of a large problem in parallel.

Geodynamic inverse modelling involves using known phenomena to constrain a model, or set of models, to improve knowledge of unknowable or unmeasurable parameters involved within a dynamic system. For example, one could use measured GPS velocities and seafloor spreading rates to infer changes to mantle convection or combine GPS velocities with near-surface temperatures to monitor the growth of magma chambers.

Since geodynamic inversion involves running multiple constrained forward models it naturally suffers from the performance issues linked to the forward models themselves, briefly: a poorly performing forward model will mean a poorly performing inverse model. More than that, the inverse model must perform multiple (the more the better) forward models with updated parameters. These requirements, a good enough forward model, and an efficient method of performing multiple models in parallel while optimizing the desired parameters, leads to an obvious conclusion: to perform forward models in parallel while searching for an optimal model.

To this end we present a geodynamic inversion model, which we call ShellSetHPC, which uses a combination of existing, well known, and robust software to model the neotectonics of planetary lithosphere. This is further combined with an efficient, de-centralised random search algorithm able to generate testable models within a user defined N-dimensional space. This algorithm is also able to launch multiple models in parallel thanks to an MPI framework. The forward model makes use of Intel’s thread safe math kernel library (MKL) to solve the linear system in parallel. This hybrid approach lends itself to use on high performance computing (HPC) machines which would allow a more complete utilisation of these features.

In this work we will present scaling results on an HPC cluster and compare these with results obtained from more typical search algorithms. All tests are performed within a realistic geological setting, the results of which will be used to gather insight into the performance of the driving model generators and their settings.

How to cite: May, J. B., Bird, P., and Carafa, M. M. C.: ShellSetHPC – Parallel dynamic neotectonic modelling, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-11777, https://doi.org/10.5194/egusphere-egu24-11777, 2024.

Landlab is an open-source Python package that streamlines the creation, combination, and reuse of 2D numerical models and is a key element of the Community Surface Dynamics Modeling System (CSDMS) Workbench. The Landlab Toolkit provides building blocks for model development such as grid data structures, input/output functions, and a library of several dozen components that each model a separate physical process. Additionally, it provides a framework for assembling integrated models from component parts. We've found that Landlab significantly accelerates model development, encourages user-developers to adopt standard practices and contribute new components to the library. It serves as a platform that nurtures a community of model developers, assisting them in creating coupled models to investigate non-linear interactions between geologic processes.

Using the Landlab toolkit, we developed a new model, Sequence, which is a modular 2D (i.e., profile) sequence stratigraphic model that incorporates key geophysical processes influencing accommodation space in both terrestrial and marine environments. These factors include tectonics and faulting, eustatic sea level changes, flexural isostatic compensation of sediment and water, sediment compaction, and hypopycnal sediment plumes. Each process is encapsulated as an individual, standalone Landlab component, providing flexibility in the construction of new models. Sequence serves not only as a distinct model, but also as a scaffold for the development of new models.

Sequence simulates the evolution of stratigraphy on a continental margin over time scales ranging from thousands to millions of years. Sediment transport and deposition primarily occur during infrequent, high-energy events like storms and floods. For these extended time frames, Sequence employs a scale-integral approach. This method utilizes differential equations to summarize the cumulative effect of sediment transport and deposition across different depositional environments over longer periods (e.g., on the order of a hundred years). The model features a moving-boundary formulation to track shoreline changes and partitions the domain into distinct areas: coastal plain, continental shelf, and upper and lower slope/rise. Submarine sediment transport and deposition are modeled through nonlinear diffusion, with a diffusion coefficient that varies inversely with water depth. The model tracks evolving stratigraphic layers and sediment lithology that is a mixtuer of two grain sizes (sand and mud) each with separate transport functions.

How to cite: Hutton, E., Steckler, M., and Tucker, G.: Sequence: A coupled sequence-stratigraphic model built using Landlab, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-14445, https://doi.org/10.5194/egusphere-egu24-14445, 2024.

We present a high-performance Material Point Method (MPM) numerical solver for simulating the landslide run-out process with high resolution. The solver has been adapted to run on multi-GPU platforms. The current version is backend-agnostic, operating efficiently across various CPU and GPU hardware from different vendors, utilizing the same codebase. We evaluate multiple performance metrics and ensure minimal data synchronization between different devices at each iteration. We validate the solver's accuracy by comparing the simulation results of an aluminum-bar collapse with the corresponding experimental outcomes. Consistency is observed between numerical and experimental results for the free and failure surfaces. The results also indicate favorable performance scalability in high-resolution models, significantly enhancing computational efficiency. Further, we simulate the slumping mechanic problem with a simplified landslide geometry. The results show that the shear band develops within the high plastic strain area. The failure surface is in good agreement with the solution reported by Huang^[1] et al., demonstrating that MPM can accurately handle failure and large deformation problems as they occur in landslides. Our multi-GPU implementation using MPI makes it possible to perform large-scale simulations that enable to tackle research in the field of geotechnical engineering.

References:
[1]. Huang, Peng, Shun-li Li, Hu Guo, and Zhi-ming Hao. “Large Deformation Failure Analysis of the Soil Slope Based on the Material Point Method.” Computational Geosciences 19, no. 4 (August 2015): 951–63. https://doi.org/10.1007/s10596-015-9512-9.

How to cite: Huo, Z., Jaboyedoff, M., Podladchikov, Y., Räss, L., and Wyser, E.: Multi-GPU Material Point Method Solver for Landslide Simulation, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-15096, https://doi.org/10.5194/egusphere-egu24-15096, 2024.

Geophysical datasets and their interpretations form the basis of geodynamic simulations of the Earth’s mantle and lithosphere. Yet, creating geodynamic models from these datasets is often non-trivial, particularly in complex regions such as active orogens. Creating self-consistent three-dimensional models from these datasets poses a challenge due to technical issues such as different data set formats and different spatial resolutions, but also due to discrepancies in the data itself. At the same time, the different datasets obtained through initiatives such as AlpArray contain a wealth of data that can help to constrain subsurface models to an unprecedented extent. Yet interpreting these different data still involves subjective steps and ideally different datasets are combined in the process.

To facilitate the joint interpretation of these datasets and the generation of geodynamic model setups, we therefore developed an open-source package - the Geophysical Model Generator (GMG) - to assist with unifying these datasets in a common data format that can then be further used to visualize, compare and interpret data. Within this package, we provide a set of routines to import different datasets, convert them to a common data format and to process them further (e.g., to create vote maps from different tomographies). These unified datasets can then be exported as vtk-files for further 3D visualization (e.g., Paraview). Moreover, with the Geophysical Model Generator it is also possible to create model setups for numerical models (such as the 3D geodynamic code LaMEM). This package thus covers the entire workflow from data import to numerical model generation. Key features of the Geophysical Model Generator include 1) the creation of 3D volumes from seismic tomography models, 2) the import of 2D data (e.g., surface or Moho topography or screenshots from published papers) and 3) the incorporation of point data such as earthquake locations or GPS measurements. Both scalar and vector data can be handled. With these tools, one can then create a consistent overview of the entire data available for a given region.

The package is written in Julia and hosted as a public open-source repository on GitHub (https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl). To assist the joint interpretation of different geophysical datasets, we furthermore provide a graphical user interface that allows to view and compare them (https://github.com/JuliaGeodynamics/GeoDataPicker.jl). The GUI provides an interactive interface, allows loading different datasets and facilitates the manual interpretation of different structures (such as subducting slabs) along profiles and visualize them in 3D while taking different data into account. 

How to cite: Thielmann, M., Kaus, B. J. P., Spang, A., Schuler, C., and Aellig, P.: The Geophysical Model Generator: a tool to unify and interpret geophysical datasets, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-15998, https://doi.org/10.5194/egusphere-egu24-15998, 2024.

Julia is a relatively new scientific computing language with which you can write code that is nearly as fast as Fortran or C, but being a dynamic language, it remains very compact. Since the language is newer than say Python, many of the shortcomings of earlier languages have been fixed. Installing and updating packages with external dependencies is straightforward, for example, and exact reproducibility of computational results can be done by uploading a single additional file.

Here we will discuss our recent experiences with using Julia for teaching and research:

• The binarybuilder package simplifies the precompilation of non-Julia open-source packages (C, C++, Fortran), allowing cross-compilation for Apple, Linux, and Windows. Uploaded centrally, these precompiled binaries can be effortlessly installed via the Julia package manager. This expedites the setup of complex research code in teaching, reducing the time spent on package installation during the initial lecture to just a matter of minutes. Examples where this was recently implemented was GMT.jl (the Julia interface to the Generic Mapping Tools) which now comes with precompiled version of GMT. Likewise, we provide precompiled packages for LaMEM (a parallel 3D geodynamic code) and MAGEMin (a new thermodynamic code).
• It is quite easy to call such binary packages, either by calling the executable itself or by automatically creating wrapper (using Clang.jl) and calling the functions within dynamic libraries directly from julia.
• There is a large and ever-growing ecosystem of existing packages that are useful for computational geosciences. Whereas it is still not as mature as in Python or MATLAB, much of what is required in our daily life is already there, including plotting, machine learning, data I/O.
• It is easy to develop additional packages, and testing is build-in the language (and CI/CD is free on GitHub for open-source packages). That encourages users to implement tests and helps maintaining packages.
• Packages can talk to each other in a straightforward manner, which facilitates building composable software stacks.
• There is great support for GPU computing.
• It is possible to create GUI’s that run in a webbrowser (for example by using Dash.jl). Recent examples where we used that is in InteractiveGeodynamics.jl (which provides user interfaces to simulate typical geodynamic problems such as convection), GeoDataPicker.jl (a 2D/3D tool to analyze and visualize geodynamic data) and MAGEMin_app (a web-based tool to compute phase diagrams).
• It can directly be used in Jupyter notebooks (as the “Ju” in its name stands for Julia).
• It is easy for students to start writing simple codes without having to think about definition of types or having to load packages to create vectors or matrixes, which makes it well-suited for teaching quantitative classes. In many cases, these simple codes already run very fast, as the Julia compiler does an excellent job in optimisations. With a few additional tricks (ensuring type stability and reducing allocations) the simple examples can be turned into a high-performance code that runs at nearly the speed of the more classical languages, while having the advantage to still be very compact and readable.

We will highlight these aspects including with interactive demos.

How to cite: Kaus, B., de Montserrat, A., Riel, N., and Aellig, P.: Using Julia for teaching and research in the Geosciences, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-16725, https://doi.org/10.5194/egusphere-egu24-16725, 2024.

The study of coupled Earth systems, and in particular the coupled interactions between the lithosphere, atmosphere, and biosphere, have received greater attention in recent years (Gerya et al. 2020). Interactions between these systems occur primarily at the surface, and are driven on the large scale by topographic and bathymetric evolution controlled by deep mantle processes. However, due to the large difference in length scales between the mantle and the surface, it is difficult to capture topographic evolution to a high degree of accuracy in existing global mantle convection models including a free surface boundary condition.

Global mantle convection models incorporating a free surface often employ a marker-in-cell technique with a layer of “sticky air” (i.e. material with the density of the air and sufficiently low viscosity, which is still much higher than that of real air) to characterise the surface. However, accurate topographic evolution using this method requires a high density of markers near the surface. This need for additional computational resources motivates alternative methods of tracking the interface between the air and rock layers, as is done frequently in existing multiphase fluid flow codes. We demonstrate the implementation of two such methods of modelling the surface. The first is a Lagrangian marker chain (or mesh in 3D models) which, when combined with an appropriate remeshing procedure, directly tracks the rock-air interface. The second is a volume of fluid approach adapted from the open source code gVOF using the unsplit volume of fluid library gVOF (López & Hernández, 2022).

We demonstrate toy models and benchmarks (based on those in Crameri, 2012) comparing the Lagrangian marker method and the volume of fluid methods as implemented in the global scale mantle convection code StagYY (Tackley, 2008). Models of global scale topography and evolution produced using StagYY may then be used as a tool for further studies on the coupling of mantle dynamics with modelling of the landscape, and the evolution of the atmosphere and biosphere. Initial applications include modelling hypsometric curves from global scale numerical models, and the tracking of sea level changes over time.

How to cite: Gray, T., Gerya, T., and Tackley, P.: Free surface methods applied to global scale numerical geodynamic models, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-18423, https://doi.org/10.5194/egusphere-egu24-18423, 2024.

Developing a deep understanding of magmatic processes, such as the ascent of magmatic melts through the lithosphere, is a notoriously complex interdisciplinary task that involves contributions from various branches of geosciences. Here, we focus on the implementation of mode-1 plasticity, which is highly relevant for the modeling of dyke propagation through the brittle crust under an excess of fluid (melt) pressure, as part of a nonlinear elasto-visco-plastic rheology that is also appropriate for modeling the ductile-brittle transition in the upper mantle.

For this, we adopt a coupled poro-elasto-visco-plastic description of the strain localization process to capture the onset and advance of vertical dykes originating from magma reservoirs located at various depth. We mostly aim to investigate the hydro-mechanical interactions of such a system (e.g. permeability increase and elastic modulii degradation with increasing plastic strain, and influence of fluid pressure on both total stresses and on the yield surface). Thermal effects are essentially ignored, apart from the background geothermal gradient. The reservoir is represented as a cavity subjected to an increased fluid pressure.

We describe the brittle deformation using multi-surface plasticity models in the framework of the flow plasticity theory. A number of challenging problems are usually encountered in this case from an algorithmic viewpoint, including a lack of convexity and continuity of both the yield surface and flow potential, spurious elastic domains, singularity points and loss of convergence. We address these issues by using a relatively simple Perzyna-type visco-plastic model that consists of a linear Drucker-Prager shear failure envelop, and a circular tensile cap. Both surfaces are combined with each other in a way that enforces dimensional consistency, convexity, and continuity throughout the entire stress space.

Finally, we discuss the algorithmic details of the model which is incorporated in an implicit finite element code and demonstrate the  application results.

How to cite: Popov, A. A., Berlie, N., and Kaus, B. J. P.: Modeling fluid-induced dyke propagation in elasto-visco-plastic rocks, EGU General Assembly 2024, Vienna, Austria, 14–19 Apr 2024, EGU24-12283, https://doi.org/10.5194/egusphere-egu24-12283, 2024.