# Quasi-reversibility method for data assimilation in models of mantle dynamics

A. Ismail-Zadeh, A. Korotkii, G. Schubert, I. Tsepelev

► To cite this version:

A. Ismail-Zadeh, A. Korotkii, G. Schubert, I. Tsepelev. Quasi-reversibility method for data assimilation in models of mantle dynamics. Geophysical Journal International, 2007, 170, pp.1381-1398. 10.1111/j.1365-246X.2007.03496.x. insu-03603179

HAL Id: insu-03603179

https://insu.hal.science/insu-03603179

Submitted on 9 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quasi-reversibility method for data assimilation in models of mantle dynamics

A. Ismail-Zadeh,<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> A. Korotkii,<fnref n="4" /> G. Schubert<fnref n="5" /> and I. Tsepelev<fnref n="4" />

Accepted 2007 May 15. Received 2007 May 14; in original form 2007 January 31

## SUMMARY

Rapid progress in imaging deep Earth structures using seismic tomography and in studies of physical and chemical properties of mantle rocks facilitates research in assimilation of data related to mantle dynamics. In this paper, we present a new numerical approach for data assimilation, which allows for incorporating observations (at present) and unknown initial conditions (in the past) for mantle temperature and flow into a 3-D dynamic model in order to determine the initial conditions. The dynamic model is described by the backward heat, motion and continuity equations. The use of the quasi-reversibility (QRV) method implies the introduction into the backward heat equation of the additional term involving the product of a small regularization parameter and a higher order temperature derivative. The data assimilation in this case is based on a search of the best fit between the forecast model state and the observations by minimizing the regularization parameter. We apply the QRV data assimilation method to restore the evolution of (i) mantle plumes (a synthetic case study) and (ii) the lithospheric slab imaged by teleseismic body-wave tomography in the southeastern Carpathians. For both models the present temperature and mantle flow are assimilated to the geological past, and the prominent features of mantle structures are recovered. We then model the evolution of the mantle structures forward in time starting from the restored state to the present state and estimate the accuracy of the model predictions. The results of the QRV data assimilation are compared to that obtained by the variational (VAR) and backward advection data assimilation. Although the accuracy of the QRV data assimilation is lower than that of the VAR data assimilation, the QRV method does not require any additional smoothing of the input data or filtering of temperature noise as the VAR method does. Based on the results and the comparison of the methods, we consider the QRV method to be a highly promising approach to assimilation of data related to mantle dynamics.

**Key words:** finite-difference method, finite-element method, inverse problem, lithospheric slab, mantle convection, plume.

## 1 INTRODUCTION

Many geophysical (geodynamic) problems can be described by mathematical models, that is, by a set of partial differential equations and boundary and/or initial conditions defined in a specific domain. A mathematical model links the causal characteristics of a geodynamic process with its effects. The causal characteristics of the process include, for example, parameters of the initial and boundary conditions, coefficients of the differential equations, and geometrical parameters of a model domain. The aim of the direct

mathematical problem is to determine the relationship between the causes and effects of the geophysical process and hence to find a solution to the mathematical problem for a given set of parameters and coefficients. 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 geophysical process exists). Inverse problems can be subdivided into time-reverse problems (e.g. to restore the development of a geodynamic process), coefficient problems (e.g. to determine the coefficients of the model equations and/or boundary conditions),

<footnote n="1">¹Geophysikalishes Institut, Universität Karlsruhe, Hertzstr. 16, Karlsruhe 76187, Germany. E-mail: alik.ismail-zadeh@gpi.uni-karlsruhe.de</footnote>

<footnote n="2">²Institut de Physique du Globe de Paris, 4 Place Jussieu, Paris 75252, France</footnote>

<footnote n="3">³International Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences, Profsoyuznaya str. 84/32, Moscow 117997, Russia</footnote>

<footnote n="4">⁴Institute of Mathematics and Mechanics, Ural Branch, Russian Academy of Sciences, S. Kovalevskoy ul. 16, Yekaterinburg 620219, Russia</footnote>

<footnote n="5">⁵Department of Earth and Space Sciences & Institute of Geophysics and Planetary Physics, University of California, 3806 Geology Building, 595 Charles Young Drive East, Los Angeles, CA 90095-1567, USA</footnote>
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.

The mantle is heated from the core and from inside due to decay of radioactive elements. Since mantle convection is described by heat advection and diffusion, one can ask: is it possible to tell, from the present temperature estimations of the Earth, something about the Earth's temperature in the geological past? Even though heat diffusion is irreversible in the physical sense, it is possible to predict accurately the heat transfer backwards in time using data assimilation techniques without contradicting the basic thermodynamic laws (e.g. Ismail-Zadeh et al. 2004a).

In this paper, we consider inverse (time-reverse) problems of thermal convection in the Earth's mantle with the aim of restoring mantle dynamics in the geological past. In other words, the present observations (mantle temperature and velocity) can be assimilated into the past to constrain the initial conditions for the mantle temperature and velocity. Data assimilation can be defined as the incorporation of present (observations) and past data (initial conditions) in an explicit dynamic model to provide time continuity and coupling among the physical fields (e.g. velocity and temperature). The basic principle of data assimilation is to consider the initial condition as a control variable and to optimize the initial condition in order to minimize the discrepancy between the observations and the solution of the model.

Inverse problems are often ill-posed. Jacques Hadamard introduced the idea of well- (and ill-)posed problems in the theory of partial differential equations (Hadamard 1902). A mathematical model for a geophysical problem has to be well posed in the sense that it has to have the properties of existence, uniqueness and stability of a solution to the problem. Problems for which at least one of these properties does not hold are called ill posed. The requirement of stability is the most important one. 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.

The inverse problem of thermal convection in the mantle is an ill-posed problem, since the backward heat problem, describing both heat advection and conduction through the mantle backwards in time, possesses the properties of ill-posedness (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 mantle temperatures in the past (e.g. see Appendix A, Ismail-Zadeh et al. 2004a).

Despite the fact that many inverse problems are ill-posed, there are methods for solving the problems. Andrei Tikhonov introduced the idea of conditionally well-posed problems and the regularization method (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. The Tikhonov regularization is essentially a trade-off between fitting the observations and reducing a norm of the solution to the mathematical model of a geophysical problem.

Three principal techniques are currently employed to assimilate data related to a thermoconvective mantle flow: (i) backward advection, (ii) sequential filtering and (iii) variational (adjoint) method.

If heat diffusion is neglected, the present mantle temperature and flow can be assimilated using the backward advection (BAD) into the geological past. A numerical approach to the solution of the inverse problem of the Rayleigh-Taylor instability was proposed by Ismail-Zadeh (1999) and was developed later for a dynamic restoration of diapiric structures to their earlier stages (Ismail-Zadeh et al. 2001a; Kaus & Podladchikov 2001; Korotkii et al. 2002; Ismail-Zadeh et al. 2004b). Steinberger & O'Connell (1998) and Conrad & Gurnis (2003) modelled the mantle flow backwards in time from present-day mantle density heterogeneities inferred from seismic observations. Both direct (forward in time) and inverse (backward in time) problems of the heat (density) advection are well posed. This is because the time-dependent advection equation has the same form of characteristics for the direct and inverse velocity field (the vector velocity reverses its direction, when time is reversed). Therefore, numerical algorithms used to solve the direct problem of the gravitational instability can also be used in studies of the time-reverse problems by replacing positive time steps with negative ones.

In sequential filtering a numerical model is computed forward in time for the interval for which observations have been made, updating the model each time where observations are available. Bunge et al. (1998, 2002) used this approach to compute mantle circulation models. Despite sequential data assimilation well adapted to mantle circulation studies, each individual observation influences the model state at later times. Information propagates from the geological past into the future, although our knowledge of the Earth's mantle at earlier times is much poor than at present.

The variational (VAR) data assimilation method has been pioneered by meteorologists and used very successfully to improve operational weather forecasts (e.g. Kalnay 2003). The use of VAR data assimilation in models of mantle dynamics (to estimate mantle temperature and flow in the geological past) has been put forward by Bunge et al. (2003) and Ismail-Zadeh et al. (2003a,b). The major differences between the two approaches are that Bunge et al. (2003) applied the VAR method to the coupled Stokes, continuity and heat equations (generalized inverse), whereas Ismail-Zadeh et al. (2003a) applied the VAR method to the heat equation only. The VAR approach by Ismail-Zadeh et al. (2003a) is computationally less expensive, because it does not involve the Stokes equation into the iterations between the direct and adjoint problems. Moreover, this approach admits the use of temperature-dependent viscosity.

The VAR data assimilation algorithm was employed for numerical restoration of models of present prominent mantle plumes to their past stages (Ismail-Zadeh et al. 2004a; Hier-Majumder et al. 2005). Effects of thermal diffusion and temperature-dependent viscosity on the evolution of mantle plumes was studied by Ismail-Zadeh et al. (2006) to recover the structure of mantle plumes prominent in the past from that of present plumes weakened by thermal diffusion. Ismail-Zadeh et al. (2006) showed also that smoothness of the input data (present mantle temperature) is a primary factor affecting the accuracy of the VAR data assimilation. An insignificant perturbation of the initial temperature can result in a large increase of the amplitude of the perturbation with time. Smoothing filters are required to reduce the noise (e.g. Samarskii et al. 1997), although the employment of such filters decreases the efficiency of the VAR data assimilation algorithm (computations become time consuming in 3-D cases). Another way to reduce the noise is to employ high-order adjoint techniques (Alekseev & Navon 2001).

In this paper, we present a new approach to assimilation of mantle related data based on a quasi-reversibility (QRV) method (Lattes & Lions 1969). We describe the QRV method in Section 2, present the mathematical statement of the direct and inverse problems of
thermal convection in the mantle in Section 3 and a numerical approach to the QRV data assimilation in Section 4. Numerical results are presented in Section 5 for two models of mantle dynamics: evolution of (i) mantle plumes (a synthetic case study) and (ii) a lithospheric slab imaged by the teleseismic body-wave tomography. In Section 6, we discuss the efficiency of the proposed numerical approach comparing it with that based on the VAR and BAD methods, and finally derive our conclusion.

## 2 QUASI-REVERSIBILITY METHOD

The principal idea of the QRV method is based on the transformation of an ill-posed problem into a well-posed problem (Lattes & Lions 1969). In the case of the backward heat equation, this implies an introduction of an additional term into the equation, which involves the product of a small regularization parameter and higher order temperature derivative. The additional term should be sufficiently small compared to other terms of the heat equation and allow for simple additional boundary conditions. The data assimilation in this case is based on a search of the best fit between the forecast model state and the observations by minimizing the regularization parameter. The QRV method is proven to be well suited for smooth and non-smooth input data (Lattes & Lions 1969; Samarskii & Vabishchevich 2004).

To explain the transformation of the problem, we consider the following boundary-value problem for the 1-D heat conduction problem

$$ \frac{\partial T(t, x)}{\partial t} = \frac{\partial^2 T(t, x)}{\partial x^2}, \quad 0 \le x \le \pi, \quad 0 \le t \le t^*, \qquad (1) $$

$$ T(t, x = 0) = T(t, x = \pi) = 0, \quad 0 \le t \le t^*, \qquad (2) $$

$$ T(t = 0, x) = \frac{1}{4n + 1} \sin[(4n + 1)x], \quad 0 \le x \le \pi. \qquad (3) $$

The analytical solution to (1)–(3) can be obtained in the following form

$$ T(t, x) = \frac{1}{4n + 1} \exp[-(4n + 1)^2 t] \sin[(4n + 1)x]. \qquad (4) $$

Fig. 1 presents the solution (red solid curves) for time interval $0 \le t \le t^* = 0.14$ and $n = 1$.

It is known that the backward heat conduction problem is ill posed (e.g. Kirsh 1996). To transform the problem into a well-posed problem, we introduce a term in eq. (1) involving the product of a small parameter $\beta > 0$ and higher order temperature derivative:

$$ \frac{\partial T_{\beta}(t, x)}{\partial t} = \frac{\partial^2 T_{\beta}(t, x)}{\partial x^2} - \beta \frac{\partial^4}{\partial x^4} \left[ \frac{\partial T_{\beta}(t, x)}{\partial t} \right], \quad 0 \le x \le \pi, \quad 0 \le t \le t^* \quad (5) $$

$$ T_{\beta}(t, x = 0) = T_{\beta}(t, x = \pi) = 0, \quad 0 \le t \le t^*, \qquad (6) $$

$$ \frac{\partial^2 T_{\beta}(t, x = 0)}{\partial x^2} = \frac{\partial^2 T_{\beta}(t, x = \pi)}{\partial x^2} = 0, \quad 0 \le t \le t^*, \qquad (7) $$

$$ T_{\beta}(t = t^*, x) = \frac{1}{4n + 1} \times \exp[-(4n + 1)^2 t^*] \sin[(4n + 1)x], \quad 0 \le x \le \pi. \qquad (8) $$

Here the initial condition is assumed to be the solution (4) to the heat conduction problem (1)–(3) at $t = t^*$. The subscript $\beta$ at $T_{\beta}$ is used to emphasize the dependence of the solution to problem (5)–(8) on the

<figure><img src="image_2.png" /><figcaption>Figure 1. Comparison of the exact solutions to the heat conduction problem (red curves; a and b) and to the regularized backward heat conduction problem (a: $\beta = 10^{-3}$ and b: $\beta = 10^{-7}$; blue dashed curves). The temperature residual between two solutions is presented in panel (c) at various values of the regularization parameter $\beta$ and $x \in [0, \pi/5]$.</figcaption></figure>

regularization parameter. The analytical solution to the regularized backward heat conduction problem (5)–(8) is represented as:

$$ T_{\beta}(t, x) = A_n \exp\left[\frac{-(4n + 1)^2 t}{1 + \beta(4n + 1)^4}\right] \sin[(4n + 1)x], \\ A_n = \frac{1}{4n + 1} \exp[-(4n + 1)^2 t^*] \exp^{-1}\left[\frac{-(4n + 1)^2 t^*}{1 + \beta(4n + 1)^4}\right], \qquad (9) $$
and the solution approaches the initial condition for the problem (1)–(3) at $t = 0$ and $\beta \to 0$. Figs 1(a) and (b) illustrates the solution to the regularized problem at two values of $\beta$ (blue dashed curves) and $n = 1$. The temperature residual (Fig. 1c) indicates that the solution (9) approaches the solution (4) with $\beta \to 0$.

Samarskii & Vabischevich (2004) estimated the stability of the solution to problem (5)–(7) with respect to the initial condition expressed in the form $T_{\beta}(t = t^*, x) = T_{\beta}^*$:

$$ \|T_{\beta}(t, x)\| + \beta \|\partial T_{\beta}(t, x)/\partial x\| \le C(\|T_{\beta}^*\| + \beta \|\partial T_{\beta}^*/\partial x\|) \exp[(t^* - t)\beta^{-1/2}], \quad (10) $$

where $C$ is a constant, and showed that the natural logarithm of errors will increase in a direct proportion to time and inversely to the root square of the regularization parameter.

The transformation to the regularized backward heat problem is not only a mathematical approach to solving ill-posed backward heat problems, but has some physical meaning: it can be explained on the basis of the concept of relaxing heat flux for heat conduction (e.g. Vernotte 1958). The classical Fourier heat conduction theory provides the infinite velocity of heat propagation in a region. The instantaneous heat propagation is unrealistic, because the heat is a result of the vibration of atoms and the vibration propagates in a finite speed (Morse & Feshbach 1953). To accommodate the finite velocity of heat propagation, a modified heat flux model was proposed by Vernotte (1958) and Cattaneo (1958).

The modified Fourier constitutive equation is expressed as $\vec{Q} = -k\nabla T - \tau \partial^2 \vec{Q} / \partial t^2$, where $\vec{Q}$ is the heat flux, and $k$ is the coefficient of thermal conductivity. The thermal relaxation time $\tau = k/(\rho c_p v^2)$ is usually recognized to be a small parameter (Yu et al. 2004), where $\rho$ is the density, $c_p$ is the specific heat and $v$ is the heat propagation velocity. The situation for $\tau \to 0$ leads to instantaneous diffusion at infinite propagation speed, which coincides with the classical thermal diffusion theory. The heat conduction equation $\partial T / \partial t = \nabla^2 T + \tau \partial^2 T / \partial t^2$ based on non-Fourier heat flux can be considered as a regularized heat equation. If the Fourier law is modified further by an addition of the second derivative of heat flux, for example, $\vec{Q} = -k\nabla T + \beta \frac{\partial^2 \vec{Q}}{\partial t^2}$, where small $\beta$ is the relaxation parameter of heat flux (Bubnov 1976, 1981), the heat conduction equation can be transformed into a higher order regularized heat equation similar to eq. (5).

### 3 MATHEMATICAL STATEMENT OF THE PROBLEM

We consider a model of thermoconvective mantle flow in the 3-D domain $\Omega = [0, x_1 = l_1] \times [0, x_2 = l_2] \times [0, x_3 = l_3 = h]$, where $\mathbf{x} = (x_1, x_2, x_3)$ are the Cartesian coordinates. We assume that the mantle behaves as a Newtonian incompressible fluid with a temperature-dependent viscosity and infinite Prandtl number. Rising mantle plumes (subducting lithosphere) are modelled as a hot (cold) viscous fluid ascending (descending) into the relatively cold (hot) ambient viscous fluid heated from below. The mantle flow is described by heat, motion and continuity equations (Chandrasekhar 1961). To simplify the governing equations, we make the Boussinesq approximation (Boussinesq 1903) keeping the density constant everywhere except for the buoyancy term in the equation of motion. We note that a temperature-dependent density, internal heating, phase transformations in the mantle and other physical complications can be also considered in the modelling of thermoconvective mantle circulation (e.g. Schubert et al. 2001; Ismail-Zadeh et al. 2003b). In

the Boussinesq approximation the motion of a viscous fluid is described by the following dimensionless formulation: the boundary value problem for the flow velocity (it includes the Stokes equation, the incompressibility equation subject to appropriate boundary conditions)

$$ \nabla P = \operatorname{div}(\eta(T)\mathbf{E}) + RaT\mathbf{e}, \quad \mathbf{x} \in \Omega, \quad (11) $$

$$ \operatorname{div}\mathbf{u} = 0, \quad \mathbf{x} \in \Omega, \quad (12) $$

$$ \mathbf{u} \cdot \mathbf{n} = 0, \quad \mathbf{x} \in \partial\Omega, \quad (13) $$

$$ \partial \mathbf{u}_t / \partial \mathbf{n} = 0, \quad \mathbf{x} \in \partial \Omega, \quad (14) $$

and the initial-boundary value problem for temperature (it includes the heat equation subject to appropriate boundary and initial conditions)

$$ \partial T / \partial t + \mathbf{u} \cdot \nabla T = \nabla^2 T + f, \quad t \in [0, \vartheta], \quad \mathbf{x} \in \Omega, \quad (15) $$

$$ \sigma_1 T + \sigma_2 \partial T / \partial \mathbf{n} = T_*, \quad t \in [0, \vartheta], \quad \mathbf{x} \in \partial \Omega, \quad (16) $$

$$ T(0, \mathbf{x}) = T_0(\mathbf{x}), \quad \mathbf{x} \in \Omega. \quad (17) $$

Here $T$, $t$, $\mathbf{u}$, $P$ and $\eta$ are dimensionless temperature, time, velocity, pressure and viscosity, respectively; $\mathbf{E} = e_{ij}(\mathbf{u}) = \{\partial u_i / \partial x_j + \partial u_j / \partial x_i\}$ is the strain rate tensor; $u_i$ are the velocity components; $\mathbf{e} = (0, 0, 1)$ is the unit vector; $\nabla$ is the gradient operator; $\text{div}$ is the divergence operator; $f$ is the heat source; $\mathbf{n}$ is the outward unit normal vector at a point on the model boundary; $\mathbf{u}_t$ is the projection of the velocity vector onto the tangent plane at the same point on the model boundary; $[t = 0, t = \vartheta]$ is the model time interval; $\sigma_1$ and $\sigma_2$ are some piecewise smooth functions or constants such that $\sigma_1^2 + \sigma_2^2 \neq 0$.

We consider the impermeability condition with perfect slip on $\partial \Omega$. Choosing $\sigma_1$, $\sigma_2$ and $T_*$ in a proper way we can specify temperature or heat flux at the model boundaries. By $\Gamma_u = \{\mathbf{x} : (\mathbf{x} \in \partial \Omega) \cap (x_3 = l_3)\}$, $\Gamma_l = \{\mathbf{x} : (\mathbf{x} \in \partial \Omega) \cap (x_3 = 0)\}$, and $\Gamma_v = \bigcup_{i=1,2} \{\mathbf{x} : (\mathbf{x} \in \Omega) \cap (x_i = 0)\} \cup \{\mathbf{x} : (\mathbf{x} \in \Omega) \cap (x_i = l_i)\}$, we denote the parts of the model boundary that $\Gamma_u \cup \Gamma_l \cup \Gamma_v = \partial \Omega$. We assume the constant temperature at the horizontal boundaries and zero heat flux at vertical boundaries of the model domain: $\sigma_1(t, \mathbf{x}) = 1$, $\sigma_2(t, \mathbf{x}) = 0$ and $T_*(t, \mathbf{x}) = 0$ at $(t, \mathbf{x}) \in [0, \vartheta] \times \Gamma_u; \sigma_1(t, \mathbf{x}) = 1$, $\sigma_2(t, \mathbf{x}) = 0$ and $T_*(t, \mathbf{x}) = 1$ at $(t, \mathbf{x}) \in [0, \vartheta] \times \Gamma_l$; and $\sigma_1(t, \mathbf{x}) = 0$, $\sigma_2(t, \mathbf{x}) = 1$ and $T_*(t, \mathbf{x}) = 0$ at $(t, \mathbf{x}) \in [0, \vartheta] \times \Gamma_v$.

The Rayleigh number is defined as $Ra = \alpha g \rho_{\text{ref}} \Delta T h^3 \eta_{\text{ref}}^{-1} \kappa^{-1}$, where $\alpha$ is the thermal expansivity, $g$ is the acceleration due to gravity, $\rho_{\text{ref}}$ and $\eta_{\text{ref}}$ are the reference typical density and viscosity, respectively; $\Delta T$ is the temperature contrast between the lower and upper boundaries of the model domain; and $\kappa$ is the thermal diffusivity. Length, temperature and time are normalized by $h$, $\Delta T$ and $h^2 \kappa^{-1}$, respectively. The physical parameters of the fluid (temperature, velocity, pressure, viscosity and density) are assumed to depend on time and on space coordinates. The viscosity (specified later) and density depend on temperature as well.

The direct problem of thermoconvective flow is formulated as follows: find the velocity $\mathbf{u} = \mathbf{u}(t, \mathbf{x})$, the pressure $P = P(t, \mathbf{x})$ and the temperature $T = T(t, \mathbf{x})$ satisfying boundary value problem (11)–(14) and initial-boundary value problem (15)–(17). We can formulate the inverse problem in this case as follows: find the velocity, pressure and temperature satisfying boundary value
problem (11)–(14) and the final-boundary value problem which includes eqs (15) and (16) and the final condition:

$$T(\vartheta, \mathbf{x}) = T_{\vartheta}(\mathbf{x}), \quad \mathbf{x} \in \Omega, \qquad (18)$$

where $T_{\vartheta}$ is the temperature at time $t = \vartheta$.

To solve the inverse problem by the QRV method we consider the following regularized backward heat problem to define temperature in the past from the known temperature $T_{\vartheta}(\mathbf{x})$ at present time $t = \vartheta$:

$$\begin{aligned} \partial T_{\beta}/\partial t - \mathbf{u}_{\beta} \cdot \nabla T_{\beta} \\ = \nabla^2 T_{\beta} + f - \beta \Lambda(\partial T_{\beta}/\partial t), \quad t \in [0, \vartheta], \quad \mathbf{x} \in \Omega, \end{aligned} \qquad (19)$$

$$\sigma_1 T_{\beta} + \sigma_2 \partial T_{\beta} / \partial \mathbf{n} = T_{*}, \quad t \in (0, \vartheta), \quad \mathbf{x} \in \partial \Omega, \qquad (20)$$

$$\sigma_1 \partial^2 T_{\beta} / \partial \mathbf{n}^2 + \sigma_2 \partial^3 T_{\beta} / \partial \mathbf{n}^3 = 0, \quad t \in (0, \vartheta), \quad \mathbf{x} \in \partial \Omega, \qquad (21)$$

$$T_{\beta}(\vartheta, \mathbf{x}) = T_{\vartheta}(\mathbf{x}), \quad \mathbf{x} \in \Omega, \qquad (22)$$

where $\Lambda(T) = \partial^4 T / \partial x_1^4 + \partial^4 T / \partial x_2^4 + \partial^4 T / \partial x_3^4$, and the boundary value problem to determine the fluid flow:

$$\nabla P_{\beta} = -\operatorname{div}[\eta(T_{\beta})\mathbf{E}(\mathbf{u}_{\beta})] + R a T_{\beta} \mathbf{e}, \quad \mathbf{x} \in \Omega, \qquad (23)$$

$$\operatorname{div} \mathbf{u}_{\beta} = 0, \quad \mathbf{x} \in \Omega, \qquad (24)$$

$$\mathbf{u}_{\beta} \cdot \mathbf{n} = 0, \quad \mathbf{x} \in \partial \Omega, \qquad (25)$$

$$\partial(\mathbf{u}_{\beta})_{\tau} / \partial \mathbf{n} = 0, \quad \mathbf{x} \in \partial \Omega, \qquad (26)$$

where the sign of the velocity field is changed ($\mathbf{u}_{\beta}$ by $-\mathbf{u}_{\beta}$) in eqs (19) and (23) to simplify the application of the total variation diminishing method for solving (19)–(22) (see Section A2 of Appendix A). Hereinafter, we refer to temperature $T_{\vartheta}$ as the input temperature for the problem (19)–(26). The core of the transformation of the heat equation is the addition of a high order differential expression $\Lambda(\partial T_{\beta}/\partial t)$ multiplied by a small parameter $\beta > 0$. Note that eq. (21) is added to the boundary conditions to properly define the regularized backward heat problem. Parameters $\sigma_1$ and $\sigma_2$ in (20) and (21) are the same as in (16). The solution to the regularized backward heat problem is stable for $\beta > 0$, and the approximate solution to (19)–(26) converges to the solution of (11)–(16), and (18) in some spaces, where the conditions of well posedness are met (Samarskii & Vabischevich 2004). Thus, the inverse problem of thermoconvective mantle flow is reduced to determination of the velocity $\mathbf{u}_{\beta} = \mathbf{u}_{\beta}(t, \mathbf{x})$, the pressure $P_{\beta} = P_{\beta}(t, \mathbf{x})$ and the temperature $T_{\beta} = T_{\beta}(t, \mathbf{x})$ satisfying (19)–(26).

# 4 NUMERICAL APPROACH

## 4.1 Numerical methods

To solve the heat problem (15)–(17) and the regularized backward heat problem (19)–(22), finite differences are used to derive discrete equations. We employ (i) the characteristic-based semi-Lagrangian (CBSL) method (Courant 1952; Staniforth & Coté 1991) to calculate the derivatives of the convective term in the heat eq. (15); (ii) the total variation diminishing (TVD) method (Harten 1983) to calculate the derivatives of the convective term in the regularized backward heat eq. (19); (iii) central differences to approximate the derivatives of the diffusion and regularizing terms in (15) and (19), respectively and (iv) the two-layered additively averaged scheme to represent the 3-D spatial discrete operators associated with the diffusion and

regularizing terms as 1-D discrete operators, and the componentwise splitting method to solve the set of the discrete equations (Samarskii & Vabischevich 1995). We describe the numerical methods in detail in Section A1 (for the heat problem) and Section A2 (for the regularized backward heat problem).

The Eulerian finite-element method is employed to solve the Stokes problems (11)–(14) and (23)–(26). The numerical approach is based on the representation of the flow velocity by a two-component vector potential (Ismail-Zadeh et al. 2001b) eliminating the incompressibility equation from the relevant boundary value problems. This potential is approximated by tri-cubic splines, which allows one to efficiently interpolate the velocity field (Ismail-Zadeh et al. 1998). Such a procedure results in a set of linear algebraic equations with a symmetric positive-definite banded matrix. We solve the set of discrete equations by the conjugate gradient method (Fletcher & Reeves 1964) using parallel processors. A description of this approach can be found in Ismail-Zadeh et al. (2004b).

## 4.2 Optimization problem

We seek a maximum of the following functional with respect to the regularization parameter $\beta$:

$$\delta - \|T(t = \vartheta, \cdot; T_{\beta_k}(t = 0, \cdot)) - \varphi(\cdot)\| \to \max_k, \qquad (27)$$

$$\beta_k = \beta_0 q^{k-1}, \quad k = 1, 2, \dots, \aleph, \qquad (28)$$

where sign $\|\cdot\|$ denotes the norm in the space $L_2(\Omega)$, the Hilbert space with the norm defined as $\|v\| = [\int_{\Omega} y^2(\mathbf{x}) d\mathbf{x}]^{1/2}$. Since in what follows the dependence of solutions on initial temperature data is important, we introduce these data explicitly into the mathematical representation of temperature. Here $T_k = T_{\beta_k}(t = 0, \cdot)$ is the solution to the regularized backward heat problem (19)–(22) at $t = 0$; $T(t = \vartheta, \cdot; T_k)$ is the solution to the heat problem (15)–(17) at the initial condition $T(t = 0, \cdot) = T_k$ at time $t = \vartheta$; $\varphi$ is the known temperature at $t = \vartheta$ (the input data on the present temperature); small parameters $\beta_0 > 0$ and $0 < q < 1$ are defined below and $\delta > 0$ is a given accuracy. When $q$ tends to unity, the computational cost becomes large; and when $q$ tends to zero, the optimal solution can be missed.

The prescribed accuracy $\delta$ is composed from the accuracy of the initial data and the accuracy of computations. When the input noise decreases and the accuracy of computations increases, the regularization parameter is expected to decrease. However, estimates of the initial data errors are usually inaccurate. Estimates of the computation accuracy are not always known, and when they are available, the estimates are coarse. In practical computations, it is more convenient to minimize the following functional with respect to (28)

$$\| T_{\beta_{k+1}}(t = 0, \cdot) - T_{\beta_k}(t = 0, \cdot) \| \to \min_k, \qquad (29)$$

where misfit between temperatures obtained at two adjacent iterations must be compared. To implement the minimization of temperature residual (27), the inverse problem (19)–(26) must be solved on the entire time interval as well as the direct problem (11)–(17) on the same time interval. This at least doubles the amount of computations. The minimization of functional (29) has a lower computational cost, but it does not rely on a priori information.

## 4.3 Numerical algorithm for QRV data assimilation

In this section we describe the numerical algorithm for solving the inverse problem of thermoconvective mantle flow using the QRV
method. For simplicity consider a uniform temporal partition $t_n = \vartheta - \tau n$, and $n$ successively takes integer values from 0 to some natural number $N = \vartheta/\tau$. We prescribe some values to parameters $\beta_0, q$ and $\Re$ (e.g. $\beta_0 = 10^{-3}$, $q = 0.1$ and $\Re = 10$). According to (28) we define a sequence of the values of the regularization parameter $\{\beta_k\}$. For each value $\beta = \beta_k$ model temperature and velocity are determined in the following way.

**Step 1:** Given the temperature $T_\beta = T_\beta(t, \cdot)$ at $t = t_n$, the velocity $\mathbf{u}_\beta = \mathbf{u}_\beta(t_n, \cdot)$ is found by solving problem (23)–(26). This velocity is assumed to be constant on the time interval $[t_{n+1}, t_n]$.

**Step 2:** Given the velocity $\mathbf{u}_\beta = \mathbf{u}_\beta(t_n, \cdot)$, the new temperature $T_\beta = T_\beta(t, \cdot)$ at $t = t_{n+1}$ is found on the time interval $[t_{n+1}, t_n]$ subject to the final condition $T_\beta = T_\beta(t_n, \cdot)$ by solving problem (19)–(22) according to the numerical method described in Section A2.

**Step 3:** Upon the completion of steps 1 and 2 for all $n = 0, 1, \dots, N$, the temperature $T_\beta = T_\beta(t_n, \cdot)$ and the velocity $\mathbf{u}_\beta = \mathbf{u}_\beta(t_n, \cdot)$ are obtained at each $t = t_n$. Based on the computed solution we can find the temperature and flow velocity at each point of time interval $[0, \vartheta]$ using interpolation.

**Step 4a:** The direct problem (15)–(17) is solved assuming that the initial temperature is given as $T_\beta = T_\beta(t = 0, \cdot)$, and the temperature residual (27) is found. If the residual does not exceed the pre-defined accuracy, the calculations are terminated, and the results obtained at step 3 are considered as the final ones. Otherwise, parameters $\beta_0, q$ and $\Re$ entering eq. (28) are modified, and the calculations are continued from step 1 for new set $\{\beta_k\}$.

**Step 4b:** The functional (29) is calculated. If the residual between the solutions obtained for two adjacent regularization parameters satisfies a pre-defined criterion (the criterion should be defined by a user, because no *a priori* data are used at this step), the calculation is terminated, and the results obtained at step 3 are considered as the final ones. Otherwise, parameters $\beta_0, q$ and $\Re$ entering eq. (28) are modified, and the calculations are continued from step 1 for new set $\{\beta_k\}$.

In a particular implementation, either step 4a or step 4b is used to terminate the computation. This algorithm allows (i) organizing a certain number of independent computational modules for various values of the regularized parameter $\beta_k$ that find the solution to the regularized problem using steps 1–3 and (ii) determining a *posteriori* an acceptable result according to step 4a or 4b.

# 5 METHOD'S APPLICATION

The reconstruction of mantle plumes and lithospheric slabs to earlier stages of their evolution is a major challenge in geodynamics. High-resolution global and regional seismic tomographic studies open possibilities for detailed observations of present-day mantle structures (e.g. Montelli et al. 2004; Martin et al. 2006) and for derivations of mantle temperature from seismic velocities or velocity anomalies (e.g. Goes et al. 2000; Ismail-Zadeh et al. 2005a). An accurate reconstruction would allow the test of geodynamic models by simulating the evolution of plumes or slabs starting from the restored state and comparing the derived forward state to observations.

We present the numerical results for the restoration of mantle plume evolution (synthetic case study) in Section 5.1 and for the assimilation of a present temperature model associated with a descending lithospheric slab imaged by high-resolution regional seismic tomography in Section 5.2.

## 5.1 Mantle plumes

Thermal plumes in the Earth's mantle plausibly originate near either the core-mantle boundary or the upper mantle-lower mantle transition due to instabilities in the hot thermal boundary layers that could exist at these locations. Although some mantle plumes appear to last for more than 150 Myr, they are nonetheless transient features (Schubert et al. 2001).

Mantle plumes evolve in three distinguishing stages: (i) imma-
ture, that is, an origin and initial rise of the plumes; (ii) mature, that
is, plume-lithosphere interaction, gravity spreading of plume head
and development of overhangs beneath the bottom of the lithosphere,
and partial melting of the plume material (e.g. Ribe & Christensen
1994; Moore et al. 1998) and (iii) overmature, that is, slowing-down
of the plume rise and fading of the mantle plumes due to thermal
diffusion (Davaille & Vatteville 2005; Ismail-Zadeh et al. 2006).
The ascent and evolution of mantle plumes depend on the proper-
ties of the source region (i.e., the thermal boundary layer) and the
viscosity and thermal diffusivity of the ambient mantle. The proper-
ties of the source region determine temperature and viscosity of the
mantle plumes. Structure, flow rate and heat flux of the plumes are
controlled by the properties of the mantle through which the plumes
rise. While properties of the lower mantle (e.g. viscosity and ther-
mal conductivity) are relatively constant during the approximately
150 Myr lifetime of most plumes, source region properties can vary
substantially with time as the thermal basal boundary layer feeding
the plume is depleted of hot material. Complete local depletion of
this boundary layer cuts the plume off from its source.

We start our simulations by computing a forward model of the
evolution of the thermal plumes and then we restore the evolved
plumes to their earlier stages. To compare the numerical results
obtained by the QRV method with that obtained by the VAR and
BAD methods, we develop the same forward model for mantle plume
evolution as presented by Ismail-Zadeh et al. (2006).

The evolution of mantle plumes originating at the core-mantle
boundary is modelled through numerical experiments of 3-D ther-
mal convection in a bottom heated box. The mantle behaves as
a Newtonian fluid on geological timescales, and a dimensionless
temperature-dependent viscosity law (Busse et al. 1993)

$$ \eta(T) = \exp[M(T + G)^{-1} - M(0.5 + G)^{-1}] \quad (30) $$

is used in the modelling, where $M = [225/\ln(r)] - 0.25 \ln(r)$, $G = [15/\ln(r)] - 0.5$ and $r$ is the viscosity ratio between the upper and lower boundaries of the model domain. We model the evolution of mantle plumes for three viscosity profiles: $r = 20, 200$ and $1000$. The temperature-dependent viscosity profiles have their minimum at the core-mantle boundary. A more realistic viscosity profile (e.g. Forte & Mitrovica 2001) will influence the evolution of mantle plumes, though it will not influence the restoration of the plumes. A mathematical model of mantle plume dynamics is described by a set of eqs (11)–(17), where many complications are omitted, for example a viscosity increase from the upper to the lower mantle, adiabatic heating/cooling, internal heating (we use $f = 0$ in the modelling), phase transformations in the mantle and others. Inclusion of these complications would refine the model.

With $\alpha = 3 \times 10^{-5} \text{ K}^{-1}$, $\rho_{\text{ref}} = 4000 \text{ kg m}^{-3}$, $\Delta T = 3000 \text{ K}$, $h = 2800 \text{ km}$, $\eta_{\text{ref}} = 8 \times 10^{22} \text{ Pa s}$ and $\kappa = 10^{-6} \text{ m}^2 \text{ s}^{-1}$, the initial Rayleigh number is $Ra = 9.5 \times 10^5$. While plumes evolve in the convecting heterogeneous mantle, at the initial time we assume that the plumes develop in a laterally homogeneous temperature field and hence consider that the mantle temperature in the model increases linearly with depth. The model domain is divided into $38 \times 38 \times$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Model of mantle plume evolution forward in time at successive times: (a–d) from 265 Myr ago to the present state of the plumes ($r = 20$). Assimilation of the mantle temperature and flow from the present state back to the geological past using the QRV (d–g; $\beta = 10^{-7}$) and BAD (d, e1–g1) methods. Verification of the QRV assimilation accuracy: Forward model of the plume evolution starting from the initial (restored) state of the plumes (g) to their present state (i). Temperature residuals between the initial temperature for the forward model and the temperature assimilated to the same age using the QRV and BAD methods are presented in panels (h) and (h1), respectively.</figcaption></figure>

30 rectangular finite elements to approximate the vector velocity potential by tricubic splines, and a uniform grid 112 × 112 × 88 is employed for approximation of temperature, velocity and viscosity. We employed 36 processors to solve the problem. Each time step for the restoration took less than 1 min, and about 200 time steps were performed.

Figs 2(a)–(d) illustrate the evolution of mantle plumes in the forward model. Mantle plumes are generated by random temperature perturbations at the top of the thermal source layer associated with the core–mantle boundary. The mantle material in the basal source layer flows horizontally toward the plumes. The reduced viscosity in this basal layer promotes the flow of the material to the plumes. Vertical upwelling of hot mantle material is concentrated in low viscosity conduits near the centrelines of the emerging plumes (Figs 2b and c). The plumes move upward through the model domain, gradually forming structures with well-developed heads and tails. The plumes diminish in size with time, and the plume tails are diffused before the plume heads (Fig. 2d).

The state of the plumes at the ‘present’ time (Fig. 2d) obtained by solving the direct problem was used as the input temperature for the inverse problem (an assimilation of the ‘present’ temperature to the past). Note that this initial state (input temperature) is given with an error introduced by the numerical algorithm used to solve the direct problem. Fig. 2 illustrates the states of the plumes restored by the QRV method (panels e–g) and the residual (panel h)

$$ \delta T(x_1, x_2) = \left\{ \int_{0}^{h} [T(x_1, x_2, x_3) - \tilde{T}(x_1, x_2, x_3)]^2 dx_3 \right\}^{1/2} \quad (31) $$

between the initial temperature for the forward model (Fig. 2a) and the temperature $\tilde{T}(x)$ assimilated to the same age (Fig. 2g). To check the stability of the algorithm, we run a forward model of the restored plumes using the solution to the inverse problem at the time of 265 Myr ago (Fig. 2g) as the initial state for the forward model. The result of this run is shown in Fig. 2(i).

To compare the accuracy of the data assimilation methods we develop a restoration model from the ‘present’ time (Fig. 2d) to the time of 265 Myr ago using the BAD method. Fig. 2 shows the model results (panels e1–g1) together with the temperature residual (panel h1) between the initial temperature (panel a) and the temperature assimilated to the same age (panel g1). According to our experience (Ismail-Zadeh et al. 2004, 2006) the VAR method cannot be used to assimilate data within the time interval of more than about 100 Myr (for $Ra \approx 10^6$) without proper filtering of the increasing noise due to non-smoothness of the input data and solution. Hence we do not apply the VAR method to assimilate the synthetic data for 265 Myr.

Fig. 3(a) presents the residual $J_1(\beta) = \|T_0(\cdot) - T_\beta(t = t_0, \cdot; T_\vartheta)\|$ between the initial temperature $T_0$ at $t_0 = 265$ Myr ago and the restored temperature (to the same time) obtained by solving the inverse problem with the input temperature $T_\vartheta$. The optimal accuracy is attained at $\beta^* = \arg\min\{J_1(\beta) : \beta = \beta_k, k = 1, 2, \dots, 10\} \approx 10^{-7}$ in the case of $r = 20$, and at $\beta^* \approx 10^{-6}$ and $\beta^* \approx 10^{-5.5}$ in the cases of $r = 200$ and 1000, respectively. Fig. 3(b) illustrates the residual $J_2(\beta) = \|T_\beta(t_0, \cdot; T_\vartheta) - T_{\tilde{\beta}}(t_0, \cdot; T_\vartheta)\|$ between the reconstructed temperature at $t_0 = 265$ Myr ago obtained for various values of $\beta$ in the range $10^{-9} \le \beta \le 10^{-3}$ and $\tilde{\beta} = \beta/2$. These results show the choice of the optimal value of the regularization parameter using step 4b of the numerical algorithm for the QRV data assimilation (Section 4.3). In the case of $r = 20$ the parameter $\beta^* = \arg\min\{J_2(\beta) : \beta = \beta_k, k = 1, 2, \dots, 12\} \approx 10^{-8}$ provides the optimal accuracy for the solution; in the cases of $r = 200$ and 1000 the optimal accuracy is achieved at $\beta^* \approx 10^{-7}$ and $\beta^* \approx 10^{-6.5}$, respectively. Comparison of the temperature residuals for three values of the viscosity ratio $r$ indicates that the residuals become larger as the viscosity ratio increases (see Fig. 3). The numerical experiments show that the algorithm for solving the inverse problem performs well when the regularization parameter is in the range
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Temperature misfit (a) *J*₁ and (b) *J*₂ as functions of the regularization parameter *β*. The minimum of the temperature misfit is achieved at *β*<sup>*</sup>, an optimal regularization parameter. Solid curves: *r* = 20; dashed curves: *r* = 200; and dash-dotted curves: *r* = 1000.</figcaption></figure>

*10*⁻⁸ ≤ *β* ≤ *10*⁻⁶. For greater values, the solution of the inverse problem retains the stability but is less accurate. For *β* < *10*⁻⁹ the numerical procedure becomes unstable, and the computations must be stopped.

Both laboratory (Davaille & Vatteville 2005) and numerical ex-
periments forward in time (Ismail-Zadeh et al. 2006) show that
thermal plumes start disappearing from bottom up due to a week
feeding of mantle plumes by the hot material from the boundary
layer. To compare how three techniques for data assimilation can
restore the prominent state of the thermal plumes in the past from
their ‘present’ weak state, we develop initially a forward model from
the prominent state of the plumes (Fig. 4a) to their diffusive state in
100 Myr (Fig. 4b) using 50 × 50 × 50 finite rectangular elements
to approximate the vector velocity potential and a finite difference
grid 148 × 148 × 148 for approximation of temperature, velocity
and viscosity. All other parameters of the model are the same as
described in Section 5.1.

We apply the QRV, VAR and BAD methods to restore the plumes
from their weak state and present the results of the restoration and
temperature residuals in Fig. 4. The VAR method (Figs 4d and g)
provides the best performance for the diffused plume restoration.
The BAD method (Figs 4e and h) cannot restore the diffused parts
of the plumes, because temperature is only advected backward in
time. The QRV method (Figs 4c and f) restores the diffused thermal
plumes, meanwhile the restoration results are not so perfect as in
the case of VAR method (compare temperature residuals in Figs 4f
and g). Although the accuracy of the QRV data assimilation is lower
compared to the VAR data assimilation, the QRV method does not
require any additional smoothing of the input data and filtering of
temperature noise as the VAR method does.

5.2 Descending lithosphere

In this section, we present the results of the assimilation of a present
temperature model for the southeastern Carpathians to the geologi-
cal past. To develop the temperature model we use the most recent
high-resolution teleseismic body-wave tomography image (map of
the anomalies of P-wave seismic velocities) of the lithosphere and
asthenosphere in the region (Martin *et al.* 2006). Smearing from
strong crustal velocity anomalies into the upper mantle is success-
fully suppressed by traveltime corrections with an *a priori* 3-D re-
gional crustal velocity model (Martin *et al.* 2005). The image shows
a high velocity body beneath the Vrancea region and the Moesian
platform interpreted by Martin *et al.* (2006) as the subducted litho-
spheric slab. The NE-part of the slab hosts the intermediate-depth
earthquakes and is known as the Vrancea region. The model of
present temperature developed by Ismail-Zadeh *et al.* (2005a) is
based on the previous seismic-tomographic model by Martin *et al.*
(2001), where variations in the crustal thickness were not considered
in traveltime calculations.

The seismic tomographic model of the region consists of eight
vertical layers of different thickness (15–70 km) starting from the
depth of 35 km down to 440 km. Each layer is subdivided hori-
zontally into 16 × 16 km² blocks (Martin et al. 2006). To restrict
numerical errors in our data assimilation we smooth the velocity
anomaly data between the blocks and the layers using a spline in-
terpolation. To convert the P-wave seismic velocity anomalies into
temperature we use the approach described by Goes et al. (2000)
and Ismail-Zadeh et al. (2005a).

Fig. 5(a) (left column) illustrates depth slices of the present tem-
perature model derived from the seismic tomography data. The mod-
elled low mantle temperatures are associated with the high-velocity
body beneath the Vrancea region (VRA) and the East European
platform (EEP). High temperatures are predicted beneath the Tran-
sylvanian Basin (TRB) at about 70–110 km depth. Two other high
temperature regions are found at 110–150 km depth below the Moe-
sian platform (MOP) and deeper than 200 km under the EEP and
the Dobrogea orogen (DOB), which might be correlated with the
regional lithosphere/asthenosphere boundary.

We assimilate the present temperature data into the geological past to restore the prominent thermal features of the lithospheric slab at shallow depths in the region. We use the following parameters in the modelling: *h* = 670 km, the aspect ratio (ratio between horizontal and vertical dimensions of the model) is 1.5, *r* = 1000, Δ*T* = 1700 K, *ρ*<sub>ref</sub> = 3400 kg m<sup>-3</sup>, *η*<sub>ref</sub> = 10<sup>21</sup> Pa s and *Ra* = 5.2 × 10<sup>5</sup>. The present temperature above 440 km depth is derived from the seismic velocity anomalies. We use the adiabatic geotherm for potential temperature 1750 K (Katsura *et al.* 2004) to define the present temperature below 440 km (where seismic tomography data have not been not available). Fig. 5(b) (central column) shows the temperature restored to 22 Myr ago. We assimilate present mantle temperature and flow to the Miocene time, because the post-Miocene descent of the slab is believed to be induced mainly by gravity forces (Csontos *et al.* 1992). The evolution of the southeastern Carpathian region earlier than Miocene, when active subduction of oceanic (or continental) lithosphere was likely to be active, is less clear. The regional horizontal movements are indeed poorly known, and hence the implementation of the poor knowledge into the assimilation model could result in incorrect scenarios of mantle and lithosphere dynamics in the region. Therefore, we have avoided the assimilation of the data beyond the Miocene time.

Early Miocene subduction beneath the Carpathian arc and sub-
sequent gentle continental collision transported cold and dense
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Model of mantle plume diffusion forward in time (a and b; *r* = 20). Assimilation of the mantle temperature and flow to the time of 100 Myr ago and temperature residuals between the present temperature model (b) and the temperature assimilated to the same age, using the QRV (c and f; *β* = 10⁻⁷), VAR (d and g) and BAD (e and h) methods, respectively.</figcaption></figure>

lithospheric material into the hotter mantle. The cold (blue to dark green) region seen at the depths of 40–220 km (see Fig. 5b) can be interpreted as the earlier evolutionary stage of the lithospheric slab. The slab is almost invisible at the shallow depths in the model of the present temperature (see relevant temperature slices in Fig. 5a). Since active subduction of the lithospheric slab in the region has ended in late Miocene times, we argue that the cold slab, descending slowly at these depths due to gravity, has been warmed up, and its thermal shape has faded away due to heat diffusion. Thermal conduction in the shallow Earth (where viscosity is high) plays a significant part in heat transfer compared to thermal convection. The deeper we look into the region, the larger are effects of thermal advection compared to diffusion: the lithosphere has moved upwards to the place where it had been in the Miocene times. At 280 km depth and down to 340 km depth the thermal shape of the slab is clearly visible at the slices of the present temperature model and nearly invisible at the slices of the restored temperature model, because the slab did not reach these depths 22 Myr ago.

The geometry of the restored slab (based on the temperature of about 900 K) shows clearly two portions of the sinking body. One of

them has NW–SE orientation, its location coincides with the boundary between the EEP and Scythian platform (SCP), and this portion of sinking body may be a relic of eastward travelled cold lithosphere. Another portion has a NE–SW orientation and is associated with the present descending slab. The geometry shows that the restored slab is laterally thin compared to the present thick slab at depths below 90 km. This can be explained by the fact that a slab descending gravitationally into the mantle thickens with depth and develops a sheath of lithospheric material with time (e.g. Ismail-Zadeh et al. 2005b).

An interesting geometrical feature of the restored slab is its curvature beneath the SE-Carpathians (Vrancea). The slab has a concave surface, which follows the curvature of the Carpathian arc down to the depths of about 60 km, but at greater depth it changes its shape to that of a convex surface and splits into two parts at a depth of about 200 km. Although such a change in the slab curvature is visible neither in the model of the present temperature nor in the seismic tomography image most likely because of the slab warming and heat diffusion, we suppose that the convex shape of the slab is preserved at the present time. We argue that this abrupt change in the geometry of the descending slab can cause stress localization due to the
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Mantle temperatures (at horizontal slices at depths from 40 to 340 km) for the model of the descending lithospheric slab beneath the southeastern Carpathians. Left panel (a): Present temperature derived from the P-wave velocity anomalies imaged by seismic tomography (Martin *et al.* 2006). Central panel (b): Temperature obtained by assimilation of the present temperature to the Miocene time (22 Myr ago) using the QRV method ($\beta = 10^{-6}$). Right panel (c): Temperature obtained by the backward advection of the present temperature to the same Miocene time using the BAD method. DOB, Dobrogea orogen; EEP, Eastern European platform; MOP, Moesian platform; SCP, Scythian platform; TRB, Transylvanian basin and VRA, Vrancea.</figcaption></figure>

slab bending and subsequent stress release resulting in earthquakes, which occur at the depths of 70–180 km in the region. Moreover the results of the assimilation of the present temperature model to the Miocene time provide a plausible explanation for the change in

the spatial orientation of the slab from NE–SW to NS beneath 200 km observed in the seismic tomography image (Martin *et al.* 2006). The origin of the slab bending is not clear, but it might be related to a complex interaction between two portions of the sinking body
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (Continued.)</figcaption></figure>

and the surrounding mantle. The sinking body displaces the mantle, which, in its turn, forces the slab to deform due to a corner (toroidal) flow different within each of two subregions (to NW and to SE from the present descending slab). Also the curvature of the descending slab can be influenced by the slab heterogeneities due to variations in its thickness and viscosity (Morra *et al.* 2006).

Low velocity anomalies NW of the present slab between 70 and 110 km depth (see the temperature slice at the depth of 90 km) are interpreted by Martin *et al.* (2006) as a shallow asthenospheric upwelling associated with possible slab rollback. Also they mention partial melting as an additional contribution to the reduction of seismic velocities at these depths. The results of our assimilation
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Temperature misfit in the model of the descending lithospheric slab beneath the southeastern Carpathians. The misfit is defined as an integral difference between the temperature assimilated to any time $t \in [\text{present}, 22 \text{ Myr ago}]$ and that predicted by the forward model (11)–(17) to the same time assuming the assimilated temperature 22 Myr ago as the initial condition for the forward model. Solid and dashed curves present the misfits for the cases of temperature assimilation using the QRV and BAD methods, respectively.</figcaption></figure>

show that the descending slab is surrounded by a narrow border of hotter rocks at depths of 70–110 km (the temperature difference between slab and its surroundings is up to 500 K). Although we do not consider the effects of slab dehydration or partial melting in the modelling, the numerical results obtained support the hypothesis of dehydration of the descending lithosphere and its partial melting as the primary source of reduction of seismic velocities at these depths and probably deeper (see temperature slices at the depths of 130–220 km). Some areas of high temperature at depths below 280 km can be associated with present mantle upwelling in the region. The areas are not visible in the slides of the restored temperatures, because the upwelling was likely not active 22 Myr ago.

To test the accuracy of the QRV data assimilation, we employ the temperature and mantle flow restored to the time of 22 Myr ago as the initial condition for a model of the slab evolution forward in time (eqs 11–17), run the model to the present, and analyse the temperature residual (the difference between the temperatures assimilated and those predicted by the forward model). Fig. 6 (solid curve) presents the integral temperature residual as a function of time. The maximum temperature residual is about 30°, which is an evidence of rather accurate data assimilation based on the QRV method.

We compare the numerical results with that obtained by the backward advection of temperature (using the BAD method). The diffusion term in the heat equation is neglected. Fig. 5(c) (right column) presents slices of the temperature restored to 22 Myr ago. It is difficult (if even possible) to observe a slab-like feature in the slices at depths of 40–130 km. The cold (green, *ca.* 1000 K) slab was advected (with no diffusion) from the deeper mantle (340 km deep) to shallow levels (where temperature is about 900–1000 K), and therefore, the shape of the slab is indistinguishable in the shallow mantle. The shape of the slab is seen at depths of 170 and 220 km, where heat

advection becomes stronger than diffusion. Fig. 6 (dashed curve) shows that the maximum temperature residual is about 360°. Therefore, we have demonstrated here that the neglect of heat diffusion leads to an inaccurate restoration of mantle temperature, especially in the areas of low temperature and high viscosity in our model. The similar results for the BAD data assimilation have been obtained in the synthetic case study (see Figs 4e and h). The VAR method was not employed to assimilate the present temperature, because computations in this case become quite time-consuming due to the unavoidable need to smooth the solution and to filter temperature noise.

## 6 DISCUSSION AND CONCLUSION

The computational approach to assimilation of mantle related data proposed in the paper is based on the QRV method by Lattes & Lions (1969) and on the introduction into the backward heat equation of an additional term involving the product of a small regularization parameter and higher order temperature derivative. The QRV method is an alternative method for assimilation of data related to mantle dynamics. The method allows employing more sophisticated mathematical models (because it does not require derivation of an adjoint problem as in the VAR data assimilation) and hence expands the scope for applications in geodynamics (e.g. thermochemical convection, phase transformations in the mantle). It does not require that the desired accuracy of computations be directly related to the parameters of the numerical algorithm. However, the regularizing operators usually used in the QRV method enhance the order of the system of differential equations to be solved.

We compare the QRV method with the VAR and BAD methods in terms of solution stability, convergence and accuracy, time interval for data assimilation, analytical and algorithmic works, and computer performance (see Tables 1–3). The VAR data assimilation assumes that the direct and adjoint problems are constructed and solved iteratively forward in time. The structure of the adjoint problem is identical to the structure of the original problem, which considerably simplifies the numerical implementation. However, the VAR method imposes some requirements for the mathematical model (i.e. a derivation of the adjoint problem). Moreover, for an efficient numerical implementation of the VAR method, the error level of the computations must be adjusted to the parameters of the algorithm, and this complicates computations.

The BAD is the simplest method for data assimilation in models of mantle dynamics, because it does not require any additional work (neither analytical nor computational). The major difference between the BAD method and two other methods (VAR and QRV methods) is that the BAD method is by design expected to work (and hence can be used) *only* in advection-dominated heat flow. In the regions of high temperature/low mantle viscosity, where heat is transferred mainly by convective flow, the use of the BAD method is justified, and the results of numerical reconstructions can be considered to be satisfactory. Otherwise, in the regions of conduction-dominated heat flow (due to either high mantle viscosity or high conductivity of mantle rocks), the use of the BAD method cannot even guarantee any similarity of reconstructed structures. If mantle structures are diffused significantly, the remaining features of the structures can be only backward advected with the flow. We summarize the comparison between the methods for data assimilation in terms of a quality of numerical results in Table 2. The quality of the results is defined here as a relative (not absolute) measure of their accuracy. The results are good, satisfactory or poor compared
<table><caption>**Table 1.** Comparison of methods for data assimilation in models of mantle dynamics.</caption><thead><tr><th></th><th>QRV method</th><th>VAR method</th><th>BAD method</th></tr></thead><tbody><tr><th>Method</th><td>Solving the regularized<br>backward heat problem<br>with respect to parameter β</td><td>Iterative sequential solving of<br>the direct and adjoint heat<br>problems</td><td>Solving of heat advection<br>equation backward in time</td></tr><tr><th>Solution's stability</th><td>Stable for parameter β to<br>numerical errors (see text;<br>also in<sup>a</sup>) and conditionally<br>stable for parameter β to<br>arbitrarily assigned initial<br>conditions (numerically<sup>b</sup>)</td><td>Conditionally stable to<br>numerical errors depending<br>on the number of iterations<br>(theoretically<sup>c</sup>) and<br>unstable to arbitrarily<br>assigned initial conditions<br>(numerically<sup>d</sup>)</td><td>Stable theoretically<br>and numerically</td></tr><tr><th>Solution's<br>convergence</th><td>Numerical solution to the<br>regularized backward heat<br>problem converges to the<br>solution of the backward<br>heat problem in the special<br>class of admissible<br>solutions<sup>e</sup></td><td>Numerical solution converges<br>to the exact solution in the<br>Hilbert space<sup>f</sup></td><td>Not applied</td></tr><tr><th>Solution's accuracy<sup>g</sup></th><td>Acceptable accuracy for both<br>synthetic and geophysical<br>data</td><td>High accuracy for synthetic<br>data.</td><td>Low accuracy for both<br>synthetic and geophysical<br>data in<br>conduction-dominated<br>mantle flow</td></tr><tr><th>Time interval for data<br>assimilation<sup>h</sup></th><td>Limited by the characteristic<br>thermal diffusion time</td><td>Limited by the characteristic<br>thermal diffusion time and<br>the accuracy of the<br>numerical solution</td><td>No specific time limitation;<br>depends on mantle flow<br>intensity</td></tr><tr><th>Analytical work</th><td>Choice of the regularizing<br>operator</td><td>Derivation of the adjoint<br>problem</td><td>No additional analytical work</td></tr><tr><th>Algorithmic work</th><td>New solver for the<br>regularized equation<br>should be developed</td><td>No new solver should be<br>developed</td><td>No new solver should be<br>developed</td></tr></tbody></table>

<table><caption>**Table 2.** Quality of the numerical results obtained by different methods for data assimilation.</caption><thead><tr><th rowspan="2">Quality</th><th colspan="2">Synthetic data</th><th colspan="2">Geophysical data</th></tr><tr><th>Advection-dominated<br>regime</th><th>Diffusion-dominated<br>region</th><th>Advection-dominated<br>regime</th><th>Diffusion-dominated<br>region</th></tr></thead><tbody><tr><td>Good</td><td>VAR</td><td>VAR</td><td>—</td><td>—</td></tr><tr><td>Satisfactory</td><td>QRV, BAD</td><td>QRV</td><td>QRV, BAD</td><td>QRV</td></tr><tr><td>Poor</td><td>—</td><td>BAD</td><td>—</td><td>BAD</td></tr></tbody></table>

to other methods for data assimilation considered in this study. The numerical results of the reconstructions for both synthetic (see Figs 2–4) and geophysical case studies (Figs 5 and 6) show the comparison quantitatively.

There is a physical limitation of the restoration of mantle structures. If a thermal feature created, let us say, a several hundred million

years ago has completely diffused away by the present, it is impossible to restore the feature, which was more prominent in the past. The time to which a present thermal structure in the upper mantle can be restored should be restricted by the characteristic thermal diffusion time, the time when the temperatures of the evolved structure and the ambient mantle are nearly indistinguishable. In fact, the time

<table><caption>**Table 3.** Performance of data assimilation methods.</caption><thead><tr><th rowspan="2">Method</th><th colspan="2">CPU time (circa, in s)</th><th rowspan="2">Total</th></tr><tr><th>Solving the Stokes problem<br>using 50 × 50 × 50 finite<br>elements</th><th>Solving the backward heat<br>problem using 148 × 148 × 148<br>finite difference mesh</th></tr></thead><tbody><tr><td>BAD</td><td>180</td><td>2.5</td><td>182.5</td></tr><tr><td>QRV</td><td>100–180</td><td>3</td><td>103–183</td></tr><tr><td>VAR</td><td>360</td><td>1.5 n</td><td>360 + 1.5 n</td></tr></tbody></table>

<footnote><sup>a</sup>Lattes & Lions 1969; <sup>b</sup>See Fig. 3 and relevant text in the paper; <sup>c</sup>Ismail-Zadeh et al. 2004a; <sup>d</sup>Ismail-Zadeh et al. 2006; <sup>e</sup>Tikhonov & Arsenin, 1977; <sup>f</sup>Tikhonov & Samarskii 1990; <sup>g</sup>See Table 2; <sup>h</sup>See text for details.</footnote>
duration for which data assimilation methods can provide reasonable
results is much shorter than the characteristic thermal diffusion time
interval. The time interval for the VAR data assimilation depends
strongly on smoothness of the input data and the solution. The time
interval for the BAD data assimilation depends on the intensity of
mantle convection: it is short for conduction-dominated heat trans-
fer and becomes longer for advection-dominated heat flow. We note
that in the absence of thermal diffusion the backwards advection
of a low-density fluid in the gravity field will finally yield a uni-
formly stratified, inverted density structure, where the low-density
fluid overlain by a dense fluid spreads across the lower boundary
of the model domain to form a horizontal layer. Once the layer is
formed, information about the evolution of the low-density fluid will
be lost, and hence any forward modelling will be useless, because
no information on initial conditions will be available (Ismail-Zadeh
et al. 2001; Kaus and Podladchikov 2001).

The QRV method can provide stable results within the character-
istic thermal diffusion time interval. However, the length of the time
interval for QRV data assimilation depends on several factors. Let
us explain this by the example of heat conduction eq. (1). Assume
that the solution to the backward heat conduction equation with the
boundary conditions (2) and the initial condition T(t = t*, x) =
T*(x) satisfies the inequality ||∂⁴T/∂x⁴|| ≤ L_d at any time t. This
strong additional requirement can be considered as the requirement
of sufficient smoothness of the solution and initial data. Considering
the regularized backward heat conduction eq. (5) with the boundary
conditions (6)–(7) and the input temperature T_β(t = t*, x) = T_β*(x)
and assuming that ||T_β* - T*|| ≤ δ, Samarskii & Vabishchevich
(2004) estimated the temperature misfit between the solution
T(t, x) to the backward heat conduction problem and the so-
lution T_β(t, x) to the regularized backward heat conduction
equation:

$$
\| T(t, x) - T_{\beta}(t, x) \| \le \tilde{C} \delta \exp[\beta^{-1/2}(t^* - t)] + \beta L_d t, \quad 0 \le t \le t^*,
$$
where constant $\tilde{C}$ is determined from the *a priori* known parame-
ters of the backward heat conduction problem. For the given reg-
ularization parameter $\beta$, errors in the input data $\delta$ and smoothness
parameter $L_d$, it is possible to evaluate the time interval $0 \le t \le$
$t^*$ of data assimilation for which the temperature misfit would not
exceed a prescribed value.

Computer performance of the data assimilation methods can be
estimated by a comparison of CPU times for solving the inverse
problem of thermal convection. Table 3 lists the CPU times required
to perform one time step computations on 16 processors. The CPU
time for the case of the QRV method is presented for a given regu-
larization parameter $\beta$; in general, the total CPU time increases by
a factor of $\aleph$, where $\aleph$ is the number of runs required to determine
the optimal regularization parameter $\beta^*$. The numerical solution of
the Stokes problem (by the conjugate gradient method) is the most
time consuming calculation: it takes about 180 s to reach a high
accuracy in computations of the velocity potential. The reduction
in the CPU time for the QRV method is attained by employing the
velocity potential computed at $\beta_i$ as an initial guess function for
the conjugate gradient method to compute the vector potential at
$\beta_{i+1}$. An application of the VAR method requires to compute the
Stokes problem twice to determine the ‘advected’ and ‘true’ veloci-
ties (Ismail-Zadeh et al. 2004a). The CPU time required to compute
the backward heat problem using the TVD solver (Section A2) is
about 3 s in the case of the QRV method and 2.5 s in the case of the
BAD method. For the VAR case, the CPU time required to solve the
direct and adjoint heat problems by the semi-Lagrangian method is

1.5 × *n*, where *n* is the number of iterations in the gradient method used to minimize the cost functional (see eq. 5 in Ismail-Zadeh *et al.* 2004a).

Apart from the errors associated with the numerical modelling
(model, discretization and iteration errors), there are at least two
sources of errors in data assimilation: (i) data misfit associated with
the uncertainties in the present temperature distribution and/or in
the surface movements and (ii) errors associated with the uncer-
tainties in initial and boundary conditions. Since there are no direct
measurements of mantle temperatures, the temperatures can be es-
timated indirectly from either seismic wave (and their anomalies),
geochemical analysis or through the extrapolation of surface heat
flow observations. Many models of mantle temperature are based on
the conversion of seismic tomography data into temperature. Mean-
while, a seismic tomography image of the Earth's mantle is a model
indeed and incorporates its own model errors. Another sources of
uncertainty comes from the choice of mantle compositions in the
modelling of mantle temperature from the seismic velocities and
from the geodetic measurements of horizontal and vertical move-
ments. Therefore, if the present mantle temperature (and movement)
models are biased, the information can be improperly propagated to
the geological past.

The conditions at the boundaries of the model domain we used in
data assimilation are, of course, an approximation to the real temper-
ature, heat flux and movements, which are practically unknown and,
what is more important, may change over time at these boundaries.
The results of data assimilation will hence depend on the model
boundary conditions. Moreover, errors associated with the knowl-
edge of the temperature (or heat flux) evolution at the core–mantle
boundary or of the regional horizontal surface movements are an-
other essential component of errors, which can be propagated into
the past during the data assimilation.

Sensitivity analysis assists in understanding the stability of the
model solution to small perturbations in input variables or parame-
ters (Cacuci 2003). For instance, if we consider mantle temperature
in the past as a solution to the backward model, what will be its
variation if there is some perturbation in the inputs of the model
(e.g. present temperature data)? Despite the theoretical proof of the
stability of the solution to the regularized backward heat equation
(Samarskii & Vabishchevich 2004), we have performed a few tests for
stability of the numerical results with respect to small perturbations
in the initial data. Namely, the initial temperature has been perturbed
randomly by 0.5–2 per cent and then assimilated to the past. The
misfit between the restored temperatures related to the disturbed and
undisturbed initial temperature is limited, and hence the solution of
the problem is stable. Therefore, the proposed numerical approach
for assimilation of data related to thermoconvective movements in
the mantle, which is based on the QRV method, shows stability of
the numerical results with respect to errors in the initial data. This
method allows reconstructing the thermal state and dynamics of the
mantle in the past given its state and dynamics at the present time.

We have applied the QRV data assimilation to two models of
mantle dynamics: mantle plumes and lithospheric slab evolution.
For both models the present temperature and mantle flow have been
assimilated to the geological past, and the prominent features of
mantle structures have been recovered. The results of the data as-
simulation have been compared to that obtained by the VAR and
BAD methods. Based on the numerical results and the comparison
of the numerical methods, we arrive at the conclusion that the QRV
method is a highly promising approach to assimilation of geophys-
ical data related to mantle dynamics.
## ACKNOWLEDGMENTS

The authors are very grateful to Thorsten Becker, Michael Gurnis, Satoru Honda and Boris Kaus for their careful reviews that improved an initial version of the manuscript and to Olivier Talagrand for his fruitful discussion on data assimilation. We are thankful to Michael Martin and Friedemann Wenzel for the seismic tomography data on the Vrancea region made available for this research. This work was supported by the German Research Foundation (project no. DFG-Wi-687/18-2), the Russian Foundation for Basic Research (project No. 05-01-00098), and the Cooperation Program of the French Ministry of Research. This article has been written during the stay of AIZ in the Earthquake Research Institute of the University of Tokyo, and the author is grateful to S. Honda and H. Kawakatsu for their hospitality.

## REFERENCES

Alekseev, A.K. & Navon, I.M., 2001. The analysis of an ill-posed problem using multiscale resolution and second order adjoint techniques, *Comput. Meth. Appl. Mech. Eng.*, **190**, 1937–1953.

Axelsson, O., 1996. *Iterative Solution Methods*, Cambridge Univ. Press, Cambridge, 667 p.

Boussinesq, J., 1903. *Theorie Analytique de la Chaleur*, Vol. 2, 172 pp., Gauthier-Villars, Paris.

Bubnov, V.A., 1976. Wave concepts in the theory of heat, *Int. J. Heat Mass Transfer*, **19**, 175–184.

Bubnov, V.A., 1981. Remarks on wave solutions of the nonlinear heat conduction equation, *J. Eng. Phys. Thermophys.*, **40**(5), 907–913.

Bunge, H.-P., Hagelberg, C.R. & Travis, B.J., 2003. Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography, *Geophys. J. Int.*, **152**, 280–301.

Bunge, H.-P., Richards, M.A. & Baumgardner, J.R., 2002. Mantle circulation models with sequential data-assimilation: inferring present-day mantle structure from plate motion histories, *Phil. Trans. Roy. Soc. A*, **360**, 2545–2567.

Bunge, H.-P., Richards, M.A., Lithgow-Bertelloni, C., Baumgardner, J.R., Grand, S.P. & Romanowicz, B., 1998. Time scales and heterogeneous structure in geodynamic earth models, *Science*, **280**, 91–95.

Busse, F.H. et al. 1993. 3D convection at infinite Prandtl number in Cartesian geometry—a benchmark comparison, *Geophys. Astrophys. Fluid Dynamics*, **7**, 39–59.

Cacuci, D.G., 2003. *Sensitivity and Uncertainty Analysis. Volume I: Theory*, Chapman & Hall/CRC, Boca Raton, 285 p.

Cattaneo, C., 1958. Sur une forme de l'équation de la chaleur elinant le paradox d'une propagation instantanée, *Comptes Rendus*, **247**, 431–433.

Chandrasekhar, S., 1961. *Hydrodynamic and Hydromagnetic Stability*, Oxford Univ. Press, Oxford, 654 p.

Conrad, C.P. & Gurnis, M., 2003. Seismic tomography, surface uplift, and the breakup of Gondwanaland: integrating mantle convection backwards in time, *Geochem. Geophys. Geosys.*, **4**(3), doi:10.1029/2001GC000299.

Courant, R., Friedrichs, K.O. & Lewy, H., 1928. Uber die partiellen differenzen-gleichungen der mathematisches physik, *Math. Ann.*, **100**, 32–74.

Courant, R., Isaacson, E. & Rees, M., 1952. On the solution of nonlinear hyperbolic differential equations by finite differences, *Commun. Pure Appl. Math.*, **5**, 243–255.

Csontos, L., Nagymarosy, A., Horvath, F. & Kovac, M., 1992. Tertiary evolution of the intra-Carpathian area: a model, *Tectonophysics*, **208**, 221–241.

Davaille, A. & Vatteville, J., 2005. On the transient nature of mantle plumes, *Geophys. Res. Lett.*, **32**, L14 309, doi:10.1029/2005GL023029.

Ewing, R.E. & Wang, H., 2001. A summary of numerical methods for time-dependent advection-dominated partial differential equations, *J. Comput. Appl. Math.*, **128**, 423–445.

Falcone, M. & Ferretti, R., 1998. Convergence analysis for a class of high-order semi-Lagrangian advection schemes, *SIAM J. Numer. Anal.*, **35**, 909–940.

Fletcher, R. & Reeves, C.M., 1964. Function minimization by conjugate gradients, *Comp. J.*, **7**, 149–154.

Forte, A.M. & Mitrovica, J. X., 2001. Deep-mantle high-viscosity flow and thermochemical structure inferred from seismic and geodynamic data, *Nature*, **410**, 1049–1056.

Goes, S., Govers, R. & Vacher, P., 2000. Shallow mantle temperatures under Europe from P and S wave tomography, *J. Geophys. Res.*, **105**, 11 153–11 169.

Golub, G.H. & Van Loan, C.F., 1989. *Matrix Computations*, 2nd ed., 476 pp., Johns Hopkins, Baltimore.

Green, D.H. & Falloon, T.J., 1998. Pyrolite: a Ringwood concept and its current expression, in *The Earth's Mantle*, pp. 311–378, ed. Jackson, I., Cambridge University Press, Cambridge.

Hadamard, J., 1902. Sur les problèmes aux dérivées partielles et leur signification physique, *Princeton Univ. Bull.*, **15**, 49–52.

Harten, A., 1983. High resolution schemes for hyperbolic conservation laws, *J. Comp. Phys.*, **49**, 357–393.

Hier-Majumder, C.A., Belanger, E., DeRosier, S., Yuen, D.A. & Vincent, A.P. (2005). Data assimilation for plume models, *Nonlinear Proc. Geophys.*, **12**, 257–267.

Ismail-Zadeh, A.T., 1999. A time-inverse problem of gravitational instability: numerical approach and applications, *Geophys. Res. Abstr.*, **1**, 61.

Ismail-Zadeh, A.T., Korotkii, A. I. & Tsepelev, I.A., 2003a. Numerical approach to solving problems of slow viscous flow backwards in time, in *Computational Fluid and Solid Mechanics*, pp. 938–941, ed. K.J. Bathe., Elsevier Science, Amsterdam.

Ismail-Zadeh, A., Mueller, B. & Schubert, G., 2005a. Three-dimensional modeling of present-day tectonic stress beneath the earthquake-prone southeastern Carpathians based on integrated analysis of seismic, heat flow, and gravity observations, *Phys. Earth Planet. Inter.*, **149**, 81–98.

Ismail-Zadeh, A.T., Mueller, B. & Wenzel, F., 2005b. Modelling of descending slab evolution beneath the SE-Carpathians: implications for seismicity, in *Perspectives in Modern Seismology, Lecture Notes in Earth Sciences, Volume 105*, pp. 205–226, ed. F. Wenzel, Springer-Verlag, Heidelberg.

Ismail-Zadeh, A.T., Talbot, C.J. & Volozh, Y.A., 2001a. Dynamic restoration of profiles across diapiric salt structures: numerical approach and its applications, *Tectonophysics*, **337**, 21–36.

Ismail-Zadeh, A.T., Korotkii, A.I., Naimark, B.M. & Tsepelev, I.A., 2001b. Numerical simulation of three-dimensional viscous flows with gravitational and thermal effects, *Comp. Math. & Math. Phys.*, **41**, 1399–1415.

Ismail-Zadeh, A.T., Korotkii, A.I., Naimark, B.M. & Tsepelev, I.A., 2003b. Three-dimensional numerical simulation of the inverse problem of thermal convection, *Comp. Math. & Math. Phys.*, **43**, 587–599.

Ismail-Zadeh, A., Schubert, G., Tsepelev, I. & Korotkii, A., 2004a. Inverse problem of thermal convection: numerical approach and application to mantle plume restoration, *Phys. Earth Planet. Inter.*, **145**, 99–114.

Ismail-Zadeh, A., Schubert, G., Tsepelev, I. & Korotkii, A., 2006. Three-dimensional forward and backward numerical modeling of mantle plume evolution: effects of thermal diffusion, *J. Geophys. Res.*, **111**, B06 401, doi:10.1029/2005JB003782.

Ismail-Zadeh, A.T., Tsepelev, I.A., Talbot, C.J. & Korotkii, A.I., 2004b. Three-dimensional forward and backward modelling of diapirism: numerical approach and its applicability to the evolution of salt structures in the Pricaspian basin, *Tectonophysics*, **387**, 81–103.

Ismail-Zadeh, A.T., Korotkii, A.I., Naimark, B.M., Suetov, A.P. & Tsepelev, I.A., 1998. Implementation of a three-dimensional hydrodynamic model for evolution of sedimentary basins, *Comp. Math. Math. Phys.*, **38**(7), 1138–1151.

Kalnay, E., 2003. *Atmospheric Modeling, Data Assimilation and Predictability*, Cambridge Univ. Press, Cambridge, 341 p.

Katsura, T. et al. 2004. Olivine-wadsleyite transition in the system (Mg, Fe)2SiO4, *J. Geophys. Res.*, **109**, B02 209, doi:10.1029/2003JB002438.

Kaus, B.J.P. & Podladchikov, Y.Y., 2001. Forward and reverse modeling of the three-dimensional viscous Rayleigh-Taylor instability, *Geophys. Res. Lett.*, **28**(6), 1095–1098.
Kirsch, A., 1996. *An Introduction to the Mathematical Theory of Inverse Problems*, Springer-Verlag, New York, 282 pp.

Korotkii, A. I., Tsepelev, I. A., Ismail-Zadeh, A. T. & Naimark, B. M., 2002. Three-dimensional backward modeling in problems of Rayleigh-Taylor instability. *Proceedings of the Ural State University*, **22**(4), 96–104.

Lattes, R. & Lions, J.L., 1969. *The Method of Quasi-Reversibility: Applications to Partial Differential Equations*, Elsevier, New York, 388 p.

Martin, M., Ritter, J.R.R. & the CALIXTO working group, 2005. High-resolution teleseismic body-wave tomography beneath SE Romania – I. Implications for three-dimensional versus one-dimensional crustal correction strategies with a new crustal velocity model, *Geophys. J. Int.*, **162**, 448–460.

Martin, M., Wenzel, F. & the CALIXTO working group, 2006. High-resolution teleseismic body wave tomography beneath SE-Romania – II. Imaging of a slab detachment scenario, *Geophys. J. Int.*, **164**, 579–595.

Martin, M., Achauer, U., Kissling, E., Mocanu, V., Musacchio, G., Radulian, M., Wenzel, F. & CALIXTO working group, 2001. First results from the tomographic experiment CALIXTO'99 in Romania, *Geophys. Res. Abst.*, **3**, SE1.02.

McDonald, A. & Bates, J. R., 1987. Improving the estimate of the departure point isosi on in a two-time level semi-Lagrangian and semi-implicit scheme, *Mon. Weather Rev.*, **115**, 737–739.

Montelli, R., Nolet, G., Dahlen, F.A., Masters, G., Engdahl, E.R. & Hung, S.-H., 2004. Finite-frequency tomography reveals a variety of plumes in the mantle, *Science*, **303**, 338–343.

Moore, W. B., Schubert, G. & Tackley, P., 1998. Three-dimensional simulations of plume–lithosphere interaction at the Hawaiian Swell, *Science*, **279**, 1008–1011.

Morra, G., Regenauer-Lieb, K. & Giardini, D., 2006. Curvature of oceanic arcs, *Geology*, **34**, 877–880.

Morse, P.M. & Feshbach, H., 1953. *Methods of Theoretical Physics*, McGraw-Hill, New York, 865 p.

Ribe, N.M. & Christensen, U., 1994. Three-dimensional modeling of plume-lithosphere interaction, *J. Geophys. Res.*, **99**, 669–682.

Samarskii, A.A. & Vabishchevich, P.N., 1995. *Computational Heat Transfer. Vol. 1. Mathematical Modelling*, John Wiley & Sons, New York, 370 p.

Samarskii, A.A. & Vabishchevich, P.N. 1998. Nonlinear monotone schemes for the advection equation, *Doklady/Transactions of the Russian Academy of Sciences*, **361**(1), 21–23.

Samarskii, A.A. & Vabishchevich, P.N., 2004. *Numerical Methods for Solving Inverse Problems of Mathematical Physics*, URSS, Moscow, 478 p. (in Russian).

Samarskii, A.A., Vabishchevich, P.N. & Vasiliev, VI., 1997. Iterative solution of a retrospective inverse problem of heat conduction, *Math. Modeling*, **9**(5), 119–127.

Schubert, G., Turcotte, D.L. & Olson, P., 2001. *Mantle Convection in the Earth and Planets*, Cambridge University Press, Cambridge, 940 p.

Staniforth, A. & Coté, J., 1991. Semi-Lagrangian integration schemes for atmospheric models—a Review, *Mon. Weather Rev.*, **119**(9), 2206–2223.

Steinberger, B. & O'Connell, R.J., 1998. Advection of plumes in mantle flow: implications for hotspot motion, mantle viscosity and plume distribution, *Geophys. J. Int.*, **132**, 412–434.

Sweby, P.K., 1984. High resolution schemes using flux limiters for hyperbolic conservation laws, *SIAM J. Numer. Anal.*, **21**(5), 995–1001.

Tikhonov, A.N., 1963, Solution of incorrectly formulated problems and the regularization method, *Doklady Akademii Nauk SSSR*, **151**, 501–504 (Engl. transl.: *Soviet Math. Dokl.*, **4**, 1035–1038, 1963).

Tikhonov, A.N. & Arsenin, V.Y., 1977. *Solution of Ill-Posed Problems*, Winston, Washington, DC., 258 p.

Tikhonov, A.N. & Samarskii, A.A., 1990. *Equations of Mathematical Physics*, Dover Publications, New York, 765 pp.

Vernotte, P., 1958. Les paradoxes de la theorie continue de l'equation de la chaleur, *Comptes Rendus*, **246**, 3154–3155.

Wang, Y. & Hutter, K., 2001. Comparison of numerical methods with respect to convectively dominated problems, *Int. J. Numer. Meth. Fluids*, **37**, 721–745.

Yu, N., Imatani, S. & Inoue, T., 2004. Characteristics of temperature field due to pulsed heat input calculated by non-Fourier heat conduction hypothesis, *JSME Int. J.*, **Series A**, **47**(4), 574–580.

## APPENDIX A: NUMERICAL METHODS EMPLOYED

In the model domain $\Omega$ we assume a uniform spatial partition $x_i^{m_i} = m_i h_i$ with the gridpoints $\omega_{m_1 m_2 m_3} = (x_1^{m_1}, x_2^{m_2}, x_3^{m_3})$, where $m_i = 0, 1, 2, \dots, n_i$, $h_i = l_i/n_i$ and $i = 1, 2, 3$. The spatial partition defines the computational grid. The diameter of the partition of the time interval is $\tau$. Let $T_{ijk}^n$ denote the value of the corresponding grid function at the gridpoint $\omega_{ijk}$ at the time $t_n \in [0, \hat{\tau}]$, where the lower case subscripts $ijk$ denote the gridpoints and the upper case superscript $n$ indicates the time step.

### A1. The numerical method for solving the heat problem

The characteristic-based semi-Lagrangian method (Courant 1952; Staniforth & Coté 1991) is used to calculate the convective derivatives of the heat eq. (15). It accounts for the Lagrangian nature of the advection process but, at the same time, it allows computations on a fixed grid. We rewrite the heat equation in the following form

$$ DT/Dt = \nabla^2 T + f, \quad DT/Dt = \partial T/\partial t + \mathbf{u} \cdot \nabla T. \quad (A1) $$

The aim of such a splitting is to solve the first equation on the characteristics of the second equation. This method has been used in advection-diffusion systems due to two useful properties of the approximations: (i) a relatively large time step may be used in a numerical simulation, and (ii) it is stable and accurate for arbitrary relations between the time and space steps (e.g. Ewing & Wang 2001). Moreover, the implementation of this method with a high-order interpolation of the space variables yields a minimum error in the variance. In particular, such an approach is used in meteorology, where the time step must be large to ensure computational efficiency (e.g. Staniforth & Coté 1991).

Eqs (A1) are approximated by finite differences in the following form

$$ \frac{T_{ijk}^{n+1} - T_{d}^{n}}{\tau} = \nabla^2 \frac{T_{ijk}^{n+1} + T_{ijk}^{n}}{2} + \frac{f_{ijk}^{n+1} + f_{ijk}^{n}}{2}, \quad (A2) $$

$$ Dz/Dt = \mathbf{u}(t, \mathbf{z}), \quad \mathbf{z}(t_{n+1}) = \mathbf{z}_a, \quad (A3) $$

where $T_d^n$ is the temperature at the point $z_d$. The point $z_d$ is obtained by solving eq. (A3) backward in time with the final condition $z_a$, which should coincide with the corresponding gridpoint $\omega_{ijk}$ at $t = t_{n+1}$. A solution to (A3) can be obtained by solving the following system of non-linear equations by an implicit method (the number of equations is equal to the number of gridpoints):

$$ \begin{aligned} z_d &= z_a - y_k, & y_{k+1} &= \tau \mathbf{u}(t_n, z_a - 0.5y_k), \\ y_0 &= \tau \mathbf{u}(t_n, z_a), & k &= 0, 1, 2, \dots \end{aligned} \quad (A4) $$

It can also be solved using the explicit predictor-corrector method

$$ \mathbf{z}^* = \mathbf{z}_a - \tau \mathbf{u}(t_n, \mathbf{z}_a), \quad \mathbf{z}_d = \mathbf{z}_a - \tau \mathbf{u}(t_n, \mathbf{z}^*). \quad (A5) $$

The point $\mathbf{z}^*$ does not necessarily coincide with a gridpoint, and the velocity at this point can be obtained by using the interpolation of the velocities at the adjacent gridpoints. The value of $T_d^n$ at the time $t = t_n$ and at the point $z_d$ can also be obtained by interpolation.

The total error of the method is estimated to be $O(\tau^2 + h^2 + \tau^s + \tau^{-1}h^{1+s})$ and is not monotonic with respect to the time step $\tau$, where
s is the order of integration of eq. (A3) backward in time, and q is the interpolation order (McDonald & Bates 1987; Falcone & Ferretti 1998). For example, s = 2 for the predictor-corrector method (A5), and s = 4 for the Runge-Kutta method. If cubic polynomials are used for interpolation, then q = 3; for linear interpolation q = 1.

A solution to (A4) can be obtained in 3-4 iterations, if Newton's method is used to solve the set of the non-linear equations and the Courant-Friedrichs-Lewy condition $\tau \|\partial u/\partial x\| < 1$ is satisfied (Courant et al. 1928). This condition guarantees that the trajectories of the characteristics do not intersect at one time step. The procedure of solving the characteristic equation forward and backward in time is unconditionally stable. Method (A5) is easier to implement, but it is inferior to method (A4) in terms of accuracy.

The 3-D spatial discrete operator associated with the diffusion term in eq. (A2) is split into 1-D operators as $\nabla^2 \approx \Delta_1 + \Delta_2 + \Delta_3$, and the latter operators are approximated by the central differences:

$$ \Delta_1 T_{ijk}^{n} = \frac{T_{i+1jk}^{n} - 2T_{ijk}^{n} + T_{i-1jk}^{n}}{h_1^2}, \quad i = 1, 2, \dots, n_1 - 1. \quad (A6) $$

At the boundary gridpoints $i=0$ and $i=n_1$, an approximation for $\Delta_1$ is obtained from (A6) with regard for the boundary conditions (16). Expressions for $\Delta_2$ and $\Delta_3$ are determined similarly. The set of difference equations for the approximation of the heat eq. (15) on a uniform rectangular grid has the form:

$$ T_{ijk}^{+} = T[t_n, \omega_{ijk} - \tau u(t_n, z_d)], \quad (A7) $$

$$ T_{ijk}^{*} = T_{ijk}^{+} + 1.5\tau \Delta_{1}(T_{ijk}^{*} + T_{ijk}^{+}) + 1.5\tau(f_{ijk}^{n+1} + f_{ijk}^{n}), \quad (A8) $$

$$ T_{ijk}^{**} = T_{ijk}^{+} + 1.5\tau \Delta_{2}(T_{ijk}^{**} + T_{ijk}^{+}), \quad (A9) $$

$$ T_{ijk}^{***} = T_{ijk}^{+} + 1.5\tau \Delta_{3}(T_{ijk}^{***} + T_{ijk}^{+}), \quad (A10) $$

$$ T_{ijk}^{n+1} = (T_{ijk}^{*} + T_{ijk}^{**} + T_{ijk}^{***})/3. \quad (A11) $$

In the numerical implementation of this scheme, $3(n_1 n_2 n_3)$ eqs (A4) or (A5) and $3(n_1 n_2 + n_1 n_3 + n_2 n_3)$ independent sets of linear algebraic eqs (A8)-(A10) with tridiagonal (diagonally dominant) matrices should be solved. To determine $T_{ijk}^{+}$, the velocity and temperature should be interpolated at the point $z_d$. Equations (A8)-(A10) can be solved independently, and hence the numerical code is designed for multiprocessor computer using the method of tridiagonal matrix factorization (e.g. Axelsson 1996).

## A2. The numerical method for solving the regularized backward heat problem

We employ the total variation diminishing (TVD) method (Harten 1983) to solve the regularized backward heat problem (19)-(22). Notice that when oscillations (e.g. due to non-smoothness of the solution or jumps in physical parameters) arise, the numerical solution will have larger total variation of temperature (that is, the sum of the variations of temperature over the whole computational domain $TV^n = \sum_i |T_{i+1jk}^n - T_{ijk}^n| + \sum_j |T_{ij+1k}^n - T_{ijk}^n| + \sum_k |T_{ijk+1}^n - T_{ijk}^n|$ will increase with oscillations). The TVD method is designed to yield well-resolved, non-oscillatory discontinuities by enforcing that the numerical schemes generate solutions with non-increasing total variations of temperature in time (that is $TV^{n+1} \le T V^n$), thus no spurious numerical oscillations are generated (Ewing & Wang 2001). The TVD method can describe convection problems with large temperature gradients very well, because it is at most first-order accurate

at local temperature extrema (Wang & Hutter 2001). Note that the characteristic-based semi-Lagrangian method described in Section A1 can be also used to solve the problem.

Consider initially an approximation of the advection term of Eq. (19):

$$ \Xi_1 = u_{1\beta} \partial T / \partial x_1 \approx (F_{x_1}^{+} - F_{x_1}^{-}) / h_1, \quad (A12) $$

$$ F_{x_1}^{+} = 0.5u_{1\beta,ijk} (T_{i+1/2jk}^{+} + T_{i+1/2jk}^{-}) - 0.5|u_{1\beta,ijk}| (T_{i+1/2jk}^{+} - T_{i+1/2jk}^{-}), \quad (A13) $$

$$ F_{x_1}^{-} = 0.5u_{1\beta,ijk} (T_{i-1/2jk}^{+} + T_{i-1/2jk}^{-}) - 0.5|u_{1\beta,ijk}| (T_{i-1/2jk}^{+} - T_{i-1/2jk}^{-}), \quad (A14) $$

$$ \begin{aligned} T_{i+1/2jk}^{-} &= T_{ijk} + 0.5\Upsilon(\xi_i)(T_{i+1jk} - T_{ijk}), \\ T_{i+1/2jk}^{+} &= T_{i+1jk} - 0.5\Upsilon(\xi_{i+1})(T_{i+2jk} - T_{i+1jk}), \end{aligned} \quad (A15) $$

$$ \begin{aligned} T_{i-1/2jk}^{-} &= T_{i-1jk} + 0.5\Upsilon(\xi_{i-1})(T_{ijk} - T_{i-1jk}), \\ T_{i-1/2jk}^{+} &= T_{ijk} - 0.5\Upsilon(\xi_{i})(T_{i+1jk} - T_{ijk}), \end{aligned} \quad (A16) $$

$$ \begin{aligned} \xi_i &= (T_{ijk} - T_{i-1jk}) / (T_{i+1jk} - T_{ijk}), \\ \Upsilon(\xi) &= \max\{0, \min\{1, 2\xi\}, \min\{\xi, 2\}\}, \end{aligned} \quad (A17) $$

where $\Upsilon(\xi)$ is a 'superbee' limiter (Sweby 1984). Expressions for $\Xi_2$ and $\Xi_3$ are determined similarly. The solution based on the TVD method gives a second-order accurate solution (Wang & Hutter 2001). Since the formula (A17) can generate logical difficulties in the case of $T_{ijk} = T_{i-1jk} = T_{i+1jk}$, the following alternative representation of (A17) is used in our computations:

$$ \begin{aligned} \Upsilon(\xi_i)(A) &= L(A, B) = \Upsilon(1/\xi_i)(B), \\ A &= T_{i+1jk} - T_{ijk}, \quad B = T_{ijk} - T_{i-1jk}, \end{aligned} \quad (A18) $$

$$ L(A, B) = 0.5(\text{sign}(A) + \text{sign}(B)) \max\{\min\{2|A|, |B|\}, \min\{|A|, 2|B|\}\}. \quad (A19) $$

This representation of the limiter $\Upsilon$ has an explicit symmetric form compared to (A15)-(A17). The TVD numerical scheme was tested using known solutions to simple advection equations and also compared to another TVD numerical scheme by Samarskii & Vabishchevich (1998).

The 3-D spatial discrete operator associated with the diffusion term in Eq. (19) is split into 1-D operators as $\nabla^2 \approx \Delta_1 + \Delta_2 + \Delta_3$, and the latter operators are approximated by the central differences (A6) as described in Section A1. The regularization term in Eq. (19) is approximated by the differences $\Lambda \approx \Theta_1 + \Theta_2 + \Theta_3$ and

$$ \Theta_1(T_{ijk}) = \frac{T_{i+2jk} - 4T_{i+1jk} + 6T_{ijk} - 4T_{i-1jk} + T_{i-2jk}}{h_1^4}, \quad i = 2, 3, \dots, n_1 - 2. \quad (A20) $$
The difference function $\Theta_1$ is approximated with regard for the boundary conditions (20) and (21) at the gridpoints $i = 0, 1$ and $i = n_1 - 1, n_1$. The expressions for $\Theta_2$ and $\Theta_3$ are determined similarly.

The system of difference equations for the approximation of the regularized backward heat equation on a uniform rectangular grid has the form

$$ \frac{T_{ijk}^* - T_{ijk}^n}{3\tau} + \beta \frac{\Theta_1(T_{ijk}^* - T_{ijk}^n)}{3\tau} = \Xi_1 T_{ijk}^n - \Delta_1 T_{ijk}^n, \quad (A21) $$

$$ \frac{T_{ijk}^{**} - T_{ijk}^{n}}{3\tau} + \beta \frac{\Theta_2(T_{ijk}^{**} - T_{ijk}^{n})}{3\tau} = \Xi_2 T_{ijk}^{n} - \Delta_2 T_{ijk}^{n} - f_{ijk}^{n}, \quad (A22) $$

$$ \frac{T_{ijk}^{***} - T_{ijk}^{n}}{3\tau} + \beta \frac{\Theta_3(T_{ijk}^{***} - T_{ijk}^{n})}{3\tau} = \Xi_3 T_{ijk}^{n} - \Delta_3 T_{ijk}^{n}, \quad (A23) $$

$$ T_{ijk}^{n+1} = (T_{ijk}^{*} + T_{ijk}^{**} + T_{ijk}^{***})/3. \quad (A24) $$

The total error of the numerical method is $O(\tau + h^2)$, and the iterations are stable at $\beta > \tau^2/16$. Considering the independence of eqs (A21)–(A23), they can be solved on a parallel computer using the method of five-diagonal matrix factorization (e.g. Axelsson 1996).