Big data, sparse information: Bayesian inference for large-scale models, with application to inverse modeling of Antarctic ice sheet dynamics; Omar Ghattas, Center for Computational Geosciences, The University of Texas at Austin, USA

Predictive models of complex geoscience systems often contain numerous uncertain parameters. Rapidly expanding volumes of observational data present opportunities to reduce these uncertainties via solution of inverse problems. Bayesian inference provides a systematic framework for inferring model parameters with associated uncertainties from (possibly noisy) data and any prior information. However, solution of the Bayesian inverse problem via conventional Markov chain Monte Carlo methods remains prohibitive for expensive models and high-dimensional parameters, as result from discretization of infinite dimensional problems with uncertain fields. Despite the large size of observational datasets, typically they can provide only sparse information on model parameters. Based on this property we design MCMC methods that adapt to the structure of the posterior probability and exploit an effectively-reduced parameter dimension, thereby making Bayesian inference tractable for some large-scale, high-dimensional inverse problems. We apply the methodology to an inverse problem for Antarctic ice sheet flow: given a non-Newtonian flow model of the ice sheet and free surface velocity observations, we seek to infer the sliding coefficient field at the base of the ice.

This work is joint with Tan Bui-Thanh, Tobin Isaac, James Martin, Noemi Petra, Georg Stadler, Hari Sundar, and Hongyu Zhu.

Opportunities and challenges in large-scale adaptive-mesh simulations; Carsten Burstedde (Institute for Numerical Simulation, University of Bonn), Tobin Isaac, Bram Metsch

Simulations in the geosciences tend to become large in terms of the number of unknowns that are needed to reach a desired accuracy. Distributing the unknowns uniformly in space leads to simple and efficient bookkeeping, but harbors the risk of placing the majority of them where they are not strictly required. This risk is pronounced when the simulated system exhibits strongly localized phenomena, for example around fractures or fronts, where the density of unknowns should be concentrated.

Distributing the unknowns adaptively is a natural approach to balance the accuracy requirements with the associated computational cost. However, it is also highly demanding in practice. Challenges arise with respect to designing efficient parallel mesh refinement algorithms for distributed-memory computers, realizing the computer's peak performance to a high percentage, and optimally scaling multilevel implicit solvers. In this talk, we outline algorithmic concepts that offer solutions to these issues, present recent large-scale results, and discuss how we can realize adaptivity with only minimal overhead compared to state-of-the-art uniform-mesh simulations.

Challenges in the Simulation of Earthquakes; Christian Pelties (LMU)

In my presentation, I will discuss several crucial aspects of the numerical simulation of earthquakes. To this end, established and state of the art approaches used in solving the wave equation and non-linear frictional sliding are briefly discussed. I will also point to current issues and open questions.

Simulating earthquakes is used for improving our understanding of earthquake physics and for preparedness purposes of potential future events. Thereby, a non-linear friction law is evaluated on meter scale that describes the frictional sliding along an interface generating seismic waves that propagate over hundreds of kilometers. For realistic simulations, many elements as geometric details, geological properties and involved physics must be appropriately considered.

Such a computer experiment poses many challenges on all parts of the software, e.g. flexibility and accuracy of the underlying numerical method, scalability on modern HPC infrastructure, as well as parallel output techniques. Namely, numerical artifacts as spurious oscillations generated during the frictional sliding need to be relaxed by the solver without disturbing the physical modeling process. The solver accuracy has to be high and dispersion errors low to correctly simulate waves propagating over large distances. Geometric complexity in form of faults, topography or material interfaces that merge often under shallow angles and contain discontinuities have to be respected. More difficulties occur in the appropriate sampling of small scale stress fluctuations at the fault and the comprehension of scatterer in the medium. Further questions are: How to prepare for petaflop performance machines to run sufficiently large meshes to achieve the desired high frequencies of the seismic wave field? What is the best way to implement modern and heterogeneous algorithms efficiently on todays CPU’s? How to deal with the diversity of existing accelerators as GPU’s or Intel MIC architecture? Are resilience and parallel I/O essential?

HPC Solutions for Geophysical Exploration Problems: Efficiency, Reliability and Maintenance for Real-Life Applications; Josep de la Puente, Mauricio Hanzich, Natalia Gutierrez, Juan E. Rodriguez, Jean Kormann (BSC)

The field of exploration geophysics has recently embraced supercomputing as an essential part of its data processing workflow. In fields such as hydrocarbon exploration, CO2 storage or geothermal reservoir monitoring, simulation-based 3D imaging tools are now commonplace both for the industry and academia. However, due to the large number of possible geophysical imaging tools, physical approximations, numerical schemes and HPC architectures, the development of a comprehensive suite that is reliable, maintainable and efficient is a challenging task. Here, I will describe our vision which involves using a code framework structure that allows us to tackle a large set of geophysical problems with a reduced development cost and high efficiency payoff. Although compromises must always be made, I will present results in large-scale 3D full-waveform inversion, elastic reverse-time migration and controlled-source electromagnetic modelling which suggest that a multiple-purpose computing code offers many advantages in terms of reliability and portability. Furthermore, this framework-based code structure allows for novel concepts to be quickly incorporated and tested in complex real-life scale scenarios, hence blurring the line between academic research and production-ready applications.

A Semismooth Newton-CG Method for Constrained Parameter Identification in Seismic Tomography; Christian Boehm, Michael Ulbrich (Technische Universität München)

Full-waveform seismic tomography can be stated as a PDE-constrained optimization problem governed by the elastic wave equation. One way to incorporate prior knowledge and to address the ill-posedness of the inverse problem is to impose additional constraints on the material parameters. We present a semismooth Newton-PCG method with a trust-region globalization that uses a Moreau-Yosida regularization to handle the constraints. Forward and adjoint simulations are carried out to compute the reduced gradient and reduced Hessian-vector products which are required by the Newton-CG method. Furthermore, we combine ideas from optimal experimental design and stochastic programming to employ randomized data reduction techniques and efficiently gather information from a large number of seismic events. In particular, we utilize randomly chosen linear combinations of the seismic sources to approximate the Hessian operator and to reduce the number of required simulations during the CG iterations. Numerical results are shown for an application in geophysical exploration on reservoir-scale and the joint inversion for both Lamé coefficients with constraints on the Poisson's ratio.

Physics-based earthquake scenarios for hazard assessment in large urban areas; Ilario Mazzieri (Politecnico di Milano), Paola F. Antonietti (Politecnico di Milano), Alfio Quarteroni (EPF-Lausanne), Marco Stupazzini (Munich-RE), Roberto Paolucci (Politecnico di Milano), Chiara Smerzini (D'Apollonia-Genova)

In recent years, boosted by the ever-increasing power of large parallel computing facilities, numerical-based deterministic methods, relying on the extensive use of 3D numerical simulations of seismic wave propagation, are becoming more and more appealing for providing realistic earthquake ground shaking maps to be used as the key ingredient to develop earthquake hazard assessment.

In this work we provide an overview of the recent advances on this subject, and present a set of realistic physics-based ground shaking scenarios in key locations worldwide obtained with SPEED, a parallel open source code for the simulation of complex wave propagation phenomena, jointly developed by the Laboratory for Modeling and Scientific Computing MOX (Department of Mathematics, Politecnico di Milano) and the Engineering Seismology Group Department of Civil and Environmental Engineering, Politecnico di Milano). From the methodological viewpoint, SPEED is based on a spectral element formulation enhanced by the Discontinous Galerkin approach for treating non-conforming approximations.

After illustrating the SPEED code, we present different case studies and we draw some interesting considerations on the pros and on the limitations of deterministic approaches.

This work has been carried out in collaboration with Munich RE (Germany).

Parallel Adaptive Simulation of Earthquake-Tsunami Events; Michael Bader (TU München), Christian Pelties (LMU München), Jörn Behrends (University of Hamburg)

We will discuss parallelization and performance issues with large-scale simulation of earthquake-tsunami events on supercomputing platforms. Our goal in the respective collaborative research project ASCETE (see www.ascete.de) is to establish a simulation pipeline that combines highly detailed simulations of the dynamic rupture process of tsunamigenic earthquakes and computation of the resulting time-dependent displacements with dynamically adaptive simulation of tsunami generation and propagation.

Earthquake simulation is based on the software package SeisSol: SeisSol solves the 3D elastic wave equations in velocity-stress formulation, using an Arbitrary high order Derivatives Discontinuous Galerkin (ADER-DG) discretization. The resulting high computational intensity and the comparably simple communication patterns due to the DG scheme make SeisSol an excellent candidate for large-scale parallelism. SeisSol works on unstructured adaptive tetrahedral meshes, which allows to handle complicated geometries and fault planes that meet the surface at shallow angles (as typically occurring for subduction zone settings). Together with non-linear dynamic rupture models applied on these fault systems, highly detailed simulations of the entire process of tsunamigenic earthquakes become possible. We will report on recent work to realize petascale simulations with SeisSol under production conditions, and will also present latest results on optimizing SeisSol for manycore architectures such as the Intel Xeon Phi.

To simulate tsunamis resulting from computed (time-dependent) displacements of the ocean floor, we use the software framework sam(oa)². sam(oa)² features parallel adaptive mesh refinement and dynamically adaptive simulations based on a Sierpinski space-filling-curve approach. To simulation wave propagation and inundation, we solve the shallow water equations on fully adaptive triangular grids using a Finite Volume approach with different approximate Riemann solvers. A discontinuous Galerkin implementation is work in progress. Mesh refinement at propagating wave fronts and at coast lines is performed after each time step. For parallelization and load balancing, we exploit the element order defined by the Sierpinski curve, which allows load balancing after each time step at a moderate overhead. Dynamically adaptive tsunami simulations with more than a billion elements are thus possible on department-size clusters.

Modelling the solid Earth plate tectonics-mantle convection system and its feedbacks with the geodynamo, atmosphere and ocean; Paul J. Tackley (ETH Zurich)

Convection of the solid, rocky mantle and related plate tectonics are the mechanisms driving very long-term (millions to billions of years) evolution of both Earth’s interior and the atmosphere-ocean system via volcanic outgassing and recycling of CO2 and H2O. Mantle convection also controls the geodynamo by determining the heat flux extracted from the core. Mantle convection and plate tectonics are, however, challenging to model due to the large range of lengthscales (faults of width ~cms to continents 1000s km wide) and timescales (seconds to billions of years), and the multi-rheological nature of rocks (they can undergo elastic, brittle, ductile and/or plastic deformation depending on the temperature, pressure and stress). Here the state of the field is reviewed, with particular examples from the research of the presenter’s group.

It is now possible to perform global 3D spherical models of solid Earth evolution that span 4.5 billion years, include self-consistently generated plate tectonics and magmatism that generates crust and chemically differentiates the mantle. Such calculations make predictions that can be compared to observations from seismology, geochemistry, paleomagnetism and geology. Some recent models also include coupling to the atmosphere (outgassing + recycling). Despite such progress at the global scale, the processes that allow weak plate boundaries to form and persist are not well understood and require detailed local and regional models to unravel. We are also applying such models to other planets, both in our own solar system (Mars, Venus, Mercury, Io) and Earth-like planets around other stars – so-called super-Earths.

Much of this current modelling is performed using conventional finite-volume and marker-in-cell methods. Continued future progress requires substantial technical development that can take advantage of modern supercomputer architectures and algorithms. Important future directions will be reviewed, including the goals of the GeoPC project.

Faster and more accurate interface tracking using a level set two-way wave equation approach; Henri SAMUEL (IRAP)

Interface tracking problems are ubiquitous in various branches of computational Earth sciences. For this purpose, the level set tracking method is a popular approach that can handle interfaces with complex topologies. However this approach, in its pure Eulerian version, tends to suffer from significant dissipation errors. As an alternative I will present a new approach to perform Eulerian level set interface tracking. It consists in advecting the level set function using a second-order two-way wave equation instead of the standard one-way wave advection equation. The resulting numerical schemes are simple to implement, more stable, and significantly less prone to dissipation errors than popular (e.g., WENO) discretizations of the one-way advection equation of higher order, in particular for long advection times. While both the two-way wave advection and the associated numerical schemes were derived previously, these approaches have never been combined with level set advection. The level set function designed to be smooth ensures the stability of the solution when using the two-way advection equation discretized with centered finite difference schemes. Numerical tests in multi-dimensional space show that the two-way wave equation approach yields more accurate results than the standard one-way wave equation for the level set advection, at a considerably smaller (e.g., one order of magnitude) computational cost. Although the development of this new approach was originally motivated for problems in solid Earth dynamics, its versatility makes it readily suitable for a broad range of applications ranging from general CFD problems to computer vision or tumor growth modelling.

Immersive experimentation in a wave propagation laboratory; Johan O. A. Robertsson, Dirk-Jan van Manen, Stewart Greenhalgh, Thomas Blum, and Marlies Vasmel (ETH Zurich)

Physical laboratory-based seismic modelling allows wave propagation to be studied in realistic media which may be beyond current computer modeling capability. Unfortunately, most physical modeling laboratories that exist to date are plagued with shortcomings that limit their use. In addition to coupling issues, the finite size and directivity patterns of transducers, avoiding both the generation and influence of undesired boundary reflections require laboratory targets that are several orders of magnitude greater than the wavelengths studied. Hence, investigating frequency-dependent velocity and attenuation behavior in porous hydrocarbon reservoir rocks is largely beyond what current physical seismic modeling laboratories can accommodate and is limited to rock deformation experiments to characterize shear and bulk moduli. By surrounding a physical modelling experiment with closed surfaces populated by discrete transducers that can transmit and record signals, it is possible to completely link a physical wave propagation experiment with a numerical experiment, such that physical waves generated in the laboratory can interact with the numerical model and vice-versa. At ETH we are in the process of building such a Wave Propagation Laboratory (WaveLab). The physical experiment is coupled to (or immersed in) a numerical experiment in real-time through discretized convolutions with pre-computed Green’s functions (characterizing the medium outside) carried out on a high performance compute system consisting of GPUs. The real-time interaction is taking place using the theory of exact boundary conditions (EBCs), by continuously recording and emitting the outgoing and incoming wavefields on surfaces completely surrounding the physical experiment and lining the walls of the laboratory tank, respectively. Both the speed of computations as well as minimizing possible latencies of the system are critical aspects for the construction of the WaveLab. FPGAs perform time-critical pre- and post-processing tasks and communicate with the acquisition system and GPUs through fast peer-to-peer connections. In this way, the walls of the tank become transparent to outward propagating waves and the physical experiment completely behaves as if embedded in a larger medium. The WaveLab, in principle, eliminates most of the shortcomings of physical modelling. In particular the new closed loop laboratory acquisition system avoids artificial and undesired reflections from the walls of the tank. More important, however, it will enable the reverberation of energy from “virtual structures” whose presence is only reflected through pre-computed (or measured) Green’s functions. We envisage a large number of applications of the WaveLab ranging from fundamental studies of the physics of wave propagation to solving particular problems related to, for instance, upscaling or wave propagation in partially fluid-saturated media. The WaveLab is particularly suited to the study of wave propagation problems involving broadband signals and/or scenarios where the physics of wave propagation is neither well understood nor linear. Finally, the WaveLab will enable revisiting acoustic and elastic time-reversal experiments performed in the 1990s (due to it being a completely closed and instrumented cavity), offering the prospect of near-perfect time reversal. It will allow verification and experimental demonstration of recently proposed data-driven focusing techniques, which have the potential to revolutionize acquisition and imaging of exploration seismic data as well as ablation and imaging of medical data.

Constraining the Rheology of the Lithosphere Through Joint Geodynamic and Gravity Inversion; Boris Kaus, Tobias Baumann, Anton Popov (Johannes-Gutenberg University Mainz, Germany)

Understanding the physics of lithospheric deformation requires good constraints on lithospheric rheology and in particular on the effective viscosity. Typically, rheology is determined from laboratory experiments on small rock samples, which are extrapolated to geological conditions - an extrapolation over 10 orders of magnitude in deformation rates. In addition, lithospheric rheology can be strongly nonlinear and varies from powerlaw viscous at higher temperatures to plastic at lower temperatures.

Ideally, we thus need a new independent method that allows constraining the effective rheology of the lithosphere directly from geophysical data, which is the aim of this work.

Our method uses the fact that the geodynamically controlling parameters of lithospheric deformation are its effective viscosity and density structure. By appropriately parametrising the rheological structure of the lithosphere we perform instantaneous forward simulations of present-day lithospheric deformation scenarios with a finite element method to compute the gravity field as well as surface velocities. The forward modelling results can be compared with observations such as Bouguer anomalies and GPS-derived surface velocities. More precisely, we automatise the forward modelling procedure with a Markov-Chain Monte Carlo method, and in fact solve a joint geodynamic and gravity inverse problem. For this purpose we have developed a parallel thermomechanical forward modelling code, which is based on a marker-and-cell approach and a staggered finite difference discretization. We have combined this with a parallel inversion approach with non-blocking communication to deal with the fact that the runtimes between individual forward models varies significantly depending on the effective viscosity structure.

We will discuss the computational approach of the forward and inverse modelling as well as show synthetic 2D and 3D examples that show that the method works even for strongly nonlinear rheologies. In addition we show an application of the method to natural data and constrain the rheology of the lithosphere at the India-Asia collision zone.

Acknowledgements Funding was provided by the ERC under the European Community's Seventh Framework Program (FP7/2007-2013) / ERC Grant agreement #258830

Towards magnetic sounding of the Earth's core by an adjoint method; Kuan Li, Andrew Jackson (ETHZ), Philip W. Livermore (University of Leeds)

Earth's magnetic field is generated and sustained by the so called geodynamo system in the core. Measurements of the geomagnetic field taken at the surface, downwards continued through the electrically insulating mantle to the core-mantle boundary (CMB), provide primary constraints on the time evolution of the velocity, magnetic field and temperature anomaly in the fluid outer core. The aim of any study in data assimilation applied to the Earth's core is to produce a time-dependent model consistent with these observations [Fournier et al, 2010]. Snapshots of these "tuned" models provide a window through which the inner workings of the Earth's core, usually hidden from view, can be probed.

We apply a variational data assimilation framework to an inertia-free magnetohydrodynamic system (MHD) [Glatzmaier and Roberts, 1995]. Such a model is close to magnetostrophic balance [Taylor, 1963], to which we have added viscosity to the dominant forces of Coriolis, pressure, Lorentz and buoyancy, believed to be a good approximation of the Earth's dynamo in the convective time scale. We chose to study the MHD system driven by a static temperature anomaly to mimic the actual inner working of Earth's dynamo system, avoiding at this stage the further complication of solving for the time dependent temperature field. At the heart of the models is a time-dependent magnetic field to which the core-flow is enslaved.

In our previous work we laid the foundation of the adjoint methodology, applied to a subset of the full equations [Li et al, 2011]. As an intermediate step towards our ultimate vision of applying the techniques to a fully dynamic mode of the Earth's core tuned to geomagnetic observations, we present the intermediate step of applying the adjoint technique to the inertia-free Navier-Stokes equation in continuous form. We use synthetic observations derived from evolving a geophysically-reasonable magnetic field profile as the initial condition of our MHD system. Based on our study, we also propose several different strategies for accurately determining the entire trajectory of Earth's geodynamo system.

We solve this MHD system and its adjoint using a spectral method, due to its rapid convergence, where the magnetic field and the incompressible flow are represented by poloidal-toroidal fields and are further expanded in terms of spherical harmonics and specially constructed radial basis functions satisfying the relevant boundary conditions associated with the magnetic field and the flow. At the current stage of the numerical study of the geodynamo system, tens of thousands to millions of spectral coefficients for the magnetic and velocity field are to be solved at each time step and therefore can only be computed on a supercomputer with processors numbering several hundred to a thousand processors. Central to our numerical model are the facts that (i) the magnetic induction equation and its adjoint are parallelized in radius and (ii) the simplified Navier-Stokes equation is linear and diagnostic and is parallelised by symmetry class. The minimization step is computed on a single processor using the L-BFGS package, where at each assimilation cycle, the misfit and the adjoint field (representing the downhill direction) are collected from each processor.