Computational Methods for Earth and Planetary Science
The Computational Methods for Earth and Planetary Science Group at the Niels Bohr Institute is engaged in inverse theory and algorithm development for the solution of complex problems in the study of the Earth's interior. Of particular interest to us is the integration of quantitative geophysical and geological models and the related computational challenges.
Our projects span from large-scale planetary studies to Earth resource exploration, including drinking water, geothermal energy and mineral resources. All these activities share common themes that continue to fascinate us: nonlinear inverse problems, the search for feasible solutions, probabilistic Earth models, and the complex flow of information from uncertain observations, through numerical modeling, to final interpretions and decisions.
Our research is about the fundamental assumptions behind data analysis algorithms, development of new computational ideas, and new ways of tackling the most difficult challenges in computational Earth and planetary science.
Our methods are widely used in the study of Earth systems, the Lunar surface and interior, Mars' interior, and Exoplanetary systems.
Compared to the controlled laboratory measurements in many fields of physics and chemistry, Earth and planetary data are sparse, noisy, biased, and have limited resolution. Yet, this data is of key importance for understanding the complexity of the macroscopic world that surrounds us.
From the core of the Earth, to the surfaces and interiors of terrestrial planets and moons, and further out to distant exoplanetary systems, we strive to understand inaccurate, insufficient and inconsistent data, using the laws of physics. This endeavor relies on computational methods whose reliability and uncertainties are not always well understood.
Informed proposal Monte Carlo
Any search or sampling algorithm for solution of inverse problems needs guidance to be efficient. Many algorithms collect and apply information about the problem on the fly, and much improvement has been made in this way.
However, as a consequence of the No-Free-Lunch Theorem, the only way we can ensure a significantly better performance of search and sampling algorithms is to build in as much external information about the problem as possible.
In the special case of Markov Chain Monte Carlo (MCMC) sampling we review how this is done through the choice of proposal distribution, and we show how this way of adding more information about the problem can be made particularly efficient when based on an approximate physics model of the problem.
A highly non-linear inverse scattering problem with a high-dimensional model space serves as an illustration of the gain of efficiency through this approach.
- Sarouyeh Khoshkholgh, Niels Bohr Institute, University of Copenhagen
- Khoshkholgh, S., Zunino, A. and Mosegaard, K., 2022. Full-waveform Inversion by Informed-Proposal Monte Carlo. Geophysical Journal International, 230, 1824–1833.
- Khoshkholgh, S., Zunino, A. and Mosegaard, K., 2021. Informed Proposal Monte Carlo. Geophysical Journal International 226 (2), 1239-1248.
Constraining rates in Earth’s sedimentary archives from their genetic history.
Sedimentary rock piles are important archives of paleoclimatic information and also serve as CO2 storage sites in the near future.
The deposition rate of sediments is of key importance to all sedimentological studies, as it reveals the sedimentary properties and preserve paleoclimatic history.
X- Ray-Fluorescence spectroscopy (XRF) provide a wealth of new data with ~30 chemical elements measured at ultra-high stratigraphic resolution (>0.1 mm) that can even preserve ‘annual layers’ in sediments of essentially any age and teach us about paleoclimatic events at time scales relevant for humans.
We have discovered periodic ‘Milankovitch’ signals dictated by solar insolation forcing on the early Earth in multiple chemical elements detected in the XRF signals, but existing signal analysis is only applicable to unidimensional data.
By developing a Bayesian inverse approach, we will take advantage of the multidimensionality of XRF data that can better constrain sedimentation rate because
- numerous elements oscillate independently in concert with climate-driven Milankovitch forcing and
- the chemical composition reveals sedimentation rates; e.g., coarse grained material (e.g. sand) deposits faster than fine-grained material (e.g. mud).
- deploy unsupervised learning methods (Self Organizing Maps; SOM) to decipher the underlying formative processes and assign sedimentation rates from XRF data and
- develop a Bayesian inversion approach to quantitatively link astronomical theory with multidimensional geochemical data.
Together, we expect these algorithms will enable rate determination in a wider range of Earth’s sedimentary archives and not only refine the geological time scale, but also set a new standard for the mapping of Earth’s sedimentary archives.
- Tais Wittchen Dahl, Globe Institute, University of Copenhagen
Mars Insight Project: Seismic studies of the Martian crust
The estimation of crustal structure and thickness is essential in understanding the formation and evolution of terrestrial planets and moons.
Initial planetary missions with seismic instrumentation on board face the additional challenge of dealing with seismic activity levels that are only poorly constrained a priori.
For example, the lack of plate tectonics on Mars leads to low seismicity, which could, in turn, hinder the application of many terrestrial data analysis techniques.
We propose using a joint inversion of receiver functions and apparent incidence angles, which contain information on absolute S-wave velocities of the subsurface.
Since receiver function inversions suffer from a velocity depth trade-off, we in addition exploit a simple relation that defines apparent S-wave velocity as a function of observed apparent P-wave incidence angles to constrain the parameter space.
We then use the Neighborhood Algorithm for the inversion of a suitable joint objective function.
The resulting ensemble of models is then resampled and used to derive unbiased uncertainty estimates for each model parameter.
In preparation for the analysis of data from the InSight mission, we show the application of our proposed method on Mars synthetics and sparse terrestrial data sets from different geological settings using both single and multiple events.
We use information-theoretic statistical tests as model selection criteria and discuss their relevance and implications in a seismological framework.
- Rakshit Joshi, Max-Planck-Institute for Solar System Research, Göttingen, and Ludwig-Maximilians-Universität, Munich
- Brigitte Knapmeyer-Endrun, Bensberg Observatory, University of Cologne
- Heiner Igel, Ludwig-Maximilians-Universität, Munich
- Ulrich R. Christensen, Max-Planck-Institute for Solar System Research, Göttingen
Joshi, R., Knapmeyer-Endrun, B., Mosegaard, K., Igel, H., and Christensen, U., 2021. Joint inversion of receiver functions and apparent incidence angles for sparse seismic data. Earth and Space Science. DOI: 10.1029/2021EA001733.
Extraterrestrial seismology saw its advent with the deployment of seismometers during the Apollo missions that were undertaken from July 1969 to December 1972.
The Apollo lunar seismic data constitute a unique resource being the only seismic data set which can be used to infer the interior structure of a planetary body besides the Earth.
Ongoing analysis and interpretation of the seismic data continues to provide constraints that help refine lunar origin and evolution.
In addition to this, lateral variations in crustal thickness (~0–80 km) are being mapped out at increasing resolution from gravity and topography data that have and continue to be collected with a series of recent lunar orbiter missions.
Many of these also carry onboard multi-spectral imaging equipment that is able to map out major-element concentration and surface mineralogy to high precision.
These results coupled with improved laboratory-based petrological studies of lunar samples provide important constraints on models for lunar magma ocean evolution, which ultimately determines internal structure.
Whereas existing constraints on initial depth of melting and differentiation from quantitative modeling suggested only partial Moon involvement (<500 km depth), more recent models tend to favor a completely molten Moon, although the former cannot be ruled out sensu stricto.
Recent geophysical analysis coupled with thermodynamical computations of phase equilibria and physical properties of mantle minerals suggest that the Earth and Moon are compositionally distinct.
Continued analysis of ground-based laser ranging data and recent discovery of possible core reflected phases in the Apollo lunar seismic data strengthens the case for a small dense lunar core with a radius of <400 km corresponding to 1–3% of lunar mass.
An inquiry into the lunar interior: A nonlinear inversion of the Apollo lunar seismic data
This study discusses in detail the inversion of the Apollo lunar seismic data and the question of how to analyze the results. The well-known problem of estimating structural parameters (seismic velocities) and other parameters crucial to an understanding of a planetary body from a set of arrival times is strongly nonlinear.
Here we consider this problem from the point of view of Bayesian statistics using a Markov chain Monte Carlo method. Generally, the results seem to indicate a somewhat thinner crust with a thickness around 45 km as well as a more detailed lunar velocity structure, especially in the middle mantle, than obtained in earlier studies.
Concerning the moonquake locations, the shallow moonquakes are found in the depth range 50–220 km, and the majority of deep moonquakes are concentrated in the depth range 850–1000 km, with what seems to be an apparently rather sharp lower boundary.
In wanting to further analyze the outcome of the inversion for specific features in a statistical fashion, we have used credible intervals, two-dimensional marginals, and Bayesian hypothesis testing.
Using this form of hypothesis testing, we are able to decide between the relative importance of any two hypotheses given data, prior information, and the physical laws that govern the relationship between model and data, such as having to decide between a thin crust of 45 km and a thick crust as implied by the generally assumed value of 60 km.
We obtain a Bayes factor of 4.2, implying that a thinner crust is strongly favored.
- Amir Khan, Dep. of Earth Sciences, ETH Zürich
- A Pommier, School of Earth and Space Exploration, Arizona State University, Tempe, USA
- GA Neumann, NASA Goddard Space Flight Center, Greenbelt, MD, USA
- A Khan, A Pommier, GA Neumann, K Mosegaard 2013. The lunar moho and the internal structure of the Moon: A geophysical perspective: Tectonophysics 609, 331-352
- A Khan, K Mosegaard 2002. An inquiry into the lunar interior: A nonlinear inversion of the Apollo lunar seismic data: Journal of Geophysical Research: Planets 107 (E6), 3-1-3-23.
High-resolution topography from planetary images and laser altimetry
Mapping landforms on the Moon is of great interest and importance for future human settlements and resources exploration.
One of the first steps is to map the topography in great detail and resolution. However, data from the Lunar Orbiter Laser Altimeter (LOLA) provide low-resolution elevation maps in comparison to the size of detailed geological features.
To improve resolution, we developed a new method to upscale topographic maps to a higher resolution using images from the Lunar Reconnaissance Orbiter Camera (LROC).
Our method exploits the relation between topographic gradients and degrees of shading of incoming sunlight.
In contrast to earlier published methods, our approach is based on probabilistic, linear inverse theory, and its computational efficiency is very high due to its formulation through the Sylvester Equation.
The method operates on multiple images and incorporates albedo variations. A further advantage of the method is that we avoid/reduce the use of arbitrary tuning parameters through a probabilistic formulation where all weighting of data and model parameters is based on prior information about data uncertainties and reasonable bounds on the model.
Our results increase the resolution of the topography from ∼60 m per pixel to 0.9 m per pixel, bringing it to the same pixel resolution as the optical images from LROC, allowing in some cases detection of craters as small as ∼3 m of diameter.
We estimate uncertainties of the topographic model due to noise in the images, and in the low-resolution (LOLA) model.
- Iris Fernandes, Niels Bohr Institute, University of Copenhagen
- Fernandes, I., and Mosegaard K., 2022. High-Resolution topography from planetary images and laser altimetry: Planetary and Space Science, Vol. 218, 2022, 105514, https://doi.org/10.1016/j.pss.2022.105514
GeoAfrica: Utilization of geothermal resources in East Africa
The project aims at providing the scientific, technological and policy basis required for the widespread utilization of geothermal resources in East Africa.
We provide the basis to evaluate the potential of low- and medium-temperature geothermal resources, and develop technological solutions for their use, providing sustainable and reliable energy (electricity and heating/cooling). Access to green, cheap energy will support the sustainable growth and decarbonization of the economy.
The proposed solution for low- and medium-temperature reservoirs consists of a novel integration of multi-generation energy technologies that maximize conversion efficiencies while minimizing costs.
We explore the potential of selected reservoirs through integrated geoscientific studies, incorporating drilling, field and production data. A prototype plant at the Olkaria field (Kenya) will be constructed.
- Ivanka Orozova-Bekkevold, Niels Bohr Institute, University of Copenhagen
Studying the rocky interior of planets in a multiplanetary system
Interior characterization traditionally relies on individual planetary properties, ignoring correlations between different planets of the same system. For multiplanetary systems, planetary data are generally correlated.
This is because the differential masses and radii are better constrained than absolute planetary masses and radii.
We explore such correlations and data specific to the multiplanetary system of TRAPPIST-1 and study their value for our understanding of planet interiors.
Furthermore, we demonstrate that the rocky interior of planets in a multiplanetary system can be preferentially probed by studying the densest planet representing a rocky interior analog.
Our methodology includes a Bayesian inference analysis that uses a Markov chain Monte Carlo scheme.
Our interior estimates account for the anticipated variability in the compositions and layer thicknesses of core, mantle, water oceans, and ice layers, as well as a gas envelope.
Our results show that
- interior estimates significantly depend on available abundance proxies and
- the importance of interdependent planetary data for interior characterization is comparable to changes in data precision by 30%.
For the interiors of TRAPPIST-1 planets, we find that possible water mass fractions generally range from 0% to 25%. The lack of a clear trend of water budgets with orbital period or planet mass challenges possible formation scenarios.
While our estimates change relatively little with data precision, they critically depend on data accuracy. If planetary masses varied within ±24%, interiors would be consistent with uniform (∼7%) or an increasing water mass fractions with orbital period (∼2%–12%).
- Caroline Dorn and Simon L. Grimm, University of Zurich
- Yann Alibert, University of Bern
- Caroline Dorn, Klaus Mosegaard, Simon L. Grimm, and Yann Alibert 2018. Interior Characterization in Multiplanetary Systems: TRAPPIST-1: The Astrophysical Journal, 865:20, 2018 September 20
Klaus Mosegaard, Professor
Tagensvej 16, 2200 København N.
Phone: +45 21664566
|Ask Hammer||MSc. Student|
|Dan Friis Tømmerby Jensen||MSc. Student|