A GEOLOGICAL DELAYED RESPONSE MODEL FOR STRATIGRAPHIC RECONSTRUCTIONS

. We are interested by a nonlinear single lithology diﬀusion model adapted from ideas originally developed by the Institut Fran¸cais du P´etrole (IFP). The geological stratigraphic modeling has to describe transports of sediments, erosion and sedimentation processes by taking into account a limited weathering condition; the method by which the history of a sedimentary basin is revealed relies on knowledge of both initial and ﬁnal data and can be generalized to multiple lithology. For this purpose, we introduce a relaxation time related to a delayed response for establishing equilibrium states; this approach introduces regularizing eﬀects according to the ideas of G.I. Barenblatt - S. Sobolev and J.-L. Lions - O. A. Oleinik. New well-posedness results are presented.

1. Introduction and geophysical motivations. The model addressed in this work inters directly in the framework of theoretical stratigraphic researches developed in the last decade at the French Petroleum Institute ( Institut Français du Pétrole, IFP), (cf. D. Granjeon et al. [18], [19], [20], [21], [26], [27]). They concern the evolution in a sedimentary basin of multilithological strata under the constraint of a maximal erosion rate depending on the composition of the mixture and of the climatic environment. Bathymetric and topographic lows are loaded by the influx of sediment, which leads to the accumulation of sediment depending on the steepness of the hillslopes. These models do not enter in a known mathematical framework and numerous scientific issues are still open. The needs to understand the processes of basin formation and evolution are obviously of great importance to aid evaluation of potential hydrocarbon reserves ( [13], [28], [29]). Indeed, sedimentary basins are the location for almost all of the world's oilfields and as such are increasingly the focus of intense academic and industrial interest. As a key among others to understand sedimentary basins as geodynamic entities, the mathematical nonlinear analysis can bring out through models for sediment dispersal and basin fill how the basin has been formed at large time and space scales, how the sediment was transported or deposited, and reveal sources of the sediment filling.
A stratigraphic model reports on the transport of erodible sediments in a sedimentary basin when the tectonic, the epirogeneis (vertical movement of a large part of the earth's crust), the eustatism (slow variation of the ocean level) and the sediment fluxes at the boundaries of the basin are known. That introduces different diffusive physical constants depending on the position considered (in marine or earth's environment). The sediment transport only takes place on surface: the fluxes interring through a part of the boundary are known whereas the outflow fluxes on the complementary part of the boundary are determined from an unilateral constraint when specifying the maximal outflow flux. The original modeling proposed couples in a compatible way a process of gravity flow for the sedimentation and erosion on one part, and a feedback control system regulated by the limitation of the erosion on the other part (the erosion velocity is controlled by an unilateral constraint). The coupling of the two models (the mass conservation equation and the constraint on the time variation) is based on a variable factor λ (flow limiter flux to be chosen "at best" according to a maximum criterion evolving in the interval [0, 1]) and on the unilateral conditions linking the double constraint of maximum erosion rate and maximum value of λ.
The general problem of petrography focuses on a sedimentation and erosion model in sedimentary basins with several lithologies of constant porosity (compaction phenomena in the process of lithification are not considered to seriate the difficulties; see [24] for specific developments)). Its main characteristics are as follows: -the material flow is proportional to the gradient of the height of sediments in the linear gravity flow model. This model may be generalized by introducing a quasi-linear diffusion term. Otherwise, the effects of a wind or sub-marine transport of sediments, possibly periodic in the case of tides, may be taken into account.
-a maximum erosion velocity is taken into account ("weathering-limited system" model).
-an instantaneous feedback control system determine the variations of the erosion inside the domain and possibly on a fixed part of the boundary.
The originality of this approach relies on the fact that the ideal static model of instantaneous equilibrium is not a Cauchy-Kowaleska system [23]. The mathematical associated model, ill-posed, presents some real difficulties of analysis. Consequently, we take into account a relaxation time modeling a short waiting period for the realization of the equilibrium, in response to an instantaneous perturbation. This leads to a double regularizing effect, through the introduction of corrective kinetic term of Barenblatt-Sobolev pseudo-parabolic type ( [2], [3], [4], [7] associated with an elliptic regularization in the sense of J.-L. Lions and O. Oleinik ([30], [31]). The formulation then becomes well-posed to reconstruct for long periods of geological time the evolution of the basin on a time internal [0, T ] from the knowledge of the assumed initial topography, of the erosion-sedimentation rate observed at the final time T and of the boundary conditions as well. In addition, these models are advantageously really active in the sense where the initial state never constitute one or "the" solution, that would be unrealistic for geologists. Moreover, some calculations let understand the effects of the vanishing regularizations to formulate in a relevant way the ill-posed stiff diffusion model.

2.
Presentation of the weathering-limited model.

2.1.
Modeling of the geologically relevant features. The horizontal projection of the basin extension is modeled by the bounded domain Ω ⊂ R d , d = 1, 2 with a Lipschitz 1 continuous boundary ∂Ω; the bedrock position is determined by a known vertical position H = H (t, x) at time t and position x, (t, x) ∈ Q = ]0, T [ × Ω, where T is both the observation time and the age of the basin. One denotes by u the thickness of the sediments and the topography is given by H + u. One considers a gravitational model where the main basin forming mechanism is sediment loading: i) the sediments flux is assumed to be proportional to κ∇Φ (H + u) where κ is a viscosity rate 2 (κ = κ (x, u) optionally [35]) and Φ a non-decreasing function.
ii) the sedimentation-erosion velocity is governed by a weathering limited process; hence the partial derivative ∂u ∂t is underestimated by a limit erosion velocity E = E (t, x), depending on the climate and the age of the sediments; significantly, the erosion velocity increases with decreasing grainsize (according to the Hjulström-Sundborg diagram [6] which shows several key concepts about empirical relationships between erosion, transportation and deposition; see also [37]). One is led to the inequality with a moving obstacle ∂u ∂t iii) in order to make compatible this unilateral constraint and a necessary conservative formulation, D. Granjeon et al. [26] proposed to correct the diffusive flux by introducing a dimensionless flux limiter 3 λ = λ (t, x), an unknown function with values in [0, 1]. With the aim of providing a mathematical modeling for the limiter λ, R. Eymard et al. ( [18], [19], [21]) propose to consider that the unknowns λ and u are limited by the following non-standard constraints: associated to the mass balance of the sediments, when the density remains constant, ∂u ∂t − div (κλ∇Φ (H + u)) = 0 a.e. in Q It means that if the erosion rate constraint is inactive, precisely in the set {(t, x) ∈ Q, λ = 1}, the flux is equal to the diffusive one; otherwise, the flux has to be corrected in order to take into account the realistic maximum erosion velocity and the constraints lead to the relation ∂u ∂t Obviously, the constraints (2) and the associated free boundary problem constitute the main characteristic difficulty for the theoretical analysis. The boundary conditions may be of different miscellaneous types at fixed regions Γ i of ∂Ω as, by way of generic illustration ( [2], [23], [26], [38]): κλ∇Φ (H + u) .n + f = 0 a.e. on ]0, T [ × Γ 3 . (the data g, E and f are fixed), completed by specifying the initial condition u (0, .) = u 0 a.e. in Ω.

2.2.
Heuristic formulation of the mathematical problem: A partial differential inequality. We present in this section the main aspect of the model so as not to cloud the presentation with secondary aspects, even the didactic simplifications proposed may surprise geologists. This model problem is correctable using easy mathematical generalizations. The unknowns are the topography h = H + u (u the thickness of the sediments) and λ the erosion limiter. We consider in the sequel that the profile H of the bedrock is stationary 4 , the maximum erosion rate 5 E is a non negative constant; we reduce the academic part of the study to homogeneous Dirichlet conditions on the boundary 6 . In fact, we have boundary problems with incomplete data, as in the works of A. Cimetière et al ( [9], [10], [11]). With a regularization term, the problem consists in recovering missing data (displacements and forces) on some parts of the boundary of the domain from the knowledge of overspecified data on the remaining parts.
Then equations 7 and inequalities (1), (2) lead to By introducing β the Heaviside maximal monotone graph (i.e. β (r) = 0 if r < 0, β (0) = [0, 1] , β (r) = 1 if r > 0), we derive the following extended formulation for solutions of constrained problem (4) via the non-standard differential inclusion: We search an element λ and a function h which satisfy in addition to the boundary conditions: It is well known that the problem of ideally instantaneous equilibrium is ill-posed ( [2], [23]). For an example, in ( [2], [38]), the model of sedimentation without erosion (i.e. when E = 0) presents both degenerate hyperbolic and parabolic behaviors (barrier effect, finite speed behavior and dead-zones) because of the terms is obvious that a selection criterion must be added when several choices of limiters are possible [23]. In addition, the mathematical treatment faces to the subtle difficulty 8 : 4 The source term ∂H ∂t is to be taken into account in case of epirogenic movement. 5 In geological practice, the diameter of Ω can be about several hundreds of kilometers, T = 1 to 10 megayears, E = 5 m/megayear [18]. E = 0 in sedimentation modeling without erosion. 6 For a suitable choice by vertical displacement of the zero altitude. 7 Let us note that a similar equation describes one dimensional vertical unsteady flow in a non-homogeneous multi-stratum soil ( [32], [34]). 8 Evidence provided by Profs François Bouchut, Univ. Paris-Est-Marne La Vallée, France (personal communication) and Luigi Ambrosio, Scuola Normale Superiore in Pisa, Italy [1].
in the sense of distributions, it is not possible to say, without complementary information, that div F = 0 a.e. in {F = 0} [1]. Under these assumptions, there exist some situations where div F = 0 in {F = 0}. This general but very disconcerting remark on equations of divergential type concerns here the case F = λ∇h. On the set with free boundary ∂h ∂t + E < 0 , a priori non negligible, where the flow limiter λ is zero, it is not possible to conclude that div (λ∇h) = 0. To our knowledge, it is out of range to prove that the condition ∂h ∂t + E ≥ 0 is implicitly satisfied in the formulation (6). Moreover, the equivalence of the formulations (4) and (6) is not clearly stated. If d = 1, this difficulty cannot happen for distributional derivatives (Saks lemma).

2.3.
Regularizations for an ill-posed diffusion equation. The initial formulation must be now enriched to obtain a well-posed problem. The general idea relies on the introduction of small parameters (relaxation time, viscosity coefficient, etc) which lead to well-posed systems. Passing to the limit when these parameters tend to zero, with enough a priori estimations, we expect to obtain a limit that may be considered as a definition of the inital ill-posed problem and, in the same time, a supplementary discriminant property (entropy law in a broad sense [17], [22], chap. 4, principle of maximum, etc).
Other approaches are possible, especially in the paper by R. Eymard and T. Gallouët [21] thanks to the use of numerical schemes; the limiter on the fluxes happens to satisfy a stationary scalar hyperbolic inequality within a double obstacle problem; the solution is shown to be the maximal element of an appropriate convex set of functions. Indeed (see also [23], p. 717), for a given t 0 ∈ ]0, T [, the limiter λ, considered as an isolated unknown, verifies for a.e.
In this study, we introduce a relaxation time τ strictly positive, necessary to establish the equilibrium. The conservation law is then considered on the form: The assumed littleness of τ with respect to T enables the kinetic correction at the first order: Since, according to the usual chain rule, valid for the weak derivatives of Sobolev functions and taking into account (4), we get: Expression (7) may be written as ("viscous" regularization for τ > 0) Equivalently, introducing a generalized Green operator related to an homogeneous Dirichlet condition, a form of the resolvent, we get: This pseudo-parabolic Sobolev-Barenblatt regularization constitutes a generalization of the Yosida approximation method [16]. In general, a second order term ∆Φ (v) is replaced by the third order term ∆Φ (v µ ) + µ∆ ∂v µ ∂t . In the light of P.I.
Plotnikov works, developed later by L.C. Evans et M. Portilheiro [16], [17], this procedure results in a remarkable complementary information of entropy-type in the broad sense, which enables to control hysteresis phenomena in phase transition models. In thermodynamics, the term w µ = Φ (v µ ) + µ ∂v µ ∂t represents the chemical potential whose decreasing governs the sense of migration of chemical species 9 . Following this idea of relaxation, a translation with respect to the time t leads to the conservation law and with a development of the velocity at the first order, we obtain a.e. in Ω with the consistency condition z T (.) + E ≥ 0 a.e. in Ω. In this formal approach, all the variables play the same role. Using a consistent change of variables, both regularizations can be combined in the following formulation. Introducing two strictly positive parameters η and τ (time dimensioned): Equivalently, introducing a backward operator for heat equation with homogeneous Dirichlet boundary conditions whose regularizing effects are well known, we obtain: 3. Definition of solution and existence results. To our knowledge, existence and uniqueness and even well-posedness results for a solution to the above differential inclusions are still widely open problems. Our purpose is to analyze a modified model where the graph β is replaced by a function b, an approximation in a suitable sense of β. Hereinafter, very generally, b is assumed to be a continuous function and whose modulus of continuity is denoted by ω ; various moduli are used as ω (r) = ω 0 r α , ω 0 > 0, 0 < α < 1 or α ≥ 1 2 (Hölder-continuous functions whose square of ω is a Osgood function 10 ) and then Let us introduce the spaces V = H 1 0 (Ω) and H = L 2 (Ω), which are Hilbert spaces with the usual scalar products 11 , the embedding of V in H being dense and compact. Classically, we identify the pivot space H with its topologic dual, such as H −1 (Ω), the dual of V , can be identified to an over-space of H. Moreover, we note < ., . > the duality bracket between V and V , that coincides with the scalar product of L 2 (Ω) for elements of H × V.
In addition, for time dependent problems in the interval [0, T ], T > 0, , we introduce the space which is a separable Hilbertian space endowed with its natural topology. Therefore, we are interested in the η−elliptic, τ −pseudo-parabolic regularized problem and we derive the following sense for a solution: Or, equivalently, integrating by parts, in a global formulation, 10 i.e.
Proof. The result is obtained by using non-constructive Schauder-Tichonoff fixed point theorem in the framework of separable Hilbertian spaces endowed with the weak topology ( [22], p. 50 for example). For any ϕ ∈ W (0, T ), we consider as resulting from a classical result [30] the unique solution h = h (ϕ) ∈ W (0, T ) to the auxiliary linearized problem (P ϕ,η,τ ): in Ω.
Since v = h − h 0 and v = ∂h ∂t are available test-functions, one gets that and another law for energy conservation So, from the consideration of the application H : ϕ ∈ W (0, T ) → h (ϕ) ∈ W (0, T ) and from easy calculations to get a priori estimates, we can prove that exist constants C 1, C 2, C 3 independent of ϕ such that the convex closed bounded nonempty K is stable by the application H . The set K is weakly compact in W (0, T ). We prove easily that the application H is weakly-weakly sequentially continuous from K to K with the help of the compacity of the embedding from (Dubinskii's theorem, [31], p. 57 and 141). It follows that the application H admits a fixed point which is by construction solution of the problem (P η,τ ).
is a test-function and it comes that 4. The singular case of the sedimentation process.

4.1.
Illustration of degeneration effects in the initial model. As it has already been highlighted in [2], p. 215-216, the model of sedimentation without erosion (4), (5) presents some singularities not encountered in models with erosion. In this case, we have E = 0 and the constraint imposes that ∂h ∂t ≥ 0 a.e. in Q. For an example, the existence of dead zones and barrier effects are highlighted.
Consider v = h + as a test function in (4 So, for any t, h + (t, .) = h + 0 a.e. in Ω. In other words, the set {h 0 ≥ 0} is a dead-zone. Furthermore, the relation λ ∂h ∂t ∇h + = 0 a.e. in Q shows that the limiter can take the value zero and the equation describing the mass balance of the sediment is then totally degenerate. Otherwise, assume that there exists a compact set K and a regular open set O with K ⊂ O ⊂ Ω and O \ K ⊂ {h 0 ≥ 0}. Then, thanks to Urysohn's lemma, a test function v such that 1 K ≤ v ≤ 1 O exists and from (4) we prove that, for any t in ]0, T [, h (t, .) = h 0 a.e. in O. This expresses precisely that any zone surrounded by a zone where h 0 ≥ 0 is stationary, out of reach by the phenomenon (a barrier effect).
Here is the key to understanding totally degeneration effects on the system and surprising behaviors specifically when E = 0: a changing type equation illustrated by regions of finite, null or infinite velocity propagation, delimited by free boundaries.
For the same reasons related to the degenerated hyperbolic type of the constrained equation, one is able to construct traveling-wave solutions and bring out properties of finite velocity of propagation (see [2] for details).

4.2.
Smoothing regularization of sedimentation models. The regularized model introduces a continuous function b in order to have an approximation (in a sense to be specified) of the maximal monotone graph of the Heaviside function. In a classical way, we introduce the Yosida approximation b ε such as b ε (r) = min r + ε , 1 for any parameter ε > 0.
The practical interest of this choice relies on the fact that the function b ε is a Lipschitz-function. Therefore, for any function v ∈ H 1 0 (Ω), we have b ε (v) ∈ H 1 0 (Ω) with calculation rules close to classical calculation for expressions of ∇b ε (v) in the sense of distributions.
However, the numerical tests are unsatisfactory for such as choice. In addition, the theoretical model is sometimes "lazy" in the sense that the only function obtained is the initial state without any activity, that is to say a parasitic function. Accordingly, we introduce here in an original way, another kind of approximation by imposing continuous functions vanishing in 0, whose graph has a vertical tangent in zero.
As an example, we way consider functions such as b ε (r) = min r + ε α , 1 for any parameter ε > 0 and α fixed, 0 < α < 1. Generally, we consider the class of continuous functions such as Such as function is no more a Lipschitz-function. So that we introduce B the primitive of the function 1 b , defined on R + and vanishing at 0. In the next, we denote: for r ≥ 0 and δ > 0. (11) Proposition 2. (a priori estimate, energy conservation law). If h is any solution to the problem of the sedimentation case in the sense of definition 1 and under hypotheses (10), h satisfies the a priori estimate and energy conservation law More precisely, we have the following energy conservation law involving the initial and final data: Proof. Considering the test function v = B δ ∂h ∂t in (9) leads to and in particular: Thanks to both the Beppo Levi's and Lebesgue's dominated convergence theorems, passing to the limits for δ tending to zero leads to

GÉRARD GAGNEUX AND OLIVIER MILLET
or equivalently To this end, we dispose of an uniform upper bounding, independent from δ > 0, Indeed, these two quantities vanish on Ξ = (t, x) ∈ Q, ∂h ∂t = 0 as b ∂h ∂t and ∇ ∂h ∂t vanish on Ξ. The latter results from the fact that ∂h ∂t ∈ L 2 [0, T ] ; H 1 0 (Ω) ; the property of pointwise convergence is trivial on the complement. The energy conservation law results from the limiting procedure applied to (14) for δ tending to zero, thanks to the Beppo Levi's and Lebesgue's dominated convergence theorems, without using lower-bounds.

4.3.
Effects of vanishing relaxation times. The important point that must be analyzed now concerns the limit behavior of the solutions when the small parameters η > 0 and τ > 0, responsible for the regularizing effects, tends to zero. Let us understand the effects of the regularization for the ill-posed initial diffusion problem. We denote h = h η,τ the corresponding solutions to underline the dependence on these short relaxation delays. A response to this question is contained in the following proposition which results immediately from (12) and which could serve as a definition of ultra-weak solution to the original problem: Proposition 3 (behaviors at vanishing relaxation times). From the sequence {h η,τ }, for η and τ tending to zero, we can extract a subsequence, denoted similarly, and there exists a limit h ∈ H 1 [0, T ] ; L 2 (Ω) such as in Ω, h η,τ (T ) converges weakly to h (T ) in H 1 0 (Ω). Moreover, by using the lower semicontinuity property of the norm for the weak corresponding topology, passing to the limits in (15) for the extracted sequence {h η,τ } yields expected property revealing a less steep final landscape than originally, and the energy estimate More precisely, according to (13) and the convexity of the function r → rB (r) , r ≥ 0, we have: Actually, T can be replaced by any t > 0 in the three previous estimates and h 5. On uniqueness questions in stratigraphy. In this section we consider the general relaxed problem (P η,τ ) with the constraint on the time variation ∂h ∂t The uniqueness of the solution is still a widely open question. To our knowledge, the existing results are based on very technical methods to control the iterative terms concerning expression Ω b( ∂h ∂t + E)∇h.∇vdx. As an illustration, to limit the calculations, we will consider in the next the case where the relaxation parameter τ is strictly positive whereas η is equal to zero (i.e. the problem (P 0,τ )). The given final condition is not necessary anymore and the initial topography will be supposed to have no discontinuity. Accounting for the dimension of the space (d = 2 in practice) and of the associated Sobolev injections 14 , the following supplementary assumption must be done: Using the approximation technique based on a semi-discrete iterative scheme with respect to time, it can be proved that the regularized problem (P 0,τ ) has one solution h in W 1,∞ [0, T ] ; W 1,p 0 (Ω) . A detailed proof may be found in [4]. The main step of the proof relies on N.G. Meyers theorem on the regularity of elliptic equations. In a first time, we obtain a result in a space W 1,p 0 (Ω) with some 2 < p ≤ p, p not explicitly given which depends a priori on each iteration and is very close to 2. In a second time, we observe that the coefficients of the elliptic system are continuous on Ω, as are defined from a composing of a function of W 1,p 0 (Ω) and of a continuous function on b. Then the classical results of J. Nečas [12], [33] on the elliptic equations with continuous coefficients and smooth data may be applied stating that p = p uniformly.
We then have the following partial uniqueness result which links the relaxation parameter τ to the more or less steep initial topography introducing a critical value τ * : Proof. Consider h 1 any suitable solution to the problem (P 0,τ ) obtained by another method and h 2 the one given in this section, resulting from the same initial value h 0 ∈ W 1,p 0 (Ω), for a given p > 2. We define u = h 1 − h 2 and w = For almost all t ∈ ]0, T [, we substract the two equations verified by each solution to obtain The integral I can be evaluated in the following way ( > 0 to be defined) by Young's inequality: and moreover, applying Cauchy-Schwarz inequality yields Hence, applying Hölder's inequality, Thanks to the Rellich-Kondrachov embedding theorem, a positive constant C (p, Ω) exists such that Hence, we have: We then define (Ω)) , ε > 0. and so, for almost all t in ]0, T [, we obtain: We remark that: Provided that ε is chosen at ε 0 , the expression (1 − C (ε 0 )) is strictly positive once (Ω)) , and then the Gronwall's lemma leads to the uniqueness property.
Consider now h 1 and h 2 solutions to the problem (P 0,τ ) resulting from the initial values h 0,1 and h 0,2 . Assume moreover that (Ω)) .
It comes, for almost all t ∈ ]0, T [, that and the above proof ensures the local stability property.
Remark 1. In practice, proposition 4 is of real relevance for the following reason. We can squeeze more information from the previous calculations and particularly deduce this estimate on h: where the norm of h 0 in W 1,p 0 (Ω) is here the norm resulting from Poincaré inequality. Equivalently, only the directional slopes are taken in to account through the term ∇h 0 .
As a consequence, the critical relaxation time can be evaluated a priori from the inital known topography. The steeper the terrain is, the longer this respond time is.
In the case of sedimentation models (h is then limited by the constraint ∂h ∂t ≥ 0 a.e. in Q), due to regularizing effects, a subdivision of the interval [0, T ] in subintervals [t i , t i+1 ] may be introduced. The uniqueness is obtained separately on each sub-interval by iterations for critical decreasing values of τ * i . This enables to loosen progressively the constraint of the critical relaxation threshold. 6. Concluding remarks and prospects. Models describing the evolution of erosion and sedimentation in sedimentary basins have been developed recently at the l'Institut Français du Pétrole (IFP). They presents a real originality and use new mathematical formulations which are not always well controlled. These mathematical formulations are based on partial differential inequalities which do not inter, to our knowledge, in the framework of G. Duvaut and J.-L. Lions works about inequalities in Mechanics and Physics [15].
The idealized problems of instantaneous equilibrium are ill-posed without a supplementary criterion, not clearly defined up to now, for the selection of the flow limiter which make compatible the mass conservation law and the maximum erosion constraint.
To overcome this difficulty, the approach proposed is based on the introduction of small parameters (relaxation times) which traduce regularizing kinetic effects and leads to well-posed systems. Equivalently, it comes to traduce the fact that the equilibrium is attempted after a period of time. Passing at the limit when these parameters tend to zero, with enough a priori estimations, according to a general method, we expect to obtain a limit which can be used as a definition of solution of the initial ill-posed problem and also a supplementary discriminatory property. By this way, we have established a general result of existence and focused on models with sedimentation only. This situation is very singular in the family of erosion-sedimentation models. We have established a result of uniqueness and local stability from a critical value of a relaxation time. This value mainly depends on the smooth or steep state of the initial topography. However, for the moment, we do not know how to consider the case where the erosion constraint and the flow limiter are governed by a graph instead of a continuous function. This results comes from the mathematical difficulty to pass at the limit in a product of two terms converging in the sense of the weak topology, in the absence of complementary information on the compactness (one may use the Young measure theory [5]). Many open issues still subsist about geological models for stratigraphic reconstructions which are increasingly used in the framework of the petroleum engineering.
At the fact that the stratigraphic models describe slow deformations that can be followed in their movement, the Lagrangian coordinates system may present an alternative (see [14] for a situation with strong formal analogies in climatology: the socalled Budyko energy balance model governed by the Heaviside maximal monotone graph; the passage to Lagrangian coordinates allows to watch the trajectory of each of the fluid particles [36]). These matters are really issues whose interest is both industrial and academic to construct forward models of the stratigraphy of sedimentary basins.