1 - Inverse Problems and Data Assimilation in Earth Sciences

Summary

Abstract: We introduce direct and inverse problems, which describe dynamical processes causing change in the Earth system and its space environment. A well-posedness of the problems is defined in the sense of Hadamard and in the sense of Tikhonov, and it is linked to the existence, uniqueness, and stability of the problem solution. Some examples of ill- and well-posed problems are considered. Basic knowledge and approaches in data assimilation and solving inverse problems are discussed along with errors and uncertainties in data and model parameters as well as sensitivities of model results. Finally, we briefly review the book’s chapters which present state-of-the-art knowledge in data assimilation and geophysical inversions and applications in many disciplines of the Earth sciences: from the Earth’s core to the near-Earth environment.

1.1 Introduction

Many problems in Earth sciences are related to dynamic processes within the planet, on its surface, and in its space environment. Geoscientists study the processes using observations and measurements of their manifestations. Each process can be presented by a model described by physical and/or chemical laws and a set of relevant parameters. The model, in its turn, can be represented by a mathematical model; that is, a set of partial differential equations or ordinary differential equations with boundary and/or initial conditions defined in a specific domain. The mathematical model links its parameters and variables with a set of data from observations and measurements and provides a connection between the causal characteristics of the dynamic process and its effects. The causal characteristics include, for example, physical parameters (such as velocity, temperature, pressure), parameters of the initial and boundary conditions, and geometrical parameters of a model domain.

The aim of the direct mathematical problem is to determine the effects of a dynamic model process based on the knowledge of its causes, and hence to find a solution to the mathematical problem for a given set of model parameters. An inverse problem is the opposite of a direct problem. An inverse problem is considered when there is a lack of information on the causal characteristics but information on the effects of the dynamic process exists (e.g. Kirsch, 1996; Kabanikhin, 2011; Ismail-Zadeh et al., 2016). For example, the seismic wave velocities inferred from seismograph’s measurements on the Earth’s surface are related to a fault rupture in the lithosphere; the rupture process and the wave propagation are described mathematically by the wave equation, which relates causal characteristics (the velocity, density, and elastic properties) with their effects. The heat flow measured at the Earth’s surface and inferred temperature can be related to the heat equation linking temperature, thermal conductivity, density, and specific heat.

For centuries, physicists searched for and discovered the causes of the effects of the geophysical processes they observed, so that they solved simplified inverse problems. Inverse problems, as formalised mathematical studies, have been initiated in the twentieth century (e.g. Weyl, 1911). These problems allow for determining model parameters or specific model conditions that cannot be directly observed. Inverse problems can be subdivided into time-reverse or retrospective problems (e.g. to determine initial conditions in the past and/or to restore the development of a dynamic process), coefficient problems (e.g. to determine the coefficients of the model and/or boundary conditions), geometrical problems (e.g. to determine the location of heat sources in a model domain or the geometry of the model boundary), and some others.

1.2 Inverse Problems and Well-Posedness

The idea of well- (and ill-)posedness in the theory of partial differential equations was introduced by Hadamard (1902). A mathematical model of a geophysical problem is well-posed if (i) a solution to this problem exists, and the solution is (ii) unique and (iii) stable. Problems for which at least one of these three properties does not hold are called ill-posed. Existence of the problem’s solution is normally proven by mathematicians, at least in the simplest cases. Meanwhile, if the solution exists, it may not be unique, allowing for multiple theoretical solutions, as in the case of potential-field interpretation. The non-uniqueness of potential-field studies is associated with the general topic of scientific uncertainty in the Earth sciences, because problems are generally addressed with incomplete and imprecise data (Saltus and Blakely, 2011). The requirement of stability is the most important in numerical modelling. If a problem lacks the property of stability, then its solution is almost impossible to compute, because computations are polluted by unavoidable errors. If the solution of a problem does not depend continuously on the initial data, then, in general, the computed solution may have nothing to do with the true solution.

Inverse problems are often ill-posed. For example, the retrospective (inverse) problem of thermal convection is an ill-posed problem, since the backward heat problem, describing both heat advection and conduction backwards in time, possesses the property of instability (e.g. Kirsch, 1996). In particular, the solution to the problem does not depend continuously on the initial data. This means that small changes in the present-day temperature field may result in large changes of predicted temperatures in the past. Following Ismail-Zadeh et al. (2016), this statement is explained using a simple problem related to the one-dimensional (1-D) diffusion equation (Example 1).

Example 1

Example 1 Consider the following problem for the 1-D backward diffusion equation

$\partial u(t,x)/\partial t={\partial }^{2}u(t,x)/\partial x^2,0\le x\le \pi ,t<\mathrm{0,}$(1.1)

with the following boundary and initial conditions

$u(t,0)=0=u(t,\pi ),t\le 0,u(0,x)={\varphi }_n\left(x\right),0\le x\le \pi .$(1.2)

At the initial time, the function ${\varphi }_n\left(x\right)$ is assumed to take the following two forms:

${\varphi }_n\left(x\right)=\frac{\sin \left(\right(4n+1\left)x\right)}{4n+1}\text{and}{\varphi }_0\left(x\right)\equiv 0.$(1.3)

We note that

$\underset{0\le x\le \pi }{\max }\left|{\varphi }_n\left(x\right)-{\varphi }_0\left(x\right)\right|\le \frac{1}{4n+1}\to 0\text{ at }n\to \infty .$(1.4)

The two solutions of the problem

$u_n(t,x)=\frac{\sin \left(\right(4n+1\left)x\right)}{4n+1}\exp (-{(4n+1)}^{2}t)\text{ at }{\varphi }_n\left(x\right)={\varphi }_n\text{and}$(1.5)

$u_0(t,x)\equiv 0\text{ at }{\varphi }_n\left(x\right)={\varphi }_0$(1.6)

correspond to the two chosen functions of ${\varphi }_n\left(x\right)$ , respectively. At $t=-1$ and $x=\pi /2$

$u_n(-1,\frac{\pi }{2})-u_0(-1,\frac{\pi }{2})=\frac{1}{4n+1}\exp \left({(4n+1)}^2\right)\to \infty \text{ at }n\to \infty .$(1.7)

At large n, two closely set initial functions ${\varphi }_n$ and ${\varphi }_0$ are associated with the two strongly different solutions at $t=-1$ and $x=\pi /2$ . Hence, a small error in the initial data (1.4) can result in very large errors in the solution to the backward problem (1.7), and therefore the solution is unstable, and the problem is ill-posed in the sense of Hadamard.

Although inverse problems are quite often unstable and hence ill-posed, there are some methods for solving them. The idea of conditionally well-posed problems and the regularisation method were introduced by Tikhonov (1963). According to Tikhonov, a class of admissible solutions to conditionally ill-posed problems should be selected to satisfy the following conditions: (i) a solution exists in this class, (ii) the solution is unique in the same class, and (iii) the solution depends continuously on the input data (i.e. the solution is stable). The Tikhonov regularisation is essentially a trade-off between fitting the observations and reducing a norm of the solution to the mathematical model of a geophysical problem. We show the differences between the Hadamard’s and Tikhonov’s approaches to ill-posed problems in Example 2.

Example 2

Example 2 Consider the problem for the 1-D backward diffusion equation (like the problem presented in Example 1)

$\partial u(t,x)/\partial t={\partial }^{2}u(t,x)/\partial x^2,0\le x\le \pi ,-T\le t<0,$(1.11)

with the boundary and initial conditions (1.2). The solution of the problem satisfies the inequality

$\Vert u(t,x)\Vert \le {\Vert u(T,x)\Vert }^{-t/T}{\Vert u(0,x)\Vert }^{1+t/T},$(1.12)

where the norm is presented as ${\Vert u(t,x)\Vert }^2\equiv \underset{0}{\overset{\pi }{\int }}{u}^2(t,x)dx$ . We note that the inequality

$\Vert u(t,x)\Vert \le M^{-t/T}{\Vert u_0\Vert }^{1+t/T}$(1.13)

is valid in the class of functions $\Vert u(t,x)\Vert \le M=const$ (Samarskii and Vabischevich, 2007). Inequality (1.13) yields a continuous dependence of the problem’s solution on the initial conditions, and hence to well-posedness of the problem in the sense of Tikhonov. Therefore, the Tikhonov approach allows for developing methods for regularisation of the numerical solution of unstable problems.

1.3 Data Assimilation

With a growth of data related to Earth observations and laboratory measurements, the enhancement of data quality and instrumental accuracy, as well as the sophistication of mathematical and numerical models, the assimilation of available information into the models to determine specific states of geophysical/geochemical dynamic processes as accurately as possible becomes an essential tool in solving inverse problems. Data assimilation can be defined as the incorporation of observations and initial/boundary conditions in an explicit dynamic model to provide time continuity and coupling among the physical characteristics of the dynamic model (e.g. Kalnay, 2003; Lahoz et al., 2010; Law et al., 2015; Asch et al., 2016; Ismail-Zadeh et al., 2016; Fletcher, 2017).

Data assimilation has been pioneered by meteorologists and successfully applied to improve operational weather forecasts (e.g. Ghil and Malanotte-Rizzoli, 1991; Kalnay, 2003). To produce forecasts, initial conditions are required for the weather prediction models resembling the current state of the atmosphere. Data assimilation starts with ‘unknown’ forecasts and applies corrections to the forecast based on observations and estimated errors in the observations and in the forecasts. The difference between the forecast and the observed data at a certain time is assessed using different methods to make new forecast to better fit the observations.

Data assimilation and geophysical inversions have also been widely used in oceanography (e.g. Ghil and Malanotte-Rizzoli, 1991; Bennett, 1992), hydrology (e.g. McLaughlin, 2002), seismology (e.g. Backus and Gilbert, 1968), geodynamics (e.g. Bunge et al., 2003; Ismail-Zadeh et al., 2003; 2004; 2016), geomagnetism (Fournier et al., 2007; Liu et al., 2007), and other Earth science disciplines (e.g. Park and Xu, 2009; Lahoz et al., 2010; Blayo et al., 2014). We note, that depending on the geoscience discipline, data assimilation is also referred to as state estimation, history matching, and data-driven analysis.

1.4 Data Assimilation and Inversions: Basic Approaches and Sensitivity Analysis

Part I of the book introduces basic knowledge and approaches in data assimilation and inversions and presents a high-order sensitivity analysis to obtain best estimate results with reduced uncertainties.

There are two basic approaches to solve inverse problems: classical and Bayesian. The classical approach considers a mathematical model as a true model describing the physical process under study, and geoscientific data as the only available data set with some measurement errors. The goal of this approach is to recover the true model (e.g. initial or boundary conditions). Another way to treat a mathematical model is the Bayesian approach, where the model is considered as a random variable, and the solution is a probability distribution for the model parameters (e.g. Aster et al., 2005).

Chapter 2 discusses these approaches in more detail. This chapter also provides an accessible general introduction to the breadth of geophysical inversions and presents similarities and connections between different approaches (Valentine and Sambridge, this volume). Chapter 3 introduces the Bayesian data assimilation providing a history of geophysical data assimilation and its current directions (Grudzien and Bocquet, this volume).

All variables in data assimilation models (e.g. state variables describing physical properties, such as velocity, pressure, or temperature; initial and/or boundary conditions; and parameters such as viscosity or thermal diffusivity) can be polluted by errors. The source of errors comes from imperfect measurements and computations. Experimental or calibration standard errors result in measurement errors. Systematic errors in numerical modelling are associated with a mathematical model, its discretisation, and iteration errors. Model errors are associated with the idealisation of Earth system dynamics by a set of conservation equations governing the dynamics. Model errors can be defined as the difference between the actual Earth system dynamics and the exact solution of the mathematical model. Discretisation errors are associated with the difference between the exact solution of the conservation equations and the exact solution of the algebraic system of equations obtained by discretising these equations. And iteration errors are defined as the difference between the iterative and exact solutions of the algebraic system of equations (Ismail-Zadeh and Tackley, 2010). Also, errors can stem from imperfectly known physical processes or geometry. Determining the changes in computed model responses that are induced by variations in the model parameters (e.g. due to errors) is the scope of sensitivity analysis, which is linked to the stability of systems to small errors.

Sensitivity analysis assists in understanding the stability of the model solution to small perturbations in input variables or parameters. For instance, consider the thermal convection in the Earth’s mantle. If the temperature in the geological past is determined from the solution of the backward thermal convection problem using present mantle temperature assimilated to the past, the following question arises: what will be the temperature variation due to small perturbations of the present temperature data? The gradient of the objective functional (representing the misfit between the model and measured temperature) with respect to the present temperature in a variational data assimilation gives the first-order sensitivity coefficients (e.g. Hier-Majumder et al., 2006). A second-order adjoint sensitivity analysis presents some challenges associated with cumbersome computations of the product of the Hessian matrix of the objective functional with a vector (Le Dimet et al., 2002).

Chapter 4 discusses higher-order sensitivity and uncertainty analysis to obtain best estimates with reduced uncertainties. The analysis is applied to an inverse radiation transmission problem, to an oscillatory dynamical model, and to a large-scale computational model involving thousands of uncertain parameters. The examples illustrate impacts of the first-, second-, and third-order response sensitivities to parameters on the expectation, variance, and skewness of the respective model responses (Cacuci, this volume).

1.5 Data Assimilation and Inverse Problems in ‘Fluid’ Earth Sciences

Part II of the book is dedicated to applications of data assimilation and inversions to problems related to the cryosphere, hydrosphere, atmosphere, and near-Earth environment (‘fluid’ Earth spheres).

Estimates of seasonal snow often bear significant uncertainties (e.g. due to error-prone forcing data and parameter estimates), and data assimilation becomes a useful tool to minimise inherent limitations that result from the uncertainty. Chapter 5 reviews current snow models, snow remote sensing methods, and data assimilation techniques that can reduce uncertainties in the characterisation of seasonal snow (Girotto et al., this volume). Although some properties at the surface of glaciers and ice sheets can be measured from remote sensing or in-situ observations, other characteristics, such as englacial and basal properties, or past climate conditions, remain difficult or impossible to observe. Data assimilation in glaciology assists in inferring unknown properties and boundary conditions to be employed in numerical models (Morlighem and Goldberg, this volume). Chapter 6 presents common applications of data assimilation in glaciology, and some of the new directions that are currently being developed.

Data assimilation in many hydrological problems shares distinct peculiarities: scarce or indirect observation of important state variables (e.g. soil moisture, river discharge, groundwater level), very incomplete or largely conceptual modelling, extreme spatial heterogeneity, and uncertainty on controlling physical parameters (Castelli, this volume). Chapter 7 discusses the peculiarities of data assimilation for state estimation and model inversion in hydrology related to the following applications: soil moisture, runoff for flood and inundation prediction, static geophysical inversion in groundwater modelling, and dynamic geophysical inversion in coupled surface water and energy balance.

Robust estimates of trace gas emissions that impact air quality and climate provide important knowledge for informed decision-making. Better monitoring and increasing data availability due to expanding observing networks provide information on the changing composition of the atmosphere. Chapter 8 discusses the use of various inverse modelling approaches to quantify emissions of environmentally important trace gases, with a focus on the use of satellite observations. It presents the inverse problem of retrieving the atmospheric trace gas information from the satellite measurements, and the subsequent use of these satellite data for inverse modelling of sources and sinks of the trace gases (Jones, this volume).

Models of volcanic cloud propagation due to volcanic eruptions assist in operational forecasts and provide invaluable information for civil protection agencies and aviation authorities during volcanic crises. Quantitative operational forecasts are challenging due to the large uncertainties that typically exist when characterising volcanic emissions in real time, and data assimilation assists in reduction of quantitative forecast errors (Folch and Mingari, this volume). Chapter 9 reviews state-of-the-art in data assimilation of volcanic clouds and its use in operational forecasts.

Energetic charged particles trapped by the Earth magnetic field present a significant hazard for Earth orbiting satellites and humans in space, and data assimilation helps to reconstruct the global state of the radiation particle environment from observations (Shprits et al. this volume). Chapter 10 describes recent studies related to data assimilation in the near-Earth electron radiation environment. Applications to the reanalysis of the radiation belts and ring current, real-time predictions, and analysis of the missing physical processes are discussed in the chapter.

1.6 Data Assimilation and Inverse Problems in Solid Earth Sciences

Part III presents methods and applications of data assimilation and inversions in problems of the solid Earth sciences: geochronology, volcanology, seismology, gravity, geodesy, geodynamics, and geomagnetism.

Chapter 11 presents applications of inverse methods, namely, trans-dimensional Markov chain Monte Carlo, to geochronological and thermochronological data to identify the number of potential source components for detrital material in sedimentary basins and to extract temperature histories of rocks over geological time (Gallagher, this volume).

Lava dome growth and lava flow are two main manifestations of effusive volcanic eruptions. Chapter 12 discusses inverse problems related to lava dynamics. One problem is related to a determination of the thermal state of a lava flow from thermal measurements at its surface using a variational data assimilation method. Another problem aims to determine magma viscosity by comparison of observed and simulated lava domes employing artificial intelligence methods (Ismail-Zadeh et al., this volume).

Chapter 13 deals with data assimilation for real-time shake-mapping and ground shaking forecasts to assist in earthquake early warning systems. The current seismic wavefield is rapidly estimated using data assimilation, and then the future wavefield is predicted based on the physics of wave propagation (Hoshiba, this volume).

Global seismic tomography using time domain waveform inversion is overviewed in Chapter 14 in the context of imaging the Earth’s whole mantle at the global scale. The chapter discusses how the tomography problem is addressed, data selection approaches, definitions of the misfit function, and computation of kernels for the inverse step of the imaging procedure, as well as the choice of optimisation method (Romanowicz, this volume). The diversity of seismic inverse problems – in terms of scientific scope, spatial scale, nature of the data, and available resources – precludes the existence of a silver bullet to solve the problems. Chapter 15 describes smart methods for solving the inverse problems, which increase computational efficiency and usable data volumes, sometimes by orders of magnitude (Gebraad et al., this volume).

Chapter 16 deals with joint inversions as a hypothesis testing tool to study the Earth’s subsurface. It presents an application of joint inversions of gravity and magnetic data with seismic constraints in the western United States. As a result of the joint inversions, high velocity structures in the crust are found to be associated with relatively low-density anomalies, potentially indicating the presence of melt in a rock matrix (Moorkamp, this volume). An application of gravity inversion of Bouguer anomalies is presented in Chapter 17 focusing on the Moho depth and crustal density structure in the Tibetan Plateau. The inversion results clearly recognise a thick Tibetan crust and Moho depths of more than 60 km (Jin and Xuan, this volume).

Chapter 18 describes geodetic inversions and applications in geodynamics. Rapid development of the Global Navigation Satellite Systems (GNSS) allows enhanced geodynamic studies providing information about global-scale plate motions, plate boundary deformation, seismotectonic deformation, volcanology, postglacial isostatic rebound, ice flow, and water mass flow. A geophysical interpretation of GNSS observations is based on rheological models used to predict surface motions related to various tectonic processes and the corresponding inversion technique permitting us to separate the processes and to evaluate their parameters (Steblov and Vladimirova, this volume).

In Chapter 19, basic methods for data assimilation used in geodynamic modelling are described: backward advection method, variational (adjoint) method, and quasi-reversibility method. To demonstrate the applicability of the methods, two models are considered: a model of restoring prominent mantle plumes from their diffused stage, and a model of reconstruction of the thermal state of the mantle beneath the Japanese islands and their surroundings during forty million years. Also, this chapter discusses challenges, advantages, and disadvantages of the data assimilation methods (Ismail-Zadeh, Tsepelev, and Korotkii, this volume). Chapter 20 deals with global mantle convection in the Earth. Variational data assimilation allows for retrodicting past states of the Earth’s mantle as optimal flow histories relative to the current state. Poorly known mantle flow parameters, such as rheology and composition, can be then tested against observations extracting information from the geologic records (Bunge et al., this volume).

Chapter 21 presents geomagnetic data assimilation, which aims to optimally combine geomagnetic observations and numerical geodynamo models to better estimate the dynamic state of the Earth’s outer core and to predict geomagnetic secular variation. It provides a comprehensive overview of recent advances in the field, as well as some of the immediate challenges of geomagnetic data assimilation, possible solutions, and pathways to move forward (Kuang et al., this volume). Chapter 22 introduces main characteristics of geomagnetic data and magnetic field models and explores the role of model and observation covariances and localisation in typical assimilation setups, focusing on the use of three-dimensional dynamo simulations as the background model (Sanchez, this volume).

Conclusion

Inverse problems and data assimilation in Earth sciences provide many benefits to science and society. Mathematical and numerical models and methods involved in solving inverse problems and in assimilating data assist in utilisation of Earth observations and add value to the observations, for example, providing insight into physical/chemical processes and their observed manifestations. At the same time, observed and measured data help to constrain and sophisticate models and hence provide more reliable model estimates and forecasts. Society benefits from the knowledge obtained by using scientific products such as weather, air quality, space weather, and other forecasts. Applications of inverse problems and data assimilation in Earth sciences are broad, and this book covers only a part of them, including applications in atmospheric, cryospheric, geochronological, geodetical, geomagnetic, hydrological, seismological, and volcanological sciences. It presents the basics of modern theory and how theoretical methods works to decipher the puzzles of nature.