GD
Geodynamics

geodynamics

Get conference ready!

Get conference ready!

It’s almost time for the AGU fall meeting 2018! Are you ready? Have you prepared your schedule and set up all your important business meetings? Here are some final tips to nail your presentation and/or poster!

Nailing your presentation
The art of the 15-minute talk: how to design the best 15-minute talk
Presentation skills – 1. Voice: how to get the most out of your presentation voice
Presentation skills – 2. Speech: how to stop staying ‘uh’

Making the best poster
Poster presentation tips: how to design the best poster layout
The rainbow colour map (repeatedly) considered harmful: how to make the best scientific figures

An industrial placement as a geodynamicist

An industrial placement as a geodynamicist

After years of trying to get a PhD, publishing papers, networking with professors, and trying to land that one, elusive, permanent job in science, it can be quite easy to forget that you actually do have career options outside of academia. To get a little taste of this, Nico Schliffke, PhD student in geodynamics at Durham University, tries out the industry life for a few weeks!

When coming close to the final stages of a PhD life, many students reconsider whether they want to stay in academia or prefer to step over to industry or other non-academic jobs. This is surely not a simple decision to take, as it could strongly shape your future. In this blog post, I would like to report my industrial placement experience during my PhD and share a few thoughts on the topic.

The taste of industry life was an opportunity I had within the frame of my PhD project. Split into two terms, I spent four weeks at a medium-sized company developing optical imaging techniques (both software and equipment) to measure flow fields and deformation. The branch I worked in was “digital image correlation” (DIC) which measures strain on given surfaces purely by comparing successive images on an object (see figure below). This technique is used in both industry (crash tests, quality assessments, etc.) as well as in academia (analogue experiments, wind tunnels, engineering..), and has the substantial advantage of measuring physical properties precisely, without using any materials or affecting dynamical processes. DIC is not directly related to or used in my PhD (I do numerical modelling of subduction zones and continental collision), but surprisingly enough I was able to contribute more than expected – but more on that later.

Basic principle of ‘digital image correlation’. A pattern on a digital image is traced through time on successive images to calculate displacements and strain rates.
Credit: LaVision

The specific project I worked on was inspired by the analogue tectonics lab at GFZ Potsdam, that uses DIC measuring systems to quantify and measure the deformation of their sandbox experiments. Typical earthquake experiments like the figure below span periods of a few minutes to several days during which individual earthquakes occur in a couple of milliseconds. The experiment is continuously recorded by cameras to both monitor deformation visually and quantify deformation by using the optical imaging technique developed by my host company. To resolve the deformation related to individual earthquakes, high imaging rates are required which in turn produce a vast amount of data (up to 2TB per experiment). However, only a small fraction (max. 5%) of the entire dataset is of interest, as there is hardly any deformation during interseismic periods. The project I was involved in tried to tackle the issue of unnecessarily cluttered hard discs: the recording frequency should be linked to a measurable characteristic within the experiment, e.g. displacement velocities in these specific experiments, and controlled by the DIC software.

Setup of the analogue experiment to model earthquakes in subduction zones (courtesy of Michael Rudolf). Cameras above the experiment measure deformation and strain rates by tracking patterns on the surface created by the contrast of black rubber and white sugar.

My general task during the internship was to develop this idea and the required software. We finally developed a ‘live-extensometer’ to calculate displacements between two given points of an image during recording and link its values to the camera’s recording frequency. Therefore, restricting high imaging rates to large (and fast) displacements of earthquakes should result in reducing the total amount of data acquired for a typical earthquake experiment by 95%. However, we needed an actual experiment to verify this. So, I met up with the team at GFZ to test the developed feature.

The main experiment the GFZ team had in mind is sketched in the figure above: a conveyor belt modelling a subducting slab continuously creates strain in the ‘orogenic wedge’ which is released by earthquakes leading to surface deformation. Cameras above the experiment monitor the surface while software computes strain rates and displacement (see figure below). The developed feature of changing frequencies during the experiment depending on slip rates was included and worked surprisingly well. Yet freshly programmed software is seldom perfect: minor issues and bugs crept up during the experiments. My final contribution during the internship was to report these problems back to the company to be fixed.

Displacement measured by ‘digital image correlation’ during an earthquake lasting ~5 ms (courtesy of Mathias Rosenau).

My geodynamical background allowed me to contribute to various fields within the company and resulted in various individual tasks throughout the internship: coding experience helped with discussing ideal software implementations and testing the latest implemented software on small (physical) experiments. My knowledge of various deformation mechanisms and geosciences in general, with its numerous subdisciplines and methods, provided a solid base for searching further applications for the developed software within academia, but also in industry. Last but not least, pursuing my own large project (my PhD) strongly facilitated discussing possible future development steps.

The atmosphere at the company in general was very pleasant and similar to what I experienced at the university: relaxed handling, pared with discussion how to improve products or use of new techniques that might be applicable to a problem. To stay competitive, the company needs to further develop their products which requires a large amount of research, developments and innovative ideas. Meetings to discuss further improvements of certain products were thus scheduled on a (nearly) daily basis. On the one hand this adds pressure to get work done as quickly as possible, but working on a project as a team with many numerous areas of expertise is also highly exciting.

This internship help reveal the variability of possible jobs that geodynamicists can have in industry besides the ‘classical’ companies linked to exploration, tunnel engineering or geological surveys. The skill set acquired in a geodynamical PhD (coding, modelling, combining numerics, physics, and geosciences) makes a very flexible and adaptive employee which is attractive to companies who are so specialised, that there is (nearly) no classical education at university level. Jobs at small to medium-sized companies are often harder to find, but it’s just as difficult for the companies to find suitable candidates for their open positions. Hence, it may be worth searching in-depth for a suitable job, if you are considering stepping out of academia and maybe even out of geoscience as well.

If PhD students are hesitant whether to stay in academia or change into industry, I would advise to do such a short internship with a company to get a taste of ‘the other side’. During a PhD, we get to know academic life thoroughly but industry mostly remains alien. Besides giving a good impression of daily life at a company and how you can contribute, an industry internship might also widen your perspective of which areas might be relevant to you, your methodology and your PhD topic. In total, this internship was definitely a valuable experience for me and will help when deciding: academia or industry?


Here are a few links for more information:
Host company
Digital Image Correlation
TecLab at GFZ Potsdam
Previous EGU blog post interviews of former geoscientists


Inversion 101 to 201 – Part 3: Accounting for uncertainty – Bayes and friends

Inversion 101 to 201 – Part 3: Accounting for uncertainty – Bayes and friends

The Geodynamics 101 series serves to showcase the diversity of research topics and methods in the geodynamics community in an understandable manner. We welcome all researchers – PhD students to professors – to introduce their area of expertise in a lighthearted, entertaining manner and touch upon some of the outstanding questions and problems related to their fields. This time, Lars Gebraad, PhD student in seismology at ETH Zürich, Switzerland, talks about inversion and how it can be used in geodynamics! Due to the ambitious scope of this topic, we have a 3-part series on this! You can find the first post here. This week, we have our final post in the series, where Lars discusses probabilistic inversion.

One integral part of doing estimations on parameters is an uncertainty analysis. The aim of a general inverse problem is to find the value of a parameter, but it is often very helpful to indicate the measure of certainty. For example in the last figure of my previous post, the measurement values at the surface are more strongly correlated to the upper most blocks. Therefore, the result of an inversion in this set up will most likely be more accurate for these parameters, compared to the bottom blocks.

In linear deterministic inversion, the eigenvalues of the matrix system provide an indication of the resolvability of parameters (as discussed in the aforementioned work by Andrew Curtis). There are classes of methods to compute exact parameter uncertainty in the solution.

From what I know, for non-linear models, uncertainty analysis is limited to the computation of second derivatives of the misfit functional in parameter space. The second derivatives of X (the misfit function) are directly related to the standard deviations of the parameters. Thus, by computing all the second derivatives of X, a non-linear inverse problem can still be interrogated for its uncertainty. However, the problem with this is its linearisation; linearising the model and computing derivatives may not be truly how the model reacts in model space. Also, for strongly non-linear models many trade-offs (correlations) exist which influence the final solution, and these correlations may very strongly depending on the model to be inverted.

Bayes’ rule

Enter reverend Thomas Bayes

This part-time mathematician (he only ever published one mathematical work) from the 18th century formulated the Bayes’ Theorem for the first time, which combines knowledge on parameters. The mathematics behind it can be easily retrieved from our most beloved/hated Wikipedia, so I can avoid getting to caught up in it. What is important is that it allows us to combine two misfit functions or probabilities. Misfits and probabilities are directly interchangeable; a high probability of a model fitting our observations corresponds to a low misfit (and there are actual formulas linking the two). Combining two misfits allows us to accurately combine our pre-existing (or commonly: prior) knowledge on the Earth with the results of an experiment. The benefits of this are two-fold: we can use arbitrarily complex prior knowledge and by using prior knowledge that is bounded (in parameter space) we can still invert underdetermined problems without extra regularisation. In fact, the prior knowledge acts as regularisation.

One probabilistic regularised bread, please

Let’s give our friend Bayes a shot at our non-linear 1D bread. We have to come up with our prior knowledge of the bread, and because we did not need that before I’m just going to conjure something up! We suddenly find the remains of a packaging of 500 grams of flour

This is turning in quite the detective story!

However, the kitchen counter that has been worked on is also royally covered in flour. Therefore, we estimate that probably this pack was used; about 400 grams of it, with an uncertainty (standard deviation) of 25 grams. Mathematically we can formulate our prior knowledge as a Gaussian distribution with the aforementioned standard deviation and combine this with our misfit of the inverse problem (often called the likelihood). The result is given here:

Prior and original misfit

Combined misfit

One success and one failure!

First, we successfully combined the two pieces of information to make an inverse problem that is no longer non-unique (which was a happy coincidence of the prior: it is not guaranteed). However, we failed to make the problem more tractable in terms of computational requirements. To get the result of our combined misfit, we still have to do a systematic grid search, or at least arrive at a (local) minimum using gradient descent methods.

We can do the same in 2D. We combine our likelihood (original inverse problem) with rather exotic prior information, an annulus in model space, to illustrate the power of Bayes’ theorem. The used misfit functions and results are shown here:

Original misfit for a bread of 500 grams

Prior knowledge misfit

Combined misfit

This might also illustrate the need for non-linear uncertainty analysis. Trade-offs at the maxima in model space (last figure, at the intersection of the circle and lines) distinctly show two correlation directions, which might not be fully accounted for by using only second derivative approximations.

Despite this ‘non-progress’ of still requiring a grid search even after applying probability theory, we can go one step further by combining the application of Bayesian inference with the expertise of other fields in appraising inference problems…

Jump around!

Up to now, using a probabilistic (Bayesian) approach has only (apparently) made our life more difficult! Instead of one function, we now have to perform a grid search over the prior and the original problem. That doesn’t seem like a good deal. However, a much used technique in statistics deals with exactly the kind of problems we are facing here: given a very irregular and high dimensional function

How do we extract interesting information (preferably without completely blowing our budget on supercomputers)?

Let’s first say that with interesting information I mean minima (not necessarily restricted to global minima), correlations, and possibly other statistical properties (for our uncertainty analysis). One answer to this question was first applied in Los Alamos around 1950. The researches at the famous institute developed a method to simulate equations of state, which has become known as the Metropolis-Hastings algorithm. The algorithm is able to draw samples from a complicated probability distribution. It became part of a class of methods called Markov Chain Monte Carlo (MCMC) methods, which are often referred to as samplers (technically they would be a subset of all available samplers).
The reason that the Metropolis-Hastings algorithm (and MCMC algorithms in general) is useful, is that a complicated distribution (e.g. the annulus such as in our last figure) does not easily allow us to generate points proportional to its misfit. These methods overcome this difficulty by starting at a certain point in model space and traversing a random path through it – jumping around – but visiting regions only proportional to the misfit. So far, we have only considered directly finding optima to misfit functions, but by generating samples from a probability distribution proportional to the misfit functions, we can readily compute these minima by calculating statistical modes. Uncertainty analysis subsequently comes virtually for free, as we can calculate any statistical property from the sample set.

I won’t try to illustrate any particular MCMC sampler in detail. Nowadays many great tools for visualising MCMC samplers exist. This blog by Alex Rogozhnikov does a beautiful job of both introducing MCMC methods (in general, not just for inversion) and illustrating the Metropolis Hastings Random Walk algorithm as well as the Hamiltonian Monte Carlo algorithm. Hamiltonian Monte Carlo also incorporates gradients of the misfit function, thereby even accelerating the MCMC sampling. Another great tool is this applet by Chi Feng. Different target distributions (misfit functions) can be sampled here by different algorithms.

The field of geophysics has been using these methods for quite some time (Malcom Sambridge writes in 2002 in a very interesting read that the first uses were 30 years ago), but they are becoming increasingly popular. However, strongly non-linear inversions and big numerical simulations are still very expensive to treat probabilistically, and success in inverting such a problem is strongly dependent on the appropriate choice of MCMC sampler.


Summary of probabilistic inversions

In the third part of this blog we saw how to combine any non-linear statistical model, and how to sample these complex functions using MCMC samplers. The resulting sample sets can be used to do an inversion and compute statistical moments of the inverse problem.


Some final thoughts

If you reached this point while reading most of the text, you have very impressively worked your way through a huge amount of inverse theory! Inverse theory is a very diverse and large field, with many ways of approaching a problem. What’s discussed here is, to my knowledge, only a subset of what’s being used ‘in the wild’. These ramblings of a aspiring seismologist might sound uninteresting to the geodynamiscists at the other side of the geophysics field. Inverse methods seem to be not nearly discussed as much in geodynamics as they are in seismology. Maybe it’s the terminology that differs, and that all these concepts are well known and studied under different names and you recognise some of the methods. Otherwise, I hope I have given an insight in the wonderful and sometimes ludicrous mathematical world of (some) seismologists.

Interested in playing around with inversion yourself? You can find a toy code about baking bread here.


Inversion 101 to 201 – Part 2: The inverse problem and deterministic inversion

Inversion 101 to 201 – Part 2: The inverse problem and deterministic inversion

The Geodynamics 101 series serves to showcase the diversity of research topics and methods in the geodynamics community in an understandable manner. We welcome all researchers – PhD students to professors – to introduce their area of expertise in a lighthearted, entertaining manner and touch upon some of the outstanding questions and problems related to their fields. This time, Lars Gebraad, PhD student in seismology at ETH Zürich, Switzerland, talks about inversion and how it can be used in geodynamics! Due to the ambitious scope of this topic, we have a 3-part series on this! You can find the first post here. This week, we have our second post in the series, where Lars discusses the inverse problem and deterministic inversion.

The idea of inversion is to literally invert the forward problem. For the aforementioned problems our knowns become unknowns, and vice versa. We wish to infer physical system parameters from measurements. Note that in our forward problem there are multiple classes of knowns; we have the forcing parameters, the material properties and the boundary conditions. All of these parameters could be inverted for in an inversion, but it is typically one only inverts for one class. Most of the time, medium parameters are the target of inversions. As we will go through the examples, I will gradually introduce some inversion methods. Every new method is marked in blue.

Estimating ingredients from a bread

Now, let’s consider our first example: the recipe for a bread. Let’s say we have a 0.5 kilogram bread. For our first case (case 1a), we will assume that the amount of water is ideal. In this case, we have one free variable to estimate (the used amount of flour), from one observable (the resulting amount of bread). We have an analytical relationship that is both invertible and one-to-one. Therefore, we can use the direct inversion of the relationship to compute the amount of flour. The process would go like this:

Direct inversion

1. Obtain the forward expression G(m) = d;
2. Solve this formula for G^-1 (d) = m;
3. Find m by plugging d into the inverse function.

Applying this direct inversion shows us that 312.5 grams of flour must have been used.

The two properties of the analytical relationship (invertible and one-to one) are very important for our choice of this inversion method. If our relationship was sufficiently complex, we couldn’t solve the formula analytically (though the limitation might then lie by the implementer 😉 ). If a function is not one-to-one, then two different inputs could have the same output, so we cannot analytically invert for a unique solution.

In large linear forward models (as often obtained in finite difference and finite element implementations), we have a matrix formula similar to

A m = d.

In these systems, m is a vector of model parameters, and d a vector of observables. If the matrix is invertible, we can directly obtain the model parameters from our observables. The invertibility of the matrix is a very important concept, which I will come back to in a bit.

Problems in direct inversion, misfits and gradient descents

Let’s see what happens when the first condition is not satisfied. I tried creating a non-linear function that at least I could not invert. Someone might be able to do it, but that is not the aim of this exercise. This new relationship is given here:

Complex flour to bread relationship, another 1D forward model.

The actual forward model is included in the script at the bottom of the post. When we cannot analytically solve the forward model, one approach would be to iteratively optimise the input parameters such that the end result becomes as close as possible to the observed real world. But, before we talk about how to propose inputs, we need a way to compare the outputs of our model to the observables!
When we talk about our bread, we see one end-product. It is very easy to compare the result of our simulation to the real world. But what happens when we make a few breads, or, similarly, when we have a lot of temperature observations in our volcanic model? One (common) way to combine observables and predictions into a single measure of discrepancy is by defining a misfit formula which takes into account all observations and predictions, and returns a single – you guessed it – misfit.
The choice of misfit directly influences which data is deemed more important in an inversion scheme, and many choices exist. One of the most intuitive is the L2-norm (or L2 misfit), which calculates the vector distance between the predicted data (from our forward model) with the observed data as if in Euclidean space. For our bread it would simply be

X = | predicted bread – observed bread |

Let’s try to create a misfit function for our 1D bread model. I use the relationship in the included scripts. Again, we have 500 grams of bread. By calculating the quantity | G(m) – 500 |, we can now make a graph of how the misfit varies as we change the amount of flour. I created a figure which shows the observed and predicted value:

Complex flour to bread relationship with observation superimposed

and a figure which shows the actual misfit at each value of m

Misfit between model and observation

Looking at the two previous figures may result in some questions and insights to the reader. First and foremost, it should be obvious that this problem does not have a unique solution. Different amounts of flour give exactly the same end result! It is thus impossible to say with certainty how much flour was used.
We have also recast our inversion problem as an optimisation function. Instead of thinking about fitting observations we can now think of our problem as a minimization of some function X (the misfit). This is very important in inversion.
Iterative optimizations schemes such as gradient descent methods (and the plethora of methods derived from it) work in the following fashion:

Gradient descent inversion/optimization

1. Pick a starting point based on the best prior knowledge (e.g., a 500 gram bread could have been logically baked using 500 grams of flour);
2. Calculate the misfit (X) at our initial guess;
3. If X is not low enough (compared to some arbitrary criterion):
• Compute the gradient of X;
• Do a step in the direction of the steepest gradient;
• Recompute X;
• Repeat from 2.
4. If X is low enough:
• Finalise optimisation, with the last point as solution.
These methods are heavily influenced by the starting point of the inversion scheme. If I would have started on the left side of the recipe-domain (e.g. 100 grams of flour), I might well have ended up in a different solution. Gradient descents often get ‘stuck’ in local solutions, which might not even be the optimal one! We will revisit this non-uniqueness in the 2D problem, and give some strategies to mitigate creating more than one solution. Extra material can be found here.

Simplistic gradient descent

Obvious solutions, grid searches, and higher dimensions

One thing that often bugged me about the aforementioned gradient descent methods is the seemingly complicated approach for such simple problems. Anyone looking at the figures could have said

Well, duh, there’s 3 solutions, here, here and here!

Why care about such an complicated way to only get one of them?
The important realisation to make here is that I have precomputed all possible solution for this forward model in the 0 – 700 grams range. This precomputation on a 1D domain was very simple; at a regular interval, compute the predicted value of baked bread. Following this, I could have also programmed my Python routine to extract all the values with a sufficiently low misfit as solutions. This is the basis of a grid search.
Let’s perform a grid search on our second model (1b). Let’s find all predictions with 500 grams of bread as the end result, plus-minus 50 grams. This is the result:

Grid search for the 2D model

The original 312.5 grams of flour as input is part of the solution set. However, the model actually has infinitely many solutions (extending beyond the range of the current grid search)! The reason that a grid search might not be effective is the inherent computational burden. When the forward model is sufficiently expensive in numerical simulation, exploring a model space completely with adequate resolution might take very long. This burden increases with model dimension; if more model parameters are present, the relevant model space to irrelevant model space becomes very small. This is known as the curse of dimensionality (very well explained in Tarantola’s textbook).
Another reason one might want to avoid grid searches is our inability to appropriately process the results. Performing a 5 dimensional or higher grid searches is sometimes possible on computational clusters, but visualizing and interpreting the resulting data is very hard for humans. This is partly why many supercomputing centers have in-house departments for data visualisation, as it is a very involved task to visualise complex data well.

Now: towards solving our physical inversions!

Non-uniqueness, regularization and linear algebra (bleh!)

One big problem in inversion is non-uniqueness: the same result can be obtained from different inputs. The go-to way to combat this is to add extra information of any form to the forward problem. In our bread recipe we could think of adding extra observables to our problem, such as the consistency of the bread, its taste, color, etc. Another option could be to add constraints on the parameters, such as using the minimum amount of ingredients. This is akin to asking the question: given this amount of bread, how much flour and water was minimally used to make it?
Diffusion type problems are notorious for their non-uniqueness. Many different subsurface heat conduction distributions might result in the observations (imagine differently inclined volcanic conduits). An often used method of regularisation (not limited to diffusion type studies!) is spatial smoothing. This method requires that among equally likely solutions, the smoothest solutions are favoured, for it is more ‘natural’ to have smoothly varying parameters. Of course, in many geoscientific settings one would definitely expect sharp contrasts. However, in ‘underdetermined’ problems (i.e., you do not have enough observations to constrain a unique solution), we favour Occam’s Razor and say

The simplest solution must be assumed

When dealing with more parameters than observables (non-uniqueness) in linear models it is interesting to regard the forward problem again. If one would parameterize our volcanic model using 9 parameters for the subsurface and combine that with the 3 measurements at the surface, the result would be an underdetermined inverse problem.

Rough parametrization for the heat conduction model

This forward model (the Laplace equation) can be discretised by using, for example, finite differences. The resulting matrix equation would be Am = d, with A a 3 x 9 matrix, m a 9 dimensional vector and d a 3 dimensional vector. As one might recall from linear algebra classes, for a matrix to have an inverse, it has to be square. This matrix system is not square, and therefore not invertible!

Aaaaahhh! But don’t panic: there is a solution

By adding either prior information on the parameters, smoothing, or extra datapoints (e.g., taking extra measurements in wells) we can make the 3 x 9 system a perfect 9 x 9 system. By doing this, we condition our system such that it is invertible. However, many times we end up overdetermining our system which could result in a 20 x 9 system, for example. Note that although neither the underdetermined nor the overdetermined systems have an exact matrix inverse, both do have pseudo-inverses. For underdetermined systems, I have not found these to be particularly helpful (but some geophysicists do consider them). Overdetermined matrix systems on the other hand have a very interesting pseudo-inverse: the least squares solution. Finding the least squares solution in linear problems is the same as minimising the L2 norm! Here, two views on inversion come together: solving a specific matrix equation is the same as minimising some objective functional (at least in the linear case). Other concepts from linear algebra play important roles in linear and weakly non-linear inversions. For example, matrix decompositions offer information on how a system is probed with available data, and may provide insights on experimental geophysical survey design to optimise coverage (see “Theory of Model-Based Geophysical Survey and Experimental Design Part A – Linear Problems” by Andrew Curtis).
I would say it is common practice for many geophysicists to pose an inverse problem that is typically underdetermined, and keep adding regularization until the problem is solvable in terms of matrix inversions. I do not necessarily advocate such an approach, but it has its advantages towards more agnostic approaches, as we will see in the post on probabilistic inversion next week!


Summary of deterministic inversions

We’ve seen how the forward model determines our inversion problem, and how many measurements can be combined into a single measure of fit (the misfit). Up to now, three inversion strategies have been introduced:

Direct inversion: analytically find a solution to the forward problem. This method is limited to very specific simple cases, but of course yields near perfect results.
Gradient descent methods: a very widely used class of algorithms that iteratively update solutions based on derivatives of the misfit function. Their drawbacks are mostly getting stuck in local minima, and medium computational cost.
Grid searches: a method that searches the entire parameter space systematically. Although they can map all the features of the inverse problem (by design), they are often much too computationally expensive.

What might be even more important, is that we have seen how to reduce the amount of possible solutions from infinitely many to at least a tractable amount using regularisation. There is only one fundamental piece still missing… Stay tuned for the last blog post in this series for the reveal of this mysterious missing ingredient!

Interested in playing around with inversion yourself? You can find a toy code about baking bread here.