Identification of a stratigraphic model with seismic constraints

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.


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 [ 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 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 R d independent of time, denoted by , and defining the horizontal extension of the basin. We shall denote by x ∈ 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 ∈ , 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 ∈ and time t. The second input parameter of the model is the boundary flux ϕ defined at the boundary ∂ of 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 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, . . . , L, and the model accounts for the conservation equations of each lithology i = 1, . . . , 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 ∈ , and the sediment thickness h(x, T ), x ∈ . 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 .
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 ϕ 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.

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 −∇(h + H ) or equivalent to ∇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 ϕ at the boundary ∂ . The input flux ϕ < 0 represents the sediment supply at river entries. The output flux ϕ 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, ϕ, b 0 , K(b), the stratigraphic model finds the solution h such that (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 0) than in marine (b > 0) environments. It is typically given by with β 0, and K c K m > 0. The model (1) can be viewed as a nonlinear parabolic equation for the bathymetry b with Neumann boundary conditions ϕ, and the rate of accommodation change ∂ 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.

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, 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 ∂ t a = 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: Due to the Neumann boundary condition ϕ set on the stratigraphic model, it must be assumed that the condition on the data always holds in the following. This condition is necessary for the existence of a solution to (4). This also explains the equation 1 (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 ω w ⊂ representing the given well sections of measure |ω w |.
Remark 2.1. In practice, the remaining sediment supply ϕ 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 ϕ. This is to be viewed as a constraint on the sediment supply ϕ that will in practice always be satisfied from the identification procedure of ϕ 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.

Linear model
Let us consider the inverse problem (4) in the linear case, for which the diffusion coefficient K is a given function in endowed with its usual norm denoted by V . One can prove the following proposition.

the inverse problem (4). This solution satisfies the bound
with C a constant independent of h d , b 0 , ϕ and b w .
Proof. Let us set and The inverse problem (4) equivalently amounts to finding and One can show that the solutionb| t=T of (8) for v given, defines a linear bounded operator L from L 2 0 ( ) to L 2 0 ( ) which is a strict contraction in L 2 0 ( ). It results that equation (9) Assuming furthermore that K is a constant, the solution of (4) is explicit in the spectral basis (v i ) i∈N as given by the following proposition.

Proposition 3.2. If K is a constant, the solution of proposition 3.1 is defined by
with C a constant independent of any parameter.

the inverse problem (4). This solution satisfies the bound
with C a constant independent of h d , b 0 , ϕ and b w .

Preconditioned Newton-Krylov algorithms for the inverse problem
Let us define H : L 2 ( ) → L 2 ( ) such that H(a T ) is the unique solution h| t=T = a T − b| t=T in L 2 ( ) of the stratigraphic model (1)-(3) for a given a T ∈ L 2 ( ) (see [LSU68]). Our objective is to build an algorithm for the solution a T ∈ L 2 ( ) of the nonlinear problem 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 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.

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 we obtain a family of solutions b such that In (15), α is an arbitrary continuous function depending only on time such that The accommodation given by a = ψ −1 (φ + α) + 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.
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 0 T of the accommodation a T . Then, for k 0, given a k T , the stratigraphic model (1) computes the sediment thickness H a k T , the bathymetry b k and the residual r a k T = H a k T − h d . Let C(b) denote the linear bounded operator from L 2 0 ( ) defined by (7) to such that C(b)r is the solution p of the following linearized reduced model: The quasi-Newton correction step is defined by a k+1 T = a k T + C(b k )r a k T , and the algorithm is iterated until a stopping criterion for the L 2 norm of the residual r a k T 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 ρ = 1 − 1 + 1 Kω 1 T (1 − α(Kω 1 T )) > 0 bounded by 1 2 independently of the diffusion coefficient K. By contrast, a simple fixed point algorithm a k+1 T = a k T − r a k T + 1 |ω w | ω w r a k T dx has a convergence rate ρ = 1−e −Kω 1 T Kω 1 T which deteriorates rapidly for small diffusion coefficients.
Nevertheless, in the nonlinear case, C(b k )r a k T 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 C(b).

Preconditioned Newton-Krylov algorithms
Let (b, a T ) be the current solution of the stratigraphic model (1)-(3). The Jacobian J (a T ) of H(a T ) is the linear operator from L 2 w ( ) to L 2 0 ( ) such that for any q ∈ L 2 w ( )J (a T )q = q − δb| t=T where δb is the unique solution of the parabolic equation The Newton step p ∈ L 2 w ( ) is the solution of the equation

J (a T )p = −r(a T ),
where r(a T ) = H(a T ) − h d is the residual. This equation is solved using a GMRes linear iterative solver with the right preconditioner C(b). The matrix vector product is computed by the finite difference formula for a suitable choice of 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 ) = 1 2 r(a T ) 2 L 2 ( ) defined in L 2 w ( ). Its gradient ∇f is the function of L 2 w ( ) defined by

(J (a T )v, r(a T )) L
i.e. the orthogonal projection on L 2 w ( ) of the gradient of the merit function f when defined on the full space L 2 ( ). This latter gradient is computed using the adjoint state of the stratigraphic model.

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 α k to move from the current iterate a k T to a new iterate a k+1 T with a lower merit function value. The new iterate is given by a k+1 T = a k T + α k p k . In our case, p k is an inexact Newton step which approximates the solution of the Newton equation J a k T p k = −r a k T using a preconditioned GMRes algorithm. The step length α k is an approximate minimizer of the one-dimensional function φ defined by φ(α) = f a k T + αp k , α > 0, which must satisfy the Armijo condition [NW99] f a k T + α k f a k T + δα k ∇f a k T , p k L 2 ( ) with 0 < δ < 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 φ using φ(0), φ(1) and φ (0). Note that the Armijo condition can be rewritten as follows: φ(α k ) φ(0) + δα k φ (0). This procedure can be seen as an enhancement of the backtracking approach [NW99].

Procedure
Step length selection the merit function f a k T 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 k T = 1 2 r a k T 2 2 , the quadratic model m k is defined by and the step p k is an approximate solution of the sub-problem where k is the radius of the trust region. The exact solution of (20) is denoted by p * ( k ).
Dogleg direction. The dogleg method (see [NW99]) finds an approximate solution to (20) by replacing the curved trajectory p * ( k ) with a path consisting of two line segments. The first segment runs from the origin to the Cauchy point p k T ) 2 ∇f a k T , and the second segment runs from the Cauchy point to the inexact Newton direction denoted by p k N ∼ −J a k T −1 r a k T . 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 k T and the quadratic model m k at the previous iteration, which is evaluated by the ratio ρ k 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ρ 1 andρ 2 such thatρ 1 <ρ 2 and three multiplicative factors β 1 , β 2 and β 3 such that β 1 < β 2 < β 3 are introduced. The new radius is computed using the trust region radius update procedure.

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 where c 0 i , 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 × (0, T ) accounting for the conservation of the lithology fraction 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 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.

Adjoint state
The computation of the gradient ∇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 s i , 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 × (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, . . . , L, only requiring the 3D storage of c h i = c i | ζ =h , c T i = 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).

Reduced model
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 s i , d i solution of the following system of PDEs, for all i = 1, . . . , L: One can show in the one-dimensional case d = 1 that this model has a unique solution

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 1 is the ratio between the diffusion coefficients of the first and second lithologies. The non-negative function K(b) is given by where K m > 0 is the marine diffusion coefficient of the second lithology, and r cm 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. 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 = (0, 2l) × (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 km 2 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(π(x + y))) m, a T = 400(1 − cos(π(x + y))) m, in the definition of the accommodation a = b 0 + t T (a T − b 0 ). The boundary fluxes vanish on y = 0 and x = 2l, and input fluxes are fixed to ϕ = −0.75 km 2 My −1 on y = l and to ϕ = −1.5 km 2 My −1 on x = 0 with input fractional flow µ 1,e = µ 2,e = 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 k T + J a k T p k k r a k T of the GMRes iterative solver used for the computation of the inexact Newton direction p k at the kth nonlinear iteration is chosen with a decreasing sequence k = 1 2 k as discussed in [BS90]. In addition the maximum number of GMRes iterations is fixed to 50. The parameter δ of the step length selection procedure is set to 10 −4 . The following parameters are used in the dogleg direction and trust region radius procedures of the trust region algorithm:ρ 1 = 0.25,ρ 2 = 0.75, β 1 = 0.5, β 2 = 1.0, β 3 = 1.5. The initial trust region radius is given by µ times the unconstrained Cauchy point length T ) 2 for µ = 1.3, and the maximum radius¯ is given by¯ = 10 1 . The nonlinear algorithm is stopped after reduction of the maximum residual r a k T ∞ 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 )p computed by finite differences. It is assumed that the complexity of the gradient ∇f computation by adjoint state is roughly n t c l , and we neglect the complexity of the reduced model assuming that n t 1. Then, the line search algorithm cost in terms of an equivalent number of multi-lithology stratigraphic model evaluations is roughly equal to and for the trust region algorithm, it is given by 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 2. NME of the preconditioned line search algorithm with (r 1 , r cm ) = (10, 10).  Table 3. NME of the preconditioned line search algorithm with (r 1 , r cm ) = (10, 100).  Table 4. NME of the preconditioned line search algorithm with (r 1 , r cm ) = (100, 10).  Table 5. NME of the preconditioned line search algorithm with (r 1 , r cm ) = (100, 100).  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 ).

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.