# Identification of a stratigraphic model with seismic constraints

Jacques Blum, Gabriela Dobranszky, Robert Eymard, Roland Masson

► To cite this version:

Jacques Blum, Gabriela Dobranszky, Robert Eymard, Roland Masson. Identification of a stratigraphic model with seismic constraints. Inverse Problems, 2006, 22, pp.1207-1225. 10.1088/0266-5611/22/4/006 . inria-00414400

HAL Id: inria-00414400

https://inria.hal.science/inria-00414400

Submitted on 8 Sep 2009

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
# Identification of a stratigraphic model with seismic constraints

Jacques Blum<fnref n="1" />, Gabriela Dobranszky<fnref n="2" />, Robert Eymard<fnref n="3" />
and Roland Masson<fnref n="2" />

E-mail: roland.masson@ifp.fr

## Abstract

Stratigraphic models are used in oil exploration in order to determine the geometry of sedimentary layers and to characterize the geological properties of the sediments. They simulate the infill of a sedimentary basin at large scales of space and time, governed by the interaction of accommodation change, sediment supply and sediment transport. One of the main difficulties in the use of such models is that all their input parameters have to be identified to fit well-log and seismic data. A basic issue is to compute the rate of accommodation change in order to fit the sediment thicknesses measured between each successive seismic marker. This leads to the solution to an inverse problem for which existence is proved in the single lithology case. Numerical algorithms are proposed and studied based on the construction of a reduced model and Newton-Krylov methods.

## 1. Introduction

The concepts of sequence stratigraphy have revolutionized our understanding of the formation and infill of sedimentary basins. The concepts essentially state that the geometry and the rock properties of sedimentary layers result from the interaction of the rate of subsidence, the rate of sea level change, the sediment supply at river entries, and the climatic environment governing the transport laws of the sediments.

Many forward numerical models have been developed to study the effects of these controlling parameters on the distributions of the rock properties and the stratal geometry of sedimentary basins (see, e.g. [TH89, AH89, TS94]). Among these, dynamic slope models

<footnote n="1">¹ Laboratoire Jean-Alexandre Dieudonné, Université de Nice Sophia Antipolis,
Faculté des Sciences, Parc Valrose, 06108 Nice Cedex 02, France</footnote>

<footnote n="2">² Institut Français du Pétrole, 1 et 4 avenue de Bois Préau 92852 Rueil Malmaison Cedex, France</footnote>

<footnote n="3">³ Département de Mathématiques, Université de Marne la Vallée, 5 boulevard Descartes,
Champs sur Marne, F-77454, Marne-la-Vallée, Cedex 2, France</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** 2D section of a basin at time *t* with the input parameters *a* and *φ*, and the unknowns *h* and *b*: basement *z* = *H*(*x*, *t*), sea level *z* = *H*<sub>m</sub>(*t*), accommodation *a*(*x*, *t*) = *H*<sub>m</sub>(*t*) − *H*(*x*, *t*), sediment supply *φ*, sediment thickness *h*(*x*, *t*), and bathymetry *b*(*x*, *t*) = *a*(*x*, *t*) − *h*(*x*, *t*).</figcaption></figure>

[AH89, TS94], based on diffusive transport, have already proved their efficiency in several petroleum applications [G97, GJD98, GJ99].

In order to describe the basic features of such models, let us introduce some notation. A basin is considered as a domain of $\mathbb{R}^{d+1}$ with $d=1$ for two-dimensional (2D) sections of a basin, and $d=2$ for three-dimensional (3D) basin models. The projection of the basin on a reference horizontal plane is assumed to be a fixed subdomain of $\mathbb{R}^d$ independent of time, denoted by $\Omega$, and defining the horizontal extension of the basin. We shall denote by $x \in \Omega$ the horizontal coordinates, by $z$ the vertical coordinate above the reference horizontal plane and by $t$ the time variable.

The sea level is denoted by $z = H_m(t)$ and the basement of the basin is the moving surface $z = H(x, t)$, $x \in \Omega$, accounting for the tectonic displacements, which are assumed to be vertical. The accommodation $a = H_m - H$ is the distance between the basement and the sea level, governing the available space for sedimentation in the basin. This is the first input parameter of the stratigraphic model.

The main unknowns of the model are the sediment thickness $h$ above the basement and the bathymetry $b = a - h$, functions of $x \in \Omega$ and time $t$. The second input parameter of the model is the boundary flux $\varphi$ defined at the boundary $\partial\Omega$ of $\Omega$ and accounting for the sediment supply at river entries of the basin.

Figure 1 exhibits these first two input parameters and unknowns of the model on a 2D section of a basin. For the sake of simplicity, it is assumed throughout the paper that the sediments have a null porosity, meaning that we do not consider their compaction. In that case, the stratigraphic model accounts for the conservation of the sediment thickness $h$ given a flux of sediments proportional to the gradient of the topography $h+H$. This diffusive transport law results from an averaging over several processes such as river transport, creep, slumps and small slides. It is valid only at the large geological scales in space and time used in oil exploration with typical discretization steps of 1 to 10 km and 10 000 to 100 000 years, and typical simulation ranges of 100 km and 10 million years (My).

The third input parameter of the model is the *diffusion coefficient* which depends on the climatic environment at the time of deposition. It is modelled as a nonlinear function of the bathymetry in order to distinguish between low marine and high continental diffusive transport.

The stratigraphic model is extended to take into account sediments defined as a mixture of L basic lithologies such as sand or shale, leading to the diffusive multi-lithology stratigraphic model introduced in [R92] (see also [G97, GJ99, EGGMT04]). The mixture is characterized
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Seismic section (x, z) of a sedimentary basin at present time with seismic markers delimitating three sedimentary layers and a well log interpreted in terms of ranges of concentrations and paleo-bathymetries.</figcaption></figure>

by concentrations of each lithology inside the basin denoted by $c_i \ge 0, i = 1, \dots, L$, such that $\sum_{i=1}^L c_i = 1$. In that case, the flux of each lithology is proportional to the gradient of the topography with a diffusion coefficient for each lithology. The unknowns of the model are the sediment thickness $h$, the bathymetry $b$ and the concentrations $c_i, i = 1, \dots, L$, and the model accounts for the conservation equations of each lithology $i = 1, \dots, L$. We shall consider in the following the simulation of a basin layer through the time interval or sequence $(0, T)$. The bathymetry and concentrations are assumed to be known at initial time $t = 0$. The main outputs of the stratigraphic model are the geometry of the layer at simulation time $T$ defined by the bathymetry $b(x, T), x \in \Omega$, and the sediment thickness $h(x, T), x \in \Omega$. The concentration as well as the paleo-bathymetry (bathymetry at deposition time) distributions inside the layer characterize the rock properties of the layer.

None of the input parameters, the accommodation, the sediment supply and the diffusion coefficients, can be directly obtained from available data. The predictivity of the stratigraphic model crucially depends on the identification, on well-log and seismic data, of these input parameters.

Seismic data record the reflections in the basin of mechanical waves emitted at the ground surface. These signals are processed in order to determine the position of seismic reflectors imaging the stratal geometry of the basin at the present time. As shown in the 2D seismic section of figure 2, one can delineate, on such an image, seismic markers representing topographic surfaces of the basin at successive times up to the vertical tectonic displacements which have occurred between the time of deposition and the present time. Since these tectonic displacements are unknown, the only data which can be used to fit the stratigraphic model parameters are the sediment thicknesses of the layers between successive markers. For a given layer, say the second layer in figure 2, the sediment thickness of the seismic data is denoted by $h^d$ and defined on $\Omega$.

Well-log data are obtained by measuring physical properties such as radioactivity or resistivity within a given well. In order to be used in a stratigraphic model, they are interpreted in terms of concentrations and paleo-bathymetry ranges of the sediments crossed by the well. This is illustrated in figure 2 which exhibits interpreted well-log data (each band stands for a class of sediments which correspond to ranges of concentrations and paleo-bathymetries). The identification of the input parameters of the stratigraphic model (accommodation $a$, sediment supply $\varphi$ and diffusion coefficients) from the seismic and well-log data is done layer-by-layer with a piecewise linear or piecewise constant interpolation in time of the parameters on each sequence.
In order to cope with the difficulty of simultaneously identifying all these input parameters, our methodology is based on the construction of a reduced model. This reduced model is a stationary approximation of the stratigraphic model, assuming a piecewise constant sedimentation rate $∂_t h = h^d/T$ on each sequence (0, T) corresponding to a given layer. Since h is no longer an unknown, it is possible to solve the stratigraphic model with the new unknown a rather than h. Moreover, one can show that the resulting model reduces to a time-independent problem.

The main advantage of the reduced model is the computation of the accommodation a as the solution of the reduced model and no longer as the solution of an inverse problem. In our methodology described in [DMAC04, GD05], this reduced model is further used to identify the remaining sediment supply and diffusion coefficient input parameters.

Nevertheless, the construction of the reduced model relies on the strong assumption that the sedimentation rate is piecewise constant in time on each sequence. This assumption does not take into account the non-stationary transition between two successive layers, and is known to be a poor approximation near the shoreline.

In order to improve the sedimentation rate interpolation in time, the objective of this paper is to study the identification of the accommodation such that the layer thickness $h|_{t=T}$ obtained from the stratigraphic model fits the seismic data $h^d$ for a given layer on the sequence (0, T). The solution of this inverse problem will provide a more accurate sedimentation rate which will later be used by the reduced model at a finer time scale to improve the identification of the remaining sediment supply and diffusion coefficient parameters. It is shown for $L=1$ that this inverse problem has a solution and is well posed in the linear case. Next, solution algorithms based on Krylov-Newton methods preconditioned by the linearized reduced model are proposed and tested.

In order to clarify the presentation, the stratigraphic model, the inverse problem, the reduced model and the solution algorithms are first derived in the case of a single lithology ($L=1$). In such a case, the stratigraphic model reduces to a nonlinear parabolic equation on the bathymetry unknown $b$, with the rate of accommodation change $∂_t a$ as a forcing source term. The extension to the multi-lithology model is postponed until section 5. In section 2, the stratigraphic model is introduced and the inverse problem is defined in the case of a single lithology model. The existence of a stable solution to the inverse problem is proved in section 3 for the single lithology model, and the uniqueness is obtained in the linear case. In section 4, the solution algorithms based on the construction of the reduced model and nonlinear Krylov methods are presented. The numerical results are exhibited in section 6, which illustrate the efficiency and the robustness of our approach with respect to the diffusion coefficients, and the time and space discretization parameters.

## 2. The stratigraphic model and the inverse problem

### 2.1. Single lithology diffusion model

It is known in sequence stratigraphy [AH89, TS94] that, at large scales of space and time, the evolution of a basin topography is well described by the solution of a nonlinear parabolic equation. This model accounts for the conservation of the sediment thickness $h$ given a flux of sediments proportional to the gradient of the topography $-\nabla(h+H)$ or equivalent to $\nabla b$. The diffusion coefficient denoted by $K(b)$ is assumed to be a non-negative function of the bathymetry. The boundary conditions prescribe the flux $\varphi$ at the boundary $∂Ω$. The input flux $\varphi < 0$ represents the sediment supply at river entries. The output flux $\varphi \ge 0$ usually either
vanishes sufficiently far away from the input boundary or has also to be identified as an input parameter of the model.

Let us consider the simulation of a stratigraphic layer on a given sequence (0, T), T > 0 with an initial bathymetry $b_0$. Given the input parameters $a, \varphi, b_0, K(b)$, the stratigraphic model finds the solution $h$ such that

$$ 
\begin{array}{lll}
\partial_t h + \operatorname{div}(K(b)\nabla b) = 0 & \text{on } \Omega \times (0, T), \\
K(b)\nabla b \cdot \mathbf{n} = \varphi & \text{on } \partial\Omega \times (0, T), \\
b = a - h & \text{on } \Omega \times (0, T), \\
b|_{t=0} = b_0 & \text{on } \Omega.
\end{array}
\quad (1) $$

The nonlinear dependence of the diffusion coefficient $K$ on the bathymetry $b$ accounts for the fact that the diffusive transport is usually much higher in continental ($b \le 0$) than in marine ($b > 0$) environments. It is typically given by

$$ K(b) = \begin{cases} K_m & \text{if } b \ge \beta, \\ \left(\frac{K_m}{K_c}\right)^{\frac{b}{2\beta}} (K_c K_m)^{\frac{1}{2}} & \text{if } -\beta < b < \beta, \\ K_c & \text{if } b \le -\beta, \end{cases} \quad (2) $$

with $\beta \ge 0$, and $K_c \ge K_m > 0$.

The model (1) can be viewed as a nonlinear parabolic equation for the bathymetry $b$ with Neumann boundary conditions $\varphi$, and the rate of accommodation change $\partial_t a$ as a forcing source term. This latter formulation will be used in the following analysis of the inverse problem defined below. Nevertheless, the formulation (1) will be more conveniently introduced in the reduced model in section 4.1, and in the extension to the multi-lithology model in section 5.

## 2.2. Definition of the inverse problem

Only the low frequencies of the accommodation change can be derived from the seismic data. In particular, with such data, one cannot image the high-frequency oscillations of the sea level $H_m$. The result is that the accommodation is parametrized on each sequence (0, T) by the following linear variation in time,

$$ a = b_0 + \frac{t}{T}(a_T - b_0), \quad (3) $$

where it has been assumed to fix the idea that $h|_{t=0} = 0$ i.e. $b_0 = a|_{t=0}$. This is always true up to a change of the coordinate system. The initial bathymetry $b_0$ is assumed to be known from the previous layer solution.

The rate of accommodation change $\partial_t a = \frac{a_T - b_0}{T}$ is identified such that the sediment thickness $h|_{t=T}$ at the end of the sequence fits the thickness measured on the seismic data denoted by $h^d$.

This leads to the definition of the following inverse problem: given $b_0, K(b), \varphi, b_w, h^d$, find $a_T$ and $b$ such that

$$ 
\left\{
\begin{array}{lll}
\partial_t b - \operatorname{div}(K(b)\nabla b) = \frac{a_T - b_0}{T} & \text{on } \Omega \times (0, T), \\
K(b)\nabla b \cdot \mathbf{n} = \varphi & \text{on } \partial\Omega \times (0, T), \\
b|_{t=0} = b_0 & \text{on } \Omega, \\
b|_{t=T} = a_T - h^d & \text{on } \Omega, \\
\frac{1}{|\omega_w|} \int_{\omega_w} b(x, T) \, dx = b_w & \text{on } \Omega.
\end{array}
\right.
\quad (4) $$

where the diffusion coefficient $K(b)$, and the sediment supply $\varphi$ are fixed input parameters of the stratigraphic model.
Due to the Neumann boundary condition $\varphi$ set on the stratigraphic model, it must be assumed that the condition on the data

$$ \int_{\Omega} \frac{h^d}{T} dx + \int_{\partial\Omega} \varphi d\sigma = 0, \qquad (5) $$

always holds in the following. This condition is necessary for the existence of a solution to (4). This also explains the equation $\frac{1}{|\omega_w|} \int_{\omega_w} b(x, T) dx = b_w$ in (4), which must be set to close the system (4) and expect its unicity. In this equation, the prescribed average bathymetry $b_w$ is obtained in practice from well-log data on a non-empty open subset $\omega_w \subset \Omega$ representing the given well sections of measure $|\omega_w|$.

**Remark 2.1.** In practice, the remaining sediment supply $\varphi$ and diffusion coefficient $K(b)$ input parameters are identified using the reduced model only (introduced later in section 4.1), and the well-log data interpreted in terms of paleo-bathymetries (see [DMAC04, GD05]).

Also, at the first layer, only the reduced model is used to compute the accommodation, since this model does not depend on the initial bathymetry $b_0$ which is unknown for the first layer.

In the case of erosion of a layer after its deposition, the sediment thickness $h^d$ is no longer known from the seismic data. In that case, the thickness of the eroded part of the layer has also to be identified using the reduced model and well-log data interpreted in terms of paleo-bathymetries (see [GD05]).

Condition (5) simply states that the volume of sediment in the layer defined by $h^d$ matches the volume of sediment brought by the sediment supply at the boundary $\varphi$. This is to be viewed as a constraint on the sediment supply $\varphi$ that will in practice always be satisfied from the identification procedure of $\varphi$ using the reduced model (see [DMAC04, GD05]).

The following section shows that the inverse problem (4) is well posed in the linear case for which the diffusion coefficient depends only on $x$ and $t$, and the existence of a stable solution is proved in the general nonlinear case.

# 3. Analysis of the single lithology model

## 3.1. Linear model

Let us consider the inverse problem (4) in the linear case, for which the diffusion coefficient $K$ is a given function in $L^\infty(\Omega \times (0, T))$ with $0 < K_m = \inf_{(x,t) \in \Omega \times (0,T)} K(x, t)$, $K_c = \sup_{(x,t) \in \Omega \times (0,T)} K(x, t) < +\infty$. Let us denote by $V$ the function space

$$ V = \{v \in L^2(0, T; H^1(\Omega)), \partial_t v \in L^2(0, T; H^{-1}(\Omega))\} $$

endowed with its usual norm denoted by $\|\cdot\|_V$. One can prove the following proposition.

**Proposition 3.1.** Assuming that $\varphi \in H^{-\frac{1}{2}}(\partial\Omega)$, $b_0$, $h^d \in L^2(\Omega)$, then, there exists a unique solution $(a_T, b) \in L^2(\Omega) \times V$ to the inverse problem (4). This solution satisfies the bound

$$ \|a_T\|_{L^2(\Omega)} + \|b\|_V \le C (\|b_0\|_{L^2(\Omega)} + \|h^d\|_{L^2(\Omega)} + \|\varphi\|_{H^{-\frac{1}{2}}(\partial\Omega)} + |b_w|), \quad (6) $$

with $C$ a constant independent of $h^d$, $b_0$, $\varphi$ and $b_w$.

**Proof.** Let us set

$$ L_0^2(\Omega) = \left\{ u \in L^2(\Omega) \middle| \int_{\Omega} u \, dx = 0 \right\}, \qquad (7) $$
and $a_T = b_0 + v - \gamma$ with $\gamma \in \mathbb{R}$, and $v \in L_0^2(\Omega)$. Let us also define $\tilde{b} = b - \hat{b} + \gamma \frac{t}{T}$, where $\hat{b}$ is the unique solution of (1) with $a = 0$.

The inverse problem (4) equivalently amounts to finding $(v, \tilde{b}, \gamma) \in L_0^2(\Omega) \times V \times \mathbb{R}$ such that

$$ \left\{ \begin{array}{ll} \partial_t \tilde{b} - \operatorname{div}(K \nabla \tilde{b}) = \frac{v}{T} & \text{on } \Omega \times (0, T), \\ K \nabla \tilde{b} \cdot \mathbf{n} = 0 & \text{on } \partial \Omega \times (0, T), \\ \tilde{b}|_{t=0} = 0 & \text{on } \Omega, \end{array} \right. \quad (8) $$

$$ \tilde{b}|_{t=T} = v + b_0 - \hat{b}|_{t=T} - h^d \quad \text{on } \Omega, \quad (9) $$

and $\gamma = -b_w + \frac{1}{|\omega_w|} \int_{\omega_w} (\hat{b} + \tilde{b})(x, T) dx$. One can show that the solution $\tilde{b}|_{t=T}$ of (8) for $v$ given, defines a linear bounded operator $\mathcal{L}$ from $L_0^2(\Omega)$ to $L_0^2(\Omega)$ which is a strict contraction in $L_0^2(\Omega)$. It results that equation (9) $(I - \mathcal{L})(v) = \hat{b}|_{t=T} + h^d - b_0$ has a unique solution in $L_0^2(\Omega)$, and that proposition 3.1 holds. $\square$

Let $(v_i)_{i \in \mathbb{N}}$, $v_i \in H^1(\Omega)$, $v_0 = |\Omega|^{-1/2}$, be the orthonormal basis of $L^2(\Omega)$ defined by the eigenvectors of the operator $-\Delta$ on the domain $\Omega$ with homogeneous Neumann boundary conditions. The eigenvalues are denoted by $\omega_i$, $i \in \mathbb{N}$, and satisfy $\omega_0 = 0$, $\omega_i \ge \omega_j > 0$ for all $i > j$, $i, j \in \mathbb{N}^*$.

Assuming furthermore that $K$ is a constant, the solution of (4) is explicit in the spectral basis $(v_i)_{i \in \mathbb{N}}$ as given by the following proposition.

**Proposition 3.2.** If $K$ is a constant, the solution of proposition 3.1 is defined by

$$ a_T - a_{T,0}v_0 = \sum_{i=1}^{\infty} \frac{1}{1 - \alpha(K\omega_i T)} (h_i^d + \alpha(K\omega_i T)T\varphi_i + (e^{-K\omega_i T} - \alpha(K\omega_i T))b_{0,i})v_i, \quad (10) $$

with $h_i^d = \int_{\Omega} h^d v_i \, dx$, $b_{0,i} = \int_{\Omega} b_0 v_i \, dx$, $\varphi_i = \int_{\partial\Omega} \varphi v_i \, d\sigma$, $\alpha(s) = \frac{1-e^{-s}}{s}$, and $a_{T,0}$ is computed from the equation $\frac{1}{|\omega_w|} \int_{\omega_w} (a_T - h^d) \, dx = b_w$.

If $h^d \in H^1(\Omega)$ the solution $a_T$ satisfies the estimates

$$ \left\{ \begin{array}{l} \|a_T - a_{T,0}v_0\|_{L^2(\Omega)} \le C \frac{1}{1-\alpha(K\omega_1 T)} \left( \|h^d\|_{L^2(\Omega)} + \frac{1}{KT\omega_1^{1/2}} \left( T\|\varphi\|_{H^{-1/2}(\partial\Omega)} + \|b_0\|_{H^{-1}(\Omega)} \right) \right), \\ \|\nabla a_T\|_{L^2(\Omega)} \le C \frac{1}{1-\alpha(K\omega_1 T)} \left( \|\nabla h^d\|_{L^2(\Omega)} + \frac{1}{KT} \left( T\|\varphi\|_{H^{-1/2}(\partial\Omega)} + \|b_0\|_{H^{-1}(\Omega)} \right) \right), \end{array} \right. \quad (11) $$

with $C$ a constant independent of any parameter.

### 3.2. Nonlinear model

It is assumed that the diffusion coefficient $K(b)$ is a Lipschitz function such that $0 < K_m \le K(b) \le K_c$ for all $b \in \mathbb{R}$.

**Proposition 3.3.** Assuming that $\varphi \in H^{-\frac{1}{2}}(\partial\Omega)$, $b_0$, $h^d \in L^2(\Omega)$, then, there exists at least one solution $(a_T, b) \in L^2(\Omega) \times V$ to the inverse problem (4). This solution satisfies the bound

$$ \|a_T\|_{L^2(\Omega)} + \|b\|_V \le C (\|b_0\|_{L^2(\Omega)} + \|h^d\|_{L^2(\Omega)} + \|\varphi\|_{H^{-\frac{1}{2}}(\partial\Omega)} + |b_w|), \quad (12) $$

with $C$ a constant independent of $h^d$, $b_0$, $\varphi$ and $b_w$.

**Proof.** Let us consider the application $\mathcal{P}$ from $L^2(\Omega \times (0, T))$ to $L^2(\Omega \times (0, T))$ such that $b = \mathcal{P}(\hat{b})$ is the solution of (4) for $K(x, t) = K(\hat{b}(x, t))$ for which existence, uniqueness and estimates in the function space $V$ are given by proposition 3.1. Since the injection from $V$ into $L^2(\Omega \times (0, T))$ is compact [Au63], the Schauder fixed point
theorem (see [GD03]) will conclude the proof once the continuity of $\mathcal{P}$ in $L^2(\Omega \times (0, T))$ is shown. Let us consider $\hat{b}_1, \hat{b}_2 \in L^2(\Omega \times (0, T))$, and $b_1 = \mathcal{P}(\hat{b}_1)$, $b_2 = \mathcal{P}(\hat{b}_2)$, and let us set $\tilde{b} = b_1 - b_2$, $\tilde{a}_T = a_T^1 - a_T^2$, $\tilde{K}(x, t) = \frac{1}{2}(K(\hat{b}_1(x, t)) + K(\hat{b}_2(x, t)))$, and $S \in L^2(0, T; H^{-1}(\Omega))$ such that $\langle S, u \rangle = \int_{\Omega} \frac{1}{2}(K(\hat{b}_2(x, t)) - K(\hat{b}_1(x, t)))\nabla(b_1+b_2) \cdot \nabla u \, dx$ for all $u \in L^2(0, T; H^1(\Omega))$. Then $\tilde{b}$ satisfy the following equation,

$$ 
\left\{
\begin{array}{ll}
\partial_t \tilde{b} - \operatorname{div}(\tilde{K} \nabla \tilde{b}) = \frac{\tilde{a}_T}{T} + S & \text{on } \Omega \times (0, T), \\
\tilde{K} \nabla \tilde{b} \cdot \mathbf{n} = 0 & \text{on } \partial\Omega \times (0, T), \\
\tilde{b}|_{t=0} = 0 & \text{on } \Omega, \\
\tilde{b}|_{t=T} = \tilde{a}_T & \text{on } \Omega, \\
\frac{1}{|\omega_w|} \int_{\omega_w} \tilde{b}(x, T) \, dx = 0 &
\end{array}
\right.
\qquad (13)
$$

which from the Lipschitz assumption on $K(b)$ together with proposition 3.1 (which extends to the case of a right-hand side $S$) proves the continuity of $\mathcal{P}$. $\square$

## 4. Preconditioned Newton-Krylov algorithms for the inverse problem

Let us define $\mathcal{H}: L^2(\Omega) \to L^2(\Omega)$ such that $\mathcal{H}(a_T)$ is the unique solution $h|_{t=T} = a_T - b|_{t=T}$ in $L^2(\Omega)$ of the stratigraphic model (1)–(3) for a given $a_T \in L^2(\Omega)$ (see [LSU68]). Our objective is to build an algorithm for the solution $a_T \in L^2(\Omega)$ of the nonlinear problem

$$ 
\left\{
\begin{array}{l}
\mathcal{H}(a_T) = h^d, \\
\frac{1}{|\omega_w|} \int_{\omega_w} a_T(x) \, dx = b_w + \frac{1}{|\omega_w|} \int_{\omega_w} h^d(x) \, dx,
\end{array}
\right.
\qquad (14)
$$

that should have good scalability properties with respect to the diffusion coefficients, and the space and time discretization parameters (mesh size and time step). This problem is highly nonlinear due to the usually large jump $\frac{K_c}{K_m}$ in (2) of the diffusion coefficient $K(b)$ at the shoreline $b=0$, and is ill-conditioned for a usually small diffusion coefficient $K_m$ as shown by (10) in the linear case. To cope with these difficulties, our method relies on the construction of a reduced model for the inverse problem (4) that will be used as a preconditioner of the GMRes solver in Newton-Krylov algorithms.

### 4.1. Reduced model and quasi-Newton algorithm

The construction of the reduced model is based on the assumption that the sediment thickness varies linearly in time, while the accommodation $a$ is considered as an unknown no longer constrained to satisfy (3). Thus, setting

$$ h = \frac{t}{T} h^d \quad \text{in (4)}, $$

we obtain a family of solutions $b$ such that

$$ \psi(b) = \phi + \alpha \quad \text{with} \quad \psi(u) := \int_{0}^{u} K(s) \, ds. \qquad (15) $$

In (15), $\alpha$ is an arbitrary continuous function depending only on time such that $\frac{1}{|\omega_w|} \int_{\omega_w} b(x, T) \, dx = b_w$, and $\phi \in H^1(\Omega)$ is a solution of $-\Delta\phi = \frac{h^d}{T}$ on $\Omega$, $\nabla\phi \cdot \mathbf{n} = \varphi$ on $\partial\Omega$.

The accommodation given by $a = \psi^{-1}(\phi + \alpha) + \frac{t}{T}h^d$ is no longer continuous between two successive sequences. This model does not take into account the non-stationary transition between two successive layers since the solution does not even depend on $b_0$. Hence, it is only used in the following to compute the final bathymetry $b_T$ and the final accommodation $a_T$, unique solutions of the following stationary problem which defines our reduced model.
Given $h^d$, $K(b)$, $\varphi$, $b_w$, find $(a_T, b_T) \in L^2(\Omega) \times H^1(\Omega)$ such that

$$
\begin{align}
-\Delta\phi &= \frac{h^d}{T} && \text{on } \Omega, \nonumber \\
\nabla\phi \cdot \mathbf{n} &= \varphi && \text{on } \partial\Omega, \nonumber \\
b_T &= \psi^{-1}(\phi) && \text{on } \Omega, \tag{16} \\
\frac{1}{|\omega_w|} \int_{\omega_w} b_T(x) \, dx &= b_w \nonumber \\
a_T &= b_T + h^d && \text{on } \Omega. \nonumber
\end{align}
$$

This reduced model is shown in [DMAC04] to be a good approximation of the inverse
model (4). It is used in [GD05] to identify the sediment supply and diffusion coefficients on
the well-log data interpreted in terms of paleo-bathymetries.

For our purpose, it can be used to solve the nonlinear system (14) by a quasi-Newton approach. First, the reduced model defines an initial value $a_T^0$ of the accommodation $a_T$. Then, for $k \ge 0$, given $a_T^k$, the stratigraphic model (1) computes the sediment thickness $\mathcal{H}(a_T^k)$, the bathymetry $b^k$ and the residual $r(a_T^k) = \mathcal{H}(a_T^k) - h^k$.

Let $\mathcal{C}(b)$ denote the linear bounded operator from $L_0^2(\Omega)$ defined by (7) to

$$
L_w^2(\Omega) = \left\{ \delta a \in L^2(\Omega) \middle| \int_{\omega_w} \delta a(x) \, dx = 0 \right\},
$$

such that $\mathcal{C}(b)r$ is the solution $p$ of the following linearized reduced model:

$$
\begin{align}
-\Delta\delta\phi &= -\frac{r}{T} && \text{on } \Omega, \nonumber \\
\nabla\delta\phi \cdot \mathbf{n} &= 0 && \text{on } \partial\Omega, \nonumber \\
\psi'(b)|_{t=T}\delta b_T &= \delta\phi && \text{on } \Omega, \tag{17} \\
\int_{\omega_w} \delta p(x) \, dx &= 0 \nonumber \\
p &= \delta b_T - r && \text{on } \Omega. \nonumber
\end{align}
$$

The quasi-Newton correction step is defined by $a_T^{k+1} = a_T^k + \mathcal{C}(b^k)r(a_T^k)$, and the algorithm is iterated until a stopping criterion for the $L^2$ norm of the residual $r(a_T^k)$ is fulfilled.

In the case of a constant diffusion coefficient $K$, using the notation of proposition 11, the convergence rate in $L^2$ norm is equal to $\rho = \left(1 - \left(1 + \frac{1}{K\omega_1 T}\right)\right)(1 - \alpha(K\omega_1 T)) > 0$ bounded by $\frac{1}{2}$ independently of the diffusion coefficient $K$. By contrast, a simple fixed point algorithm $a_T^{k+1} = a_T^k - r(a_T^k) + \frac{1}{|\omega_w|} \int_{\omega_w} r(a_T^k) dx$ has a convergence rate $\rho = \frac{1-e^{-K\omega_1 T}}{K\omega_1 T}$ which deteriorates rapidly for small diffusion coefficients.

Nevertheless, in the nonlinear case, $\mathcal{C}(b^k)r(a_T^k)$ is not guaranteed to be a descent direction and the algorithm can lack robustness even using a relaxation of the quasi-Newton step. A better strategy developed in the next section is to use a Newton-Krylov algorithm (see [BS90, BS94]) preconditioned by the linearized reduced model $\mathcal{C}(b)$.

4.2. Preconditioned Newton–Krylov algorithms

Let $(b, a_T)$ be the current solution of the stratigraphic model (1)-(3). The Jacobian $J(a_T)$ of $\mathcal{H}(a_T)$ is the linear operator from $L_w^2(\Omega)$ to $L_0^2(\Omega)$ such that for any $q \in L_w^2(\Omega)J(a_T)q = q - \delta b|_{t=T}$ where $\delta b$ is the unique solution of the parabolic equation

$$
\begin{align}
\partial_t \delta b - \Delta(\psi'(b)\delta b) &= \frac{q}{T} && \text{on } \Omega \times (0, T), \nonumber \\
\nabla \psi'(b) \delta b \cdot \mathbf{n} &= 0 && \text{on } \partial\Omega \times (0, T), \nonumber \\
\delta b|_{t=0} &= 0 && \text{on } \Omega. \tag{18}
\end{align}
$$
The Newton step $p \in L_w^2(\Omega)$ is the solution of the equation

$$J(a_T)p = -r(a_T),$$

where $r(a_T) = \mathcal{H}(a_T) - h^d$ is the residual. This equation is solved using a GMRes linear iterative solver with the right preconditioner $\mathcal{C}(b)$. The matrix vector product is computed by the finite difference formula

$$J(a_T)v \sim \frac{\mathcal{H}(a_T + \epsilon v) - \mathcal{H}(a_T)}{\epsilon}$$

for a suitable choice of $\epsilon$ discussed in [BS90] and related to the expected accuracy of the stratigraphic model simulation code. To improve the robustness of Newton's methods a globalization strategy must be used. In this paper, we compare a line search method and a trust region dogleg algorithm. The merit function is the usual least-squares function $f(a_T) = \frac{1}{2}\|r(a_T)\|_{L^2(\Omega)}^2$ defined in $L_w^2(\Omega)$. Its gradient $\nabla f$ is the function of $L_w^2(\Omega)$ defined by

$$(J(a_T)v, r(a_T))_{L^2(\Omega)} = (\nabla f, v)_{L^2(\Omega)}$$

for all $v$ in $L_w^2(\Omega)$, i.e. the orthogonal projection on $L_w^2(\Omega)$ of the gradient of the merit function $f$ when defined on the full space $L^2(\Omega)$. This latter gradient is computed using the adjoint state of the stratigraphic model.

### 4.2.1. Newton-Krylov algorithm with backtracking line search globalization.
The line search algorithm chooses at each iteration $k$ a descent direction $p^k$ and searches along this direction a step length $\alpha^k$ to move from the current iterate $a_T^k$ to a new iterate $a_T^{k+1}$ with a lower merit function value. The new iterate is given by $a_T^{k+1} = a_T^k + \alpha^k p^k$.

In our case, $p^k$ is an inexact Newton step which approximates the solution of the Newton equation $J(a_T^k)p^k = -r(a_T^k)$ using a preconditioned GMRes algorithm. The step length $\alpha^k$ is an approximate minimizer of the one-dimensional function $\phi$ defined by $\phi(\alpha) = f(a_T^k + \alpha p^k), \alpha > 0$, which must satisfy the Armijo condition [NW99] $f(a_T^k + \alpha^k) \le f(a_T^k) + \delta\alpha^k(\nabla f(a_T^k), p^k)_{L^2(\Omega)}$ with $0 < \delta < 1$. This condition ensures a sufficient decrease of the merit function $f$.

**Step length selection.** We use a procedure based on a quadratic interpolation of the function $\phi$ using $\phi(0), \phi(1)$ and $\phi'(0)$. Note that the Armijo condition can be rewritten as follows: $\phi(\alpha^k) \le \phi(0) + \delta\alpha^k\phi'(0)$. This procedure can be seen as an enhancement of the backtracking approach [NW99].

<table><thead><tr><th>Procedure</th><th>Step length selection</th></tr></thead><tbody><tr><td>choose $\alpha, \delta$;</td><td></td></tr><tr><td>while $\phi(\alpha) > \phi(0) + \delta\alpha\phi'(0)$ do</td><td></td></tr><tr><td>&nbsp;&nbsp;&nbsp;&nbsp;$\alpha \leftarrow -\frac{\phi'(0)\alpha^2}{2[\phi(\alpha_i) - \phi(0) - \phi'(0)\alpha]}$;</td><td></td></tr><tr><td>end</td><td></td></tr></tbody></table>

### 4.2.2. Newton-Krylov algorithm with dogleg trust region globalization.
The trust region strategy defines a region around the current iterate $a_T^k$ within which the quadratic model $m^k$ of
the merit function $f(a_T^k)$ is supposed to be a 'good' representation of the function, and then chooses the step direction $p^k$ to be an approximate minimizer of the quadratic model in the trust region. For the least-squares merit function $f(a_T^k) = \frac{1}{2}\|r(a_T^k)\|_2^2$, the quadratic model $m^k$ is defined by

$$m^k(p) = \frac{1}{2}\|r(a_T^k) + J(a_T^k)p\|_2^2, \quad (19)$$

and the step $p^k$ is an approximate solution of the sub-problem

$$\min_{\|p\| \le \Delta^k} m^k(p) \qquad (20)$$

where $\Delta^k$ is the radius of the trust region. The exact solution of (20) is denoted by $p^*(\Delta^k)$.

**Dogleg direction.** The dogleg method (see [NW99]) finds an approximate solution to (20) by replacing the curved trajectory $p^*(\Delta^k)$ with a path consisting of two line segments. The first segment runs from the origin to the Cauchy point $p_C^k = -\min(\frac{\Delta^k}{\|\nabla f(a_T^k)\|}, \frac{\|\nabla f(a_T^k)\|^2}{\|J(a_T^k)\|\|\nabla f(a_T^k)\|^2})\nabla f(a_T^k)$, and the second segment runs from the Cauchy point to the inexact Newton direction denoted by $p_N^k \sim -J(a_T^k)^{-1}r(a_T^k)$. Then, the path $p^k$ is obtained using the *dogleg direction* procedure stated below.

**Trust region update.** The choice of the trust region radius at each iteration $k$ is based on the agreement between the merit function $f(a_T^k)$ and the quadratic model $m^k$ at the previous iteration, which is evaluated by the ratio $\rho_k$

$$\rho_k = \frac{\|r(a_T^k)\|^2 - \|r(a_T^k + p^k)\|^2}{\|r(a_T^k)\|^2 - \|r(a_T^k) + J(a_T^k)p^k\|^2} \quad (21)$$

between the actual reduction and the predicted reduction. In order to define the success level (failure, success, high success) of the current dogleg step and to update the trust region radius according to the success level, two thresholds parameters $\bar{\rho}_1$ and $\bar{\rho}_2$ such that $\bar{\rho}_1 < \bar{\rho}_2$ and three multiplicative factors $\beta_1, \beta_2$ and $\beta_3$ such that $\beta_1 < \beta_2 < \beta_3$ are introduced. The new radius is computed using the *trust region radius update* procedure.

<table><tr><td><b>Procedure</b></td><td>Dogleg direction</td></tr><tr><td>given</td><td>Δ<sup>k</sup>;</td></tr><tr><td>calculate</td><td>p<sub>N</sub><sup>k</sup>;</td></tr><tr><td><b>if</b></td><td><b>||p<sub>N</sub><sup>k</sup>|| ≤ Δ<sup>k</sup></b> <b>then</b></td></tr><tr><td>|</td><td>p<sup>k</sup> &larr; p<sub>N</sub><sup>k</sup></td></tr><tr><td><b>else</b></td><td></td></tr><tr><td>calculate</td><td>p<sub>C</sub><sup>k</sup>;</td></tr><tr><td><b>if</b></td><td><b>||p<sub>C</sub><sup>k</sup>|| = Δ<sup>k</sup></b> <b>then</b></td></tr><tr><td>|</td><td>p<sup>k</sup> &larr; p<sub>C</sub><sup>k</sup></td></tr><tr><td><b>else</b></td><td></td></tr><tr><td>||</td><td>p<sup>k</sup> &larr; p<sub>C</sub><sup>k</sup> + τ(p<sub>N</sub><sup>k</sup> - p<sub>C</sub><sup>k</sup>) where<br/>τ is the largest value in [0, 1]<br/>such that ||p<sup>k</sup>|| ≤ Δ<sup>k</sup>;</td></tr><tr><td><b>end</b></td><td></td></tr><tr><td><b>end</b></td><td></td></tr></table>

<table><tr><td><b>Procedure</b></td><td>Trust region radius update</td></tr><tr><td>given</td><td>Δ̄, Δ<sup>k</sup>, ρ<sub>k</sub>;</td></tr><tr><td>given</td><td>ρ̄<sub>1</sub>, ρ̄<sub>2</sub>;</td></tr><tr><td>given</td><td>β<sub>1</sub>, β<sub>2</sub>, β<sub>3</sub>;</td></tr><tr><td><b>if</b></td><td><b>ρ<sub>k</sub> &lt; ρ̄<sub>1</sub></b> <b>then</b></td></tr><tr><td>|</td><td>Δ<sub>k+1</sub> = β<sub>1</sub>Δ<sup>k</sup>;</td></tr><tr><td><b>else</b></td><td></td></tr><tr><td><b>if</b></td><td><b>ρ<sub>k</sub> &lt; ρ̄<sub>2</sub></b> <b>then</b></td></tr><tr><td>|</td><td>Δ<sub>k+1</sub> = min(β<sub>2</sub>Δ<sup>k</sup>, Δ̄);</td></tr><tr><td><b>else</b></td><td></td></tr><tr><td>|</td><td>Δ<sub>k+1</sub> = min(β<sub>3</sub>Δ<sup>k</sup>, Δ̄);</td></tr><tr><td><b>end</b></td><td></td></tr><tr><td><b>end</b></td><td></td></tr></table>
# 5. Extension to the multi-lithology stratigraphic model

## 5.1. Multi-lithology model

The sediments are modelled as a mixture of L lithologies characterized by their grain size population and considered as incompressible materials of constant grain density and null porosity (the compaction is not considered for the sake of the simplicity of the presentation).

The basin sediments are described by their L concentrations $c_i \ge 0$ in lithology i defined on the domain $\mathcal{B} = \{(x, \zeta, t)| (x, t) \in \Omega \times (0, T), \zeta \in (-\infty, h(x, t))\}$ and satisfying $\sum_{i=1,...,L} c_i = 1$. The sediments transported at the surface of the basin, i.e. deposited in the case of sedimentation or passing through the surface in the case of erosion, are described by their L surface concentrations $c_i^s \ge 0$ defined on $\Omega \times (0, T)$ and also verifying $\sum_{i=1,...,L} c_i^s = 1$.

Since the compaction is not considered, and the sediment transport only occurs at the surface of the basin, the composition inside the basin remains constant: $\partial_t c_i = 0$ on $\mathcal{B}$. The evolution of $c_i$ is only governed by the boundary condition at the top of the basin stating that $c_i|_{\zeta=h} = c_i^s$ in the case of sedimentation ($\partial_t h > 0$). This yields

$$ \left\{ \begin{array}{ll} \partial_t c_i = 0 & \text{on } \mathcal{B}, \\ c_i|_{\zeta=h} = c_i^s & \text{on } \{(x, t) \in \Omega \times (0, T) | \partial_t h(x, t) > 0\}, \\ c_i|_{t=0} = c_i^0 & \text{on } \{(x, \zeta) | x \in \Omega, \zeta \in (-\infty, h|_{t=0}(x))\}, \end{array} \right. \qquad (22) $$

where $c_i^0$, $i = 1, ..., L$ are the initial concentrations of the initial basin.

These linear hyperbolic equations are coupled to the conservation equations of each lithology on the domain $\Omega \times (0, T)$ accounting for the conservation of the lithology fraction $h_i(x, t) = \int_0^{h(x,t)} c_i(x, \zeta, t) d\zeta$:

$$ \partial_t h_i + \text{div} \mathbf{f}_i = 0, $$

on $\Omega \times (0, T)$ with $\partial_t h_i = c_i|_{\zeta=h} \partial_t h$ since $\partial_t c_i = 0$ on $\mathcal{B}$. The surficial transport model considered here is the multi-lithology diffusive model introduced in [R92] for which the flux in lithology i is taken proportional to the surface concentration $c_i^s$ with a diffusion coefficient depending on the lithology i: $\mathbf{f}_i = -c_i^s r_i K(b) \nabla b$ with $r_i > 0$. A total flux $\sum_{j=1}^L \mathbf{f}_j \cdot \mathbf{n} = \varphi \in H^{-\frac{1}{2}}(\partial \Omega)$ is prescribed on $\partial \Omega$, and input fractional fluxes $\mu_{i,e} \in L^\infty(\partial \Omega_−)$ on the subset $\partial \Omega_− = \{x \in \partial \Omega / \varphi(x) < 0\}$, with $\mu_{i,e} \ge 0$ for all $i = 1, ..., L$, and $\sum_{j=1}^L \mu_{i,e} = 1$.

Finally, we end up with the following multi-lithology diffusive model: given $\varphi$, $\mu_{i,e}$, $a$, $r_i$, $K(b)$, find $c_i^s$, $c_i$ and $b$ such that for all $i = 1, ..., L$

$$ \left\{ \begin{array}{ll} c_i|_{\zeta=h} \partial_t h + \text{div}(-c_i^s r_i K(b) \nabla b) = 0, & \text{on } \Omega \times (0, T), \\ \sum_{j=1}^L c_j^s = 1 & \text{on } \Omega \times (0, T), \\ \frac{r_i c_j^s}{\sum_{j=1}^L r_j c_j^s} = \mu_{i,e} & \text{on } \partial \Omega_− \times (0, T), \\ \sum_{j=1}^L \mathbf{f}_j \cdot \mathbf{n} = \varphi & \text{on } \partial \Omega \times (0, T), \\ b = a - h & \text{on } \Omega \times (0, T), \\ b|_{t=0} = b_0 & \text{on } \Omega. \end{array} \right. \qquad (23) $$

$$ \left\{ \begin{array}{ll} \partial_t c_i = 0 & \text{on } \mathcal{B}, \\ c_i|_{\zeta=h} = c_i^s & \text{on } \{(x, t) \in \Omega \times (0, T) | \partial_t h(x, t) > 0\}, \\ c_i|_{t=0} = c_i^0 & \text{on } \{(x, \zeta) | x \in \Omega, \zeta \in (-\infty, h|_{t=0}(x))\}. \end{array} \right. \qquad (24) $$
In the case of equal diffusion coefficients for all the lithologies, the existence of a unique weak solution has been proved for this system of PDEs (see [EGGM03, GM04]). A more detailed presentation of the model including an additional constraint on the erosion rate can be found in [EGGMT04].

The inverse problem is still defined by (14) for which the operator $\mathcal{H}(a_T) = a_T - b|_{t=T}$ is given by the solution $b$ of (23), (24) with the accommodation $a$ satisfying (3). As before, it is always assumed that the necessary condition (5) for the existence of a solution to (14) holds.

## 5.2. Adjoint state

The computation of the gradient $\nabla f$ of the merit function $f$ is done using the adjoint state of the multi-lithology stratigraphic model leading to an optimal complexity. The main difficulty of adjoint methods is to circumvent the full storage of the stratigraphic model solution, resulting in our case in a 4D storage of the concentration variables $c_i$, while the multi-lithology stratigraphic model output is usually restricted to a 3D storage ($c_i|_{t=T}, c_i^s, h$).

To avoid this, we first note that equations (22) do not depend on the accommodation $a$. The result is that only the 3D storage of the adjoint state variables related to equations (23) on $\Omega \times (0, T)$ is required. Second, in order to avoid the 4D storage of the solutions $c_i$ of the stratigraphic model, they are recomputed backwards in time as the solutions of the following hyperbolic equations for each $i = 1, \dots, L$,

$$ \left\{ \begin{array}{ll} -\partial_t c_i = 0 & \text{on } \mathcal{B}, \\ c_i|_{\zeta=h} = c_i^h & \text{on } \{(x, t) \in \Omega \times (0, T) | \partial_t h(x, t) < 0\}, \\ c_i|_{t=T} = c_i^T & \text{on } \{(x, \zeta) | x \in \Omega, \zeta \in (-\infty, h|_{t=T}(x))\}, \end{array} \right. \quad (25) $$

only requiring the 3D storage of $c_i^h = c_i|_{\zeta=h}, c_i^T = c_i|_{t=T}$ and $h$. Note that, in practice, this is the adjoint state of the discrete system which is computed following these ideas (see [GD05] for details).

## 5.3. Reduced model

The reduced model is still obtained assuming that the sediment thickness varies linearly in time: $h = \frac{t}{T}h^d$ and considering the accommodation as an unknown in (23), (24).

It is assumed in addition, in case of erosion of the layer ($h^d < 0$), that the initial concentration of the basin does not depend on $\zeta$. Let us denote by $d_i$ the concentration at the top of the basin $c_i|_{\zeta=h}$, and by $d_i^0$ the initial concentration of the basin.

The system (23), (24) has formally a family of ‘stationary’ solutions $b \in V$, $d_i^s$, $d_i \in L^\infty(\Omega)$, $i = 1, \dots, L$ such that $\psi(b) = \phi + \alpha$ with $\alpha$ an arbitrary continuous function depending only on time, and $\phi, c_i^s, d_i$ are solutions of the system of PDEs on $\Omega$: given $\psi, \varphi, c_i^0, r_i, \mu_i^e$, $i = 1, \dots, L$, find $\phi, d_i, d_i^s, i = 1, \dots, L$ such that for all $i = 1, \dots, L$

$$ \begin{array}{ll} -\mathrm{div}(d_i^s r_i \nabla \varphi) = d_i \frac{h^d}{T} & \text{on } \Omega, \\ \displaystyle\sum_{j=1}^L d_j^s = 1 & \text{on } \Omega, \\ \mu_i^s := \frac{r_i d_i^s}{\sum_{j=1}^L r_j d_j^s} = \mu_{i,e} & \text{on } \Gamma_e, \\ \displaystyle\sum_{j=1}^L d_j^s r_j \nabla \varphi \cdot \mathbf{n} = \varphi & \text{on } \partial\Omega, \\ d_i = \begin{cases} d_i^s & \text{if } h^d > 0 \\ c_i^0 & \text{if } h^d \le 0 \end{cases} & \text{on } \Omega. \end{array} \quad (26) $$
As in the single lithology case, only the bathymetry $b_T$ at time $t = T$ is retained which leads to the reduced model. Find $(b_T, a_T, d_i^s, d_i)$ solution of the following system of PDEs, for all $i = 1, \dots, L$:

$$ -\operatorname{div}(d_i^s r_i \nabla \phi) = d_i \frac{h^d}{T} \quad \text{on } \Omega, $$

$$ \sum_{j=1}^{L} d_j^s = 1 \quad \text{on } \Omega, $$

$$ \mu_i^s := \frac{r_i d_i^s}{\sum_{j=1}^{L} r_j d_j^s} = \mu_{i,e} \quad \text{on } \Gamma_e, $$

$$ \sum_{j=1}^{L} d_j^s r_j \nabla \phi \cdot \mathbf{n} = \varphi \quad \text{on } \partial\Omega, \qquad (27) $$

$$ d_i = \begin{cases} d_i^s & \text{if } h^d > 0 \\ d_i^0 & \text{if } h^d \le 0 \end{cases} \quad \text{on } \Omega, $$

$$ b_T = \psi^{-1}(\phi) \quad \text{on } \Omega, $$

$$ a_T = b_T - h^d \quad \text{on } \Omega, $$

$$ \frac{1}{|\omega_w|} \int_{\omega_w} b_T \, dx = b_w. $$

One can show in the one-dimensional case $d = 1$ that this model has a unique solution $(b_T, d_i)$ in $W^{1,\infty}(\Omega) \times L^\infty(\Omega)$ (see [GD05]).

Let $(b, c_i^s, c_i)$, $h = a - b$, be the solution of the multi-lithology stratigraphic model. The linearization of the reduced model for the preconditioning of the nonlinear Krylov algorithm is done around the point $d_i^s = c_i^s|_{t=T}$, $b_T = b|_{t=T}$, and

$$ d_i^0 = \frac{1}{h(., T)} \int_0^{h(., T)} c_i^0(., \zeta, T) \, d\zeta \quad \text{if } h(., T) \le 0. $$

## 6. Numerical experiments

The Newton–Krylov algorithms for the solution of (14) are expected to achieve good scalability properties with respect to the diffusion coefficients and the time and space discretization parameters. To evaluate numerically these properties, we shall consider a 3D basin test case using the multi-lithology model described in section 5. In this example, the sediments are a mixture of two lithologies ($L = 2$) characterized by their diffusion coefficients $r_1 K(b)$ and $r_2 K(b)$ with $r_2 = 1$, and $r_1 \ge 1$ is the ratio between the diffusion coefficients of the first and second lithologies. The non-negative function $K(b)$ is given by

$$ K(b) = \begin{cases} K_m & \text{if } b > 0, \\ r_{cm} K_m & \text{if } b \le 0, \end{cases} $$

where $K_m > 0$ is the marine diffusion coefficient of the second lithology, and $r_{cm} \ge 1$ is the ratio between the continental and marine diffusion coefficients. Basically, the higher the ratio $r_1$ the higher the coupling between the lithologies 1 and 2. Note also that for $r_1 = 1$ the model reduces to a single lithology model for the unknown bathymetry, and $r_1 = r_{cm} = 1$ corresponds to the linear parabolic model for the unknown bathymetry.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Test case: the initial topography $-b_0$ with input boundary fluxes $\varphi$, the sediment thickness $h^d$, and the exact solution $-b|_{t=T}$ of the inverse problem (14) for $r_{cm} = 10, r_1 = 10, n = 20$.</figcaption></figure>

The models are discretized by a Euler implicit time integration and a finite volume method with cell-centred variables as described in [EGGMT04] for the multi-lithology stratigraphic model (23), (24), and in [DMAC04] for the reduced model (27). For the sake of simplicity, the rectangular domain $\Omega = (0, 2l) \times (0, l)$, $l = 50$ km, is endowed with a regular Cartesian mesh of size $l/n$, and the time discretization is uniform of step $T/n_t$ with $T = 10$ My. The following tests investigate the scalability of the Newton–Krylov algorithms with respect to the parameters $r_{cm}, r_1, n$ and $n_t$ for a geometric mean marine diffusion coefficient $(r_1)^{1/2}K_m$ fixed to $10^{-2.5} \text{ km}^2 \text{ ky}^{-1}$.

For each value of the parameters, the sediment thickness $h^d = h|_{t=T}$ is computed using the discretized multi-lithology stratigraphic model (23), (24) with $b_0 = 50(1 - \cos(\pi(x+y)))$ m, $a_T = 400(1 - \cos(\pi(x+y)))$ m, in the definition of the accommodation $a = b_0 + \frac{l}{T}(a_T - b_0)$. The boundary fluxes vanish on $y = 0$ and $x = 2l$, and input fluxes are fixed to $\varphi = -0.75 \text{ km}^2 \text{ My}^{-1}$ on $y = l$ and to $\varphi = -1.5 \text{ km}^2 \text{ My}^{-1}$ on $x = 0$ with input fractional flow $\mu_{1,e} = \mu_{2,e} = \frac{1}{2}$.

The initial topography $-b_0$, the sediment thickness $h^d$, and the final topography $-b|_{t=T}$ are plotted in figure 3 for the parameter values $r_{cm} = 10, r_1 = 10, n = 20$.

The Newton–Krylov algorithms solve the inverse problem (14) of exact solution $a_T$ given the sediment thickness data $h^d = h|_{t=T}$, and the bathymetry $b_w$ actually fixed at a given cell located on the shoreline. In the following tests, the initial guess will always be given by the reduced model solution.

*Newton–Krylov algorithms parameters.* The stopping criterion $\|r(a_T^k) + J(a_T^k)p^k\| \le \epsilon_k \|r(a_T^k)\|$ of the GMRes iterative solver used for the computation of the inexact Newton direction $p^k$ at the $k$th nonlinear iteration is chosen with a decreasing sequence $\epsilon_k = \frac{1}{2^k}$ as discussed in [BS90]. In addition the maximum number of GMRes iterations is fixed to 50. The parameter $\delta$ of the step length selection procedure is set to $10^{-4}$.
<table><caption>**Table 1.** NME of the line search algorithm with $n = 20$, $n_t = 10$ and different values of the ratios $(r_1, r_{cm})$.</caption><thead><tr><th>(r<sub>1</sub>, r<sub>cm</sub>)</th><th>(10, 10)</th><th>(10, 100)</th><th>(100, 10)</th><th>(100, 100)</th></tr></thead><tbody><tr><td>No preconditioning</td><td>37</td><td>36</td><td>69</td><td>81</td></tr><tr><td>Preconditioning</td><td>17</td><td>23</td><td>28</td><td>35</td></tr></tbody></table>

The following parameters are used in the dogleg direction and trust region radius procedures of the trust region algorithm: $\bar{\rho}_1 = 0.25$, $\bar{\rho}_2 = 0.75$, $\beta_1 = 0.5$, $\beta_2 = 1.0$, $\beta_3 = 1.5$. The initial trust region radius is given by $\mu$ times the unconstrained Cauchy point length $\Delta_1 = \mu \frac{\|\nabla f(a_t^k)\|^3}{\|J(a_t^k)\| \|\nabla f(a_t^k)\|^2}$ for $\mu = 1.3$, and the maximum radius $\bar{\Delta}$ is given by $\bar{\Delta} = 10\Delta_1$.

The nonlinear algorithm is stopped after reduction of the maximum residual $\|r(a_T^k)\|_{\infty}$ by a factor $10^5$.

The complexity of the algorithms is measured in terms of the number of multi-lithology stratigraphic model evaluations. Let $c_l$ denote the average complexity of one linear system solution of the multi-lithology model Newton step, and let $N_{IL}$ denote the average number of Newton steps per time step in the multi-lithology model. The complexity of the model computation is $n_t N_{IL} c_l$ as well as the complexity of the Jacobian vector product $J(a_T^k) p$ computed by finite differences. It is assumed that the complexity of the gradient $\nabla f$ computation by adjoint state is roughly $n_t c_l$, and we neglect the complexity of the reduced model assuming that $n_t \gg 1$. Then, the line search algorithm cost in terms of an equivalent number of multi-lithology stratigraphic model evaluations is roughly equal to

$$NME = \frac{N_{\text{INL}}(n_t N_{\text{IL}} + n_t + N_{\text{GMRes}}n_t N_{\text{IL}}) + N_B n_t N_{\text{INL}}}{n_t N_{\text{IL}}}$$

and for the trust region algorithm, it is given by

$$NME = \frac{N_{\text{INL}}((N_{\text{IL}} + 1)n_t + N_{\text{GMRes}}n_t N_{\text{IL}}) + 2N_{\rho}n_t N_{\text{IL}}}{n_t N_{\text{IL}}}$$

where

* $N_{\text{INL}}$ denotes the number of Newton-Krylov iterations,

* $N_{\text{GMRes}}$ denotes the average number of GMRes iterations by Newton-Krylov step,

* $N_B$ denotes the number of additional cost function evaluations in the line search algorithm, and

* $N_{\rho}$ denotes the number of evaluations of $\rho$ in the trust region algorithm.

**Preconditioning.** In order to evaluate the efficiency of the preconditioning by the reduced model, we compare the line search algorithm with and without preconditioning for $n = 20$, $n_t = 10$ and different values of the ratios $(r_1, r_{cm})$. Table 1 shows that the preconditioning provides on these test cases an average gain of roughly a factor 2 compared with the algorithm without preconditioning.

**Robustness.** Tables 2-5 exhibit the NME of the line search algorithm for different values of the discretization parameters $n, n_t$, and of the diffusion coefficient ratios $(r_1, r_{cm})$. The dependence of the complexity of the algorithm on these parameters is clearly weak, which illustrates its good scalability properties.

**Globalization.** Tables 6 and 7 compare the preconditioned line search and trust region algorithms with $n = 20$, $n_t = 10$ and different values of the ratios $(r_1, r_{cm})$. For this basin
<table><caption>**Table 2.** NME of the preconditioned line search algorithm with $(r_1, r_{cm}) = (10, 10)$.</caption><thead><tr><th>$n_t$</th><th>10</th><th>20</th><th>40</th><th>80</th></tr></thead><tbody><tr><td>$n = 10$</td><td>13</td><td>18</td><td>18</td><td>18</td></tr><tr><td>$n = 20$</td><td>17</td><td>17</td><td>17</td><td>18</td></tr><tr><td>$n = 40$</td><td>18</td><td>38</td><td>38</td><td>27</td></tr></tbody></table>

<table><caption>**Table 3.** NME of the preconditioned line search algorithm with $(r_1, r_{cm}) = (10, 100)$.</caption><thead><tr><th>$n_t$</th><th>10</th><th>20</th><th>40</th><th>80</th></tr></thead><tbody><tr><td>$n = 10$</td><td>21</td><td>29</td><td>25</td><td>24</td></tr><tr><td>$n = 20$</td><td>23</td><td>24</td><td>25</td><td>24</td></tr><tr><td>$n = 40$</td><td>24</td><td>24</td><td>31</td><td>36</td></tr></tbody></table>

<table><caption>**Table 4.** NME of the preconditioned line search algorithm with $(r_1, r_{cm}) = (100, 10)$.</caption><thead><tr><th>$n_t$</th><th>10</th><th>20</th><th>40</th><th>80</th></tr></thead><tbody><tr><td>$n = 10$</td><td>25</td><td>38</td><td>40</td><td>30</td></tr><tr><td>$n = 20$</td><td>35</td><td>20</td><td>18</td><td>36</td></tr><tr><td>$n = 40$</td><td>20</td><td>96</td><td>25</td><td>42</td></tr></tbody></table>

<table><caption>**Table 5.** NME of the preconditioned line search algorithm with $(r_1, r_{cm}) = (100, 100)$.</caption><thead><tr><th>$n_t$</th><th>10</th><th>20</th><th>40</th><th>80</th></tr></thead><tbody><tr><td>$n = 10$</td><td>32</td><td>36</td><td>52</td><td>46</td></tr><tr><td>$n = 20$</td><td>35</td><td>43</td><td>39</td><td>33</td></tr><tr><td>$n = 40$</td><td>20</td><td>36</td><td>23</td><td>36</td></tr></tbody></table>

<table><caption>**Table 6.** NME for the preconditioned line search and trust region algorithm for $n = 20, n_t = 10$, and different values of the ratios $(r_1, r_{cm})$.</caption><thead><tr><th>$(r_1, r_{cm})$</th><th>(10, 10)</th><th>(10, 100)</th><th>(100, 10)</th><th>(100, 100)</th></tr></thead><tbody><tr><td>Line search</td><td>17</td><td>24</td><td>20</td><td>43</td></tr><tr><td>Trust region</td><td>32</td><td>35</td><td>60</td><td>46</td></tr></tbody></table>

<table><caption>**Table 7.** NME for the preconditioned line search and trust region algorithm for $n = 40, n_t = 20$, and different values of the ratios $(r_1, r_{cm})$.</caption><thead><tr><th>$(r_1, r_{cm})$</th><th>(10, 10)</th><th>(10, 100)</th><th>(100, 10)</th><th>(100, 100)</th></tr></thead><tbody><tr><td>Line search</td><td>38</td><td>24</td><td>96</td><td>36</td></tr><tr><td>Trust region</td><td>31</td><td>48</td><td>56</td><td>62</td></tr></tbody></table>

test case, we conclude that both algorithms exhibit a good scalability. As expected, the trust region algorithm seems more robust as seen for $(r_1, r_{cm}) = (100, 10)$ but it is also more costly for most of the other values of the ratios $(r_1, r_{cm})$.
# 7. Conclusion

Stratigraphic models at large geological space and time scales are used in oil exploration to obtain a quantitative description of the evolution of basin geometry and rock properties. Such stratigraphic models must integrate seismic and well-log data in order to be predictive.

In this paper, a new ‘stationary’ stratigraphic model has been introduced which is used as a reduced model to solve the inverse stratigraphic modelling problems. This reduced model can be used as a good approximation of the inverse problem solution (4) in order to identify the diffusion coefficients and boundary fluxes from well-log data (see [DMAC04, GD05]).

In order to obtain a more accurate estimation of the accommodation, the reduced model is also used as a preconditioner combined with Newton–Krylov approaches to fit the thicknesses of the sedimentary layers measured from the seismic data between each successive seismic marker. Numerical experiments show that the preconditioned Newton–Krylov algorithms with line search or trust region globalization have very good scalability properties with respect to the diffusion coefficients and the time and space discretization parameters.

# Acknowledgments

The authors would like to thank Jean Charles Gilbert from INRIA and Delphine Sinoquet from the Institut Français du Pétrole for fruitful discussions on the Newton–Krylov algorithms.

# References

[Au63] Aubin T 1963 *Un théorème de compactité* *Note C. R. Acad. Sci.* **256** 5042-4

[AH89] Anderson R S and Humphrey N F 1989 *Interaction of weathering and transport processes in the evolution of Arid landscapes* *Quantitative Dynamics Stratigraphy* ed T A Cross (Englewood Cliffs, NJ: Prentice-Hall) pp 349–61

[BS90] Brown P N and Saad Y 1990 *Hybrid Krylov methods for nonlinear systems of equations* *SIAM J. Sci. Stat. Comput.* **11** 450–81

[BS94] Brown P N and Saad Y 1994 *Convergence theory of nonlinear Newton Krylov algorithms* *SIAM J. Opt.* **4** 297–330

[DMAC04] Dobranszky G, Masson R, Albouy E and Cacas M C 2004 A 4D stratigraphic model fitting seismic data for sedimentary basin modelling *Proc. ECMOR IX, EAGE*

[GD05] Dobranszky G 2005 *Systèmes d'aide à l'inversion de modèles stratigraphiques* *PhD Thesis* Université de Nice

[EGGM03] Eymard R, Gallouët T, Gervais V and Masson R 2005 *Convergence of a numerical scheme for stratigraphic modeling* *SIAM J. Numer. Anal.* **43** 474–501

[EGH00] Eymard R, Gallouët T and Herbin R 2000 *The finite volume method* *Handbook of Numerical Analysis* vol 7 ed P G Ciarlet and J L Lions (Amsterdam: Elsevier) pp 715–1022

[EGGMT04] Eymard R, Gallouët T, Granjeon D, Masson R and Tran Q H 2004 *Multi-lithology stratigraphic models under maximum erosion rate constraint* *Int. J. Numer. Methods Eng.* **60** 527–48

[GM05] Gervais V and Masson R 2005 *Numerical simulation of a stratigraphic model at geological time scales* *Comput. Geosci.* submitted

[GM04] Gervais V and Masson R 2004 *Mathematical and numerical analysis of a stratigraphic model* *M2AN* **38** 585–612

[GD03] Granas A and Dugundji J 2003 *Fixed Point Theory* (New York: Springer)

[G97] Granjeon D 1997 *Modélisation stratigraphique déterministe conception et applications d'un modèle diffusif 3D multilithologique* *PhD Dissertation* Géosciences Rennes, Rennes, France, p 189

[GJD98] Granjeon D, Joseph P and Doligez B 1998 *Using a 3-D stratigraphic model to optimise reservoir description* *Hart's Petroleum Engineer International (Nov.)* pp 51–8

[GJ99] Granjeon D and Joseph P 1999 *Concepts and applications of a 3D multiple lithology, diffusive model in stratigraphic modelling* *Numerical Experiments in Stratigraphy* (SEPM Special Publication vol 62) ed J W Harbaugh et al (Tulsa, OK: SEPM)
[LSU68] Ladyzenskaja O, Solonnikov V and Ural'ceva N 1968 Linear and quasilinear equations of parabolic type *Transl. Math. Monogr.* **23** 257

[NW99] Nocedal J and Wright S J 1999 *Numerical Optimization (Springer Series in Operations Research)* (Berlin: Springer) p 293

[R92] Rivenaes Jan C 1992 Application of a dual lithology, depth-dependent diffusion equation in stratigraphic simulation *Basin Res.* **4** 133–46

[TH89] Teztlaff D M and Harbaugh J W 1989 *Simulating Clastic Sedimentation* (New York: Van Nostrand-Reinhold)

[TS94] Tucker G E and Slingerland R L 1994 Erosional dynamics, flexural isostasy, and long-lived escarpments: a numerical modeling study *J. Geophys. Res.* **99** 12229–43