Aurélien Rizk

Grégory Batt

François Fages

Sylvain Soliman

Continuous valuations of temporal logic

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Temporal logics [START_REF] Emerson | Temporal and modal logic[END_REF][START_REF] Clarke | Model Checking[END_REF] have proven useful as specification languages for describing the behavior of a broad variety of systems ranging from electronic circuits to software programs, and more recently biological systems in either boolean [START_REF] Eker | Pathway logic: Symbolic analysis of biological signaling[END_REF][START_REF] Chabrier | Symbolic model checking of biochemical networks[END_REF], discrete [START_REF] Bernot | A fruitful application of formal methods to biological regulatory networks: Extending thomas' asynchronous logical approach with temporal logic[END_REF], stochastic [START_REF] Calder | Analysis of signalling pathways using the continuous time Markov chains[END_REF][START_REF] Heath | Probabilistic model checking of complex biological pathways[END_REF] or continuous [START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF][START_REF] Antoniotti | Model building and model checking for biochemical processes[END_REF][START_REF] Chabrier | Symbolic model checking of biochemical networks[END_REF][START_REF] Batt | Validation of qualitative models of genetic regulatory networks by model checking : Analysis of the nutritional stress response in Escherichia coli[END_REF] settings.

Because temporal logics allow us to express both qualitative (e.g. some protein is eventually produced) and quantitative (e.g. a concentration exceeds 10) information about time and systems variables, they provide a powerful specification language in comparison with the essentially qualitative properties considered in dynamical systems theory (e.g. multistability, existence of oscillations) or with the exact quantitative properties considered in optimization theory (e.g. curve fitting). In particular, these logics are well suited to the increasingly quantitative, yet incomplete, uncertain and imprecise information now accumulated in the field of quantitative systems biology.

Such a use of temporal logics relies on a logical paradigm for systems biology which consists in making the following identifications [START_REF] Fages | Formal cell biology in BIOCHAM[END_REF]: biological model = transition system biological properties = temporal logic formulae biological validation = model-checking

This paradigm has been used in many applications, first for querying large interaction maps such as Kohn's map of the cell cycle [START_REF] Chabrier | Symbolic model checking of biochemical networks[END_REF][START_REF] Kohn | Molecular interaction map of the mammalian cell cycle control and DNA repair systems[END_REF] or gene regulatory networks [START_REF] Batt | Validation of qualitative models of genetic regulatory networks by model checking : Analysis of the nutritional stress response in Escherichia coli[END_REF], then for specifying the biological properties known or inferred [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF] from experiments, validating models, discriminating between models and proposing new biological experiments [START_REF] Bernot | A fruitful application of formal methods to biological regulatory networks: Extending thomas' asynchronous logical approach with temporal logic[END_REF], finding parameter values [START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF], and estimating robustness [START_REF] Batt | Robustness analysis and tuning of synthetic gene networks[END_REF]. An important limitation of this approach is however due to the logical nature of temporal logic specifications and their boolean interpretation. A yes/no answer to a temporal logic query does not provide indeed any information on how far we are from satisfaction, or how to guide the search to satisfy a formula. A measure of how close a model is to satisfy a property is needed.

In this paper, we define a continuous violation degree that quantifies how far from satisfaction an LTL formula is in a given model. In order to accommodate the various kinds of quantitative models defined by either ordinary or stochastic differential equation systems [START_REF] Segel | Modeling dynamic phenomena in molecular and cellular biology[END_REF][START_REF]System Modeling in Cellular Biology: From Concepts to Nuts and Bolts[END_REF], rule-based languages like SBML [START_REF] Hucka | The systems biology markup language (SBML): A medium for representation and exchange of biochemical network models[END_REF] or BIOCHAM [START_REF] Fages | Formal cell biology in BIOCHAM[END_REF], hybrid Petri nets [START_REF] Gilbert | A unifying framework for modelling and analysing biochemical pathways using Petri nets[END_REF][START_REF] Matsuno | Hybrid Petri net representation of gene regulatory network[END_REF], stochastic process calculi [START_REF] Priami | Application of a stochastic name passing calculus to representation and simulation of molecular processes[END_REF][START_REF] Phillips | A correct abstract machine for the stochastic pi-calculus[END_REF], etc..., we represent the behavior of the system simply by numerical traces [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF][START_REF] Nickovic | Amt: a property-based monitoring tool for analog systems[END_REF][START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF][START_REF] Antoniotti | Model building and model checking for biochemical processes[END_REF], so our method is rather general. Our notion of violation degree is then used for two applications in systems biology: the search of kinetic parameter values in a model, and the quantitative estimation of the robustness of a model by adapting the general framework of Kitano [START_REF] Kitano | Towards a theory of biological robustness[END_REF] to our temporal logic setting.

Section 2 presents the quantifier free fragment of first-order linear time logic with constraints over the reals, QFLTL(R), studied in [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF] and used in this paper. Section 4 defines a real-valued degree of satisfaction of a (QF)LTL formula on a trace using an aggregation function which composes the distances between the validity domains of some variables in a QFLTL formula and their objective value in the given specification.

Section 5 shows how such a continuous degree of satisfaction of an LTL formula can be used as a fitness function in local search methods for finding kinetic parameter values satisfying a temporal logic specification. We describe a general local search method and use the state-of-the-art Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] to evaluate the method on models of the budding yeast cell cycle with 8 parameters and of the MAPK signaling cascade with 30 parameters and 7 unknown initial conditions.

In section 6 we propose a definition for the robustness degree of a property w.r.t. a model and a set of perturbations weighted by probabilities. This definition is inspired by the abstract definition of robustness proposed by Kitano for systems biology [START_REF] Kitano | Towards a theory of biological robustness[END_REF]. We develop it here in our temporal logic setting and illustrate its relevance by applying it to the previous model of the cell cycle. In [START_REF] Rizk | A general computational method for robustness analysis with applications to synthetic gene networks[END_REF], this method is applied with some more details to the design of synthetic gene networks.

Preliminaries on Temporal Logics with Constraints over the Reals

LTL(R)

The Linear Time Logic LTL is a temporal logic [START_REF] Clarke | Model Checking[END_REF] that extends classical logic with modal operators for qualifying when a formula is true in trace, i.e. an infinite sequence of timed states. Temporal operators include X ("next", for at the next time point), F ("finally", for at some time point in the future), G ("globally", for at all time points in the future), U ("until") and W (" weak until"). These operators enjoy some simple duality properties,

¬ X φ = X ¬φ, ¬ F φ = G ¬φ, ¬ G φ = F ¬φ, ¬(ψ U φ) = (¬φ W ¬ψ).
These properties make it possible to eliminate negations by descending them to atomic propositions. Furthermore, we have φ W ψ = (φ U ψ) ∨ G φ, hence LTL can be defined with { X F , G , U } as basis of temporal operators.

A version of LTL with constraints over the reals, named LTL(R), has been proposed in [START_REF] Antoniotti | Model building and model checking for biochemical processes[END_REF][START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF] to express temporal properties about molecular concentrations. The atomic formulae of LTL(R) are formed with inequality relations and arithmetic operators over the real values of molecular concentrations (denoted by [A], ...), their derivatives (denoted by d[A]/dt, ...) and time. Since the relation symbols are closed by negations, and negations and implications can be eliminated by propagating the negations down to the atomic constraints in the formula, we assume that all LTL(R) formulae are in negation free normal form. The precise syntax of LTL(R) is given in Table 1.

Formula ::= Atom | Formula ∧ Formula | Formula ∨ Formula | X Formula | F Formula | G Formula | Formula U Formula Atom ::= Term Op Term Op ::

= < | > | ≤ | ≥ Term Real | Term + Term | Term -Term | -Term | Term × Term | Term / Term | Term ˆTerm Real ::= float | [molecule] | d[molecule]/dt | Time Table 1 Syntax of LTL(R) formulae.
For instance, F([A]>10) expresses that the concentration of A eventually gets above the threshold value 10. G([A]+[B]<[C]) expresses that the concentration of C is always greater than the sum of the concentrations of A and B. Oscillation properties, abbreviated as oscil(M,K), are defined as a change of sign of the derivative of M at least K times:

F((d[M]/dt > 0) ∧ F((d[M]/dt < 0) ∧ F((d[M]/dt > 0)...)))
LTL(R) formulae are interpreted over infinite traces of the form

(< t 0 , x 0 , dx 0 /dt >, < t 1 , x 1 , dx 1 /dt >, ...)
which gives the concentration values x i of the molecules, and the values of their first derivatives dx i /dt, at discrete time points t i . We shall consider infinite traces obtained from finite traces by completing them with a loop on the last state. For instance, in a model described by a system of ordinary differential equations (ODE), and under the hypothesis that the initial state is completely defined, numerical integration methods (such as Runge-Kutta or Rosenbrock method for stiff systems) compute a finite simulation trace. To extend it to an infinite trace, a loop is added on the last state with the assumption that the finite time horizon considered for the numerical integration is sufficiently large to check the properties at hand.

It is worth noticing that the semantics of the "next" operator X refers to the next time point on the trace, and that in adaptive step size integration methods of ODE systems, the step size t i+1 -t i is not constant but determined through an estimation of the error made by the discretization. The X operator can be used to detect local optima in a trace for instance.

Formally, the truth value of an LTL(R) formula in a trace π is given in Table 2.

π |= φ iff s 0 |= φ s i |= α
iff α is an atomic formula and α is true at state s i ,

s i |= φ ∧ ψ iff s i |= φ and s i |= ψ, s i |= φ ∨ ψ iff s i |= φ or s i |= ψ, s i |= X φ iff s i+1 |= φ, s i |= F φ iff there exists j ≥ i s.t. s j |= φ, s i |= G φ iff for all j ≥ i, s j |= φ, s i |= φ U ψ iff there exists j ≥ i s.t. s j |= ψ and s k |= φ for all i ≤ k < j.
Table 2 Inductive definition of the truth value of an LTL(R) formula in a trace π = (s 0 , s 1 , . . .).

These truth values can be computed on traces by a model-checking algorithm [START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF].

QFLTL(R)

In [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF], the quantifier free fragment of first-order LTL(R), named QFLTL(R), has been considered for the purpose of analyzing numerical data time series in temporal logic and computing automatically LTL(R) specifications from experimental traces. Syntactically, QFLTL(R) simply adds variables to the terms:

Real ::= variable | float | [molecule] | d[molecule]/dt | ...
For instance, the QFLTL(R) formula G([A] < x) expresses the constraint that x is always greater than the concentration of A, or equivalently, that x is greater than the maximum of [A]. Obviously, an LTL(R) formula is a QFLTL(R) formula without variable.

As usual, the semantics of a QFLTL(R) formula containing variables is defined by its ground instances which are LTL(R) formulae. Given a trace π and a QFLTL(R) formula φ over a vector x of k real-valued variables, the constraint satisfaction problem, ∃v ∈ R k φ[v/x], is the problem of determining the valuations v of the variables for which the formula φ is true. The domain of validity of a formula φ over variables x with respect to the trace π is defined as

D π,φ = {v ∈ R k | π |= φ[v/x]}
In [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF], LTL(R) model-checking has been generalized to a QFLTL(R) constraint solving algorithm which computes the exact domain of validity D π,φ for any QFLTL(R) formula φ. This algorithm uses the following characterization of validity domains: Proposition 2.1 (Inductive definition of validity domain) Let π = (s 0 , s 1 , . . . , s n ) be a finite trace and φ be a QFLTL formula with k variables, we have

• D π,φ = D s 0 ,φ , • D s i ,α = {v ∈ R k | s i |= α[v/x]} for an atomic proposition α, • D s i ,φ∧ψ = D s i ,φ ∩ D s i ,ψ , • D s i ,φ∨ψ = D s i ,φ ∪ D s i ,ψ , • D s i ,Xφ = D s i+1 ,φ , • D s i ,Fφ = ∪ j≥i D s j ,φ , • D s i ,Gφ = ∩ j≥i D s j ,φ , • D s i ,φUψ = ∪ j≥i (D s j ,ψ ∩ ∩ k∈[i,j-1] D s k ,φ ).
The QFLTL(R) constraint solving algorithm computes the validity domains of the variables of the formula in each time point of the trace, by double induction on the subformulae of φ and on the states of π.

Implementation and Complexity

The atomic propositions considered so far for QFLTL(R) include non-linear inequalities which would be too general to be computed efficiently without approximations. To get a handle on the complexity of QFLTL constraint solving, we first describe a restricted fragment of QFLTL formulae leading to finite unions of orthotopes as domains, and then describe our implementation of the linear constraint fragment of QFLTL leading to finite unions of polytopes as domains.

Orthotopes

We consider here the fragment of QFLTL(R), named QFLTL(R box ), in which each atomic formula contains at most one variable:

Atom ::= Term Op Term | Term Op variable Real ::= float | [molecule] | d[molecule]/dt | ...
The validity domain of a conjunction of such constraints in one variable is an interval over the reals, and in higher-dimensions over v variables, is an orthotope, i.e. a cartesian product of v intervals over R. These orthotopes are computed in the base case of the inductive definition 2.1 of validity domains. For the other cases, the QFLTL(R box ) constraint solving algorithm computes finite intersections and unions of validity domains. As a finite intersection of orthotopes is a finite union of orthotopes, the domains computed by the algorithm are always finite unions of orthotopes.

In [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF], we have shown that the number of orthotopes needed to represent the validity domain of a QFLTL(R box ) formula may be exponential in the number of variables, but is polynomial in the size of the formula and in the size of the trace: Proposition 3.1 The validity domain of a QFLTL(R box ) formula of size f containing k variables, on a trace of length n is a union of at most (nf ) 2k orthotopes.

One example necessitating (nf ) k orthotopes is given in [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF].

Polytopes

The linear fragment QFLTL(R lin ) generalizes QFLTL(R box ) by allowing linear constraints as atomic propositions. The validity domains of such linear constraints are polytopes, and since the intersection of two polytopes is a polytope, the QFLTL(R lin ) constraint solving algorithm computes finite unions of polytopes.

Proposition 3.2

The validity domain of a QFLTL(R lin ) formula of size f on a trace of length n is a union of at most 2 nf polytopes.

Proof. As a polytope is represented by a set of linear constraints, let us consider the number of possible linear constraints appearing in the validity domain D π,φ of a formula φ of size f on a trace π of size n. These constraints are only generated by the base case of the inductive definition 2.1 of validity domains. Each constraint in φ can be evaluated on each time point of the trace, thus creating at most n different constraints. As intersections and unions do not generate new constraints and as the number of constraints in φ is less than f , the number of constraints in D π,φ is less than nf .

Because a polytope is defined by an arbitrary number of constraints, with at most nf constraints one can form at most 2 nf polytopes. Therefore, validity domain D π,φ is a union of at most 2 nf polytopes. This fragment QFLTL(R lin ) is what is currently implemented in BIOCHAM v8 [START_REF] Fages | BIOCHAM v8 user's manual[END_REF] using the Parma Polyhedra Library [START_REF] Bagnara | The Parma Polyhedra Library: Toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems[END_REF] for computing union and intersection of sets of polytopes. The performance evaluation given in Section 5 refers to this implementation, and indicates that the upper complexity bound is not reached in practice where most domains are orthotopes.

Continuous Satisfaction Degree of Temporal Logic Formulae

In order to evaluate numerically the fitness of a model w.r.t. a temporal logic specification, we show in this section how QFLTL(R) formulae can be evaluated with a continuous satisfaction degree in the interval [0, 1], instead of a Boolean value. This continuous valuation of temporal logic formulae opens up the field of model-checking to optimization, as illustrated in the following sections with applications to parameter optimization and robustness evaluation in systems biology.

Continuous Satisfaction of QFLTL formulae

Let π be a trace and φ be a QFLTL(R) formula that we wish to valuate continuously in the interval [0, 1]. For this, we introduce a QFLTL(R) pattern formula ψ(x 1 , ..., x k ) obtained by replacing some constants in φ by new variables {x 1 , ..., x k }: we have φ = ψ(v 1 , ..., v k ) for some instantiation of the variables by real values v 1 , ..., v k . The satisfaction degree of φ in π is then defined using the distance between the validity domain of the variables x 1 , ..., x k in ψ and the objective values (v 1 , ..., v k ) in R k . Definition 4.1 The violation degree vd(π, φ, ψ) of a QFLTL formula φ in a numerical trace π with respect to a pattern formula ψ(y) such that φ = ψ(v) for some real values v, is the Euclidean distance between v and the projection on variables y of domain Now, given a mathematical model of our system, let us suppose that the QFLTL constraint solving algorithm applied to ψ on simulation trace π computes D x π,ψ =] -∞, 15] as domain for variable x. Since v = 20 we get vd(π, φ, ψ) = 5, i.e. the violation degree is 5 since the compound reaches a maximum of 15 whereas the formula expresses that the threshold 20 be reached. The formula

(D π,ψ ), or +∞ if D π,ψ = ∅: vd(π, φ, ψ) = min v ∈D y π,ψ d(v , v) The satisfaction degree of φ in π w.r.t. ψ is sd(π, φ, ψ) = 1 1 + vd(π, φ, ψ)
ψ = F ([A] ≥ x) ∧ F ([A] ≤ y) ∧
x -y ≥ amp will allow us to reason on the amplitude of these variations. Since amp is the only variable of interest (the only one instantiated from ψ to φ) it is the only one used to define the objective that the amplitude should be at least 10. Now suppose that the constraint solving computes the following validity domains for x and y: D x π,ψ =] -∞, 15] and D y π,ψ = [10, +∞[. For formula φ, the maximum value of D x π,ψ represents the maximum value of [A] and the minimum value of D y π,ψ its minimum value in the trace. The domain for variable amp is D amp π,ψ =] -∞, 5] since we know that amp ≥ x -y, and thus, since the objective value is 10, we obtain vd(π, φ, ψ) = 5, i.e. the amplitude of the curve is 5 whereas we wanted it to be at least 10.

Note that if the trace π is such that φ is satisfied then vd(π, φ, ψ) = 0 since the objective value belongs to the validity domain D π,ψ . When φ is not valid on π, the violation degree vd provides a quantitative measure of its degree of non satisfaction. The use of this measure is illustrated in the following sections to improve parameter search for biological models and to define a quantitative notion of robustness of a system w.r.t. a temporal logic formula.

Abstraction of constants by variables in LTL Formulae

Since an LTL(R) formula is a special case of QFLTL(R) formula, containing no variable, the definitions of the previous section directly apply to it. However, as it is more usual to write specifications in LTL rather than QFLTL, we present here an abstraction α from LTL to QFLTL formulae providing a default choice for the pattern formula ψ associated to an LTL specification formula φ. This abstraction α thus provides a definition of the satisfaction degree of an LTL formula in a trace, without referring explicitly to any QFLTL pattern formula.

An LTL specification φ (of the expected behavior of a system), can be transform in a QFLTL formula α(φ) by mapping the occurrences of its constants c 1 , . . . , c k (i.e. real numbers corresponding to concentration thresholds, amplitudes, etc.) appearing in φ, to distinct variables x 1 , . . . , x k . The reverse mapping σ = α -1 is an instantiation of the variables x 1 , . . . , x k to real values v 1 , . . . , v k such that φ = σ(ψ).

φ = F ([A] > 20) we obtain α(φ) = F ([A] > x)
, which corresponds to the ψ that was chosen since the threshold is the only occurrence of a constant in φ.

The same abstraction process can also be applied to QFLTL formulae, resulting in QFLTL formulae having more variables. We have φ

= F ([A] ≥ v) ∧ F ([A] ≤ w) ∧ v -w ≥ 10 and ψ = F ([A] ≥ v) ∧ F ([A] ≤ w) ∧ v -w ≥ amp.
φ is not an LTL formula but the above abstraction still provides α(φ) = ψ.

Obviously, such a systematic abstraction of constants by variables may produce nonlinear constraints and need be restricted in some cases. For instance in example S5 of section 5.2, some constants are kept in ψ 4 in order to fit a curve to given time points. The abstraction α may also be inadequate when one wants to rescale the weight between variables. For instance if φ = F ([A] > 20) ∧ F ([B] < 0.001), α will put the same weight on both distances, whereas a direct formulation with ψ would permit a rescaling, as in ψ = F ([A]/20 < x) ∧ F ([B]/0.001 > y). For these reasons, the continuous valuation degree has been defined for QFLTL specifications w.r.t. a QFCTL pattern formula in a completely general form in the previous section.

Parameter Optimization w.r.t. QFLTL(R) Properties

The violation degree provides a measure of how far a given numerical trace is from satisfying a temporal logic specification. This continuous measure of satisfaction opens up the field of model-checking to optimization, and can be used to guide the search for models satisfying a given specification. In particular, in dynamical models of molecular interaction networks, the search of kinetic parameters w.r.t. LTL specifications by scanning the space of parameter values described in [START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF], can be replaced by a much more efficient local search and evolutionary methods which explore neighborhoods of parameter sets and progress by choosing the ones which minimize the violation measure.

Optimization Algorithm

Let us consider a (QF)LTL formula φ, a QFLTL formula ψ, such that φ = σ(ψ) as in section 4.1, an SBML/BIOCHAM reaction model with initial conditions and known parameter values, a set of unknown parameters to explore and for each of those an interval of search. We consider the problem of finding a set of values of the unknown parameters such that the violation degree of the corresponding trace π obtained by numerical simulation is vd(π, φ, ψ) = 0.

A generic optimization algorithm for parameter search can be described as follows:

Algorithm 5.1 1. Set the current point in the parameter space to a random point belonging to the provided search box, compute a numerical simulation with trace π and the corresponding violation degree vd(π, φ, ψ); 2. if vd = 0 jump to 5.

3. for each point in a defined neighborhood of the current point, compute a trace and its violation degree;

4. based on the violation degrees of the neighbors, determine the next point of the iteration, set the current point to this point, update current vd and go to 2.

Return the current point in the parameter space.

This procedure can be interrupted after a given number of steps, returning the best parameter set (minimizing the violation degree). It can also be restarted with a new initial point (step 1) several times in order to diversify the search. A naive method would be to define as neighborhood of the current parameter state the parameter sets obtained by modifying one parameter by values ±δ ; and to choose as next parameter set the best neighbor.

More efficient instances of the previous algorithm can be obtained however, by combining state-of-the-art nonlinear optimization methods with the computation of our violation degree used as a black box fitness function. In the following sections, we use the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) of Hansen and Ostermeier [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]. This method uses a probabilistic definition of the neighborhood, and stores information in a covariance matrix in order to replace the approximate gradient and Hessian of a quasi-Newton method by an evolutionary algorithm.

Besides a specified fitness function, CMA-ES performs a parameter search given an initial solution, stop and restart criteria and a given search space. The search stops either when a given number of violation degrees have been computed or when the violation degree gets below a given threshold. When differences between successive evaluation of the violation degree are below a given value, i.e. the violation degree landscape is too flat, the search restarts from a randomly chosen solution inside the search space. The search space can be defined by intervals or by sets of linear constraints on the searched parameters. In this case, the Parma Polyhedra library [START_REF] Bagnara | The Parma Polyhedra Library: Toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems[END_REF] is used to compute the polytope representing the search space, check that it is not empty, and randomly pick a value in this polytope for restarts.

Evaluation on Cell Cycle Models

In this section we present the application of the parameter search method outlined above to the budding yeast cell cycle model of [START_REF] Chen | Kinetic analysis of a molecular model of the budding yeast cell cycle[END_REF]. This model displays how proteins cdc2 and cyclin interact to form the heterodimer Cdc2-Cyclin~{p1,p2} known as maturation promoting factor (MPF) and playing a key role in the control of mitotic cycles. The model consists of the following BIOCHAM [START_REF] Fages | BIOCHAM v8 user's manual[END_REF] reaction rules with mass action (MA) kinetics and an autocatalytic kinetics for the dephosphorylation of the mitosis promotion factor Cdc2-Cyclin p1:

AUTOCAT(k4)= k4*[Cdc2-Cyclin~{p1,p2}]*[Cdc2-Cyclin~{p1}]^2.
We use as reference point k T yson the values of the kinetic parameters determined in [START_REF] Chen | Kinetic analysis of a molecular model of the budding yeast cell cycle[END_REF]. The simulation for k T yson of the system of ODEs extracted from these rules, given in appendix, is displayed in Figure 1. The total amount of cyclin presents oscillations of period 35 while MPF exhibits activity peaks with same period.

Using the optimization method CMA-ES together with our violation degree as a parameter search method we wonder whether it is possible to find values of the kinetic parameters corresponding to higher MPF peaks or oscillations with higher amplitudes or shorter periods.

Search problem S1 : higher MPF peaks (2 parameters unknown)

Two parameters, k4 and k6, have been found in [START_REF] Chen | Kinetic analysis of a molecular model of the budding yeast cell cycle[END_REF] to play a particular role for the existence of oscillations. Depending on their values the system exhibits either a steady state behavior or limit cycle oscillations. We wonder whether it is possible to obtain higher MPF peaks by changing values of k4 and k6 only, all other parameters remaining at the value k T yson chosen in [START_REF] Chen | Kinetic analysis of a molecular model of the budding yeast cell cycle[END_REF]. More precisely, we want to reach MPF peaks of 0.3 at least, the maximum amount of MPF for k T yson being 0.19.

Therefore we define the LTL specification :

φ 1 = F ([MPF ] > 0.
3) with the corresponding QFLTL formula being :

ψ 1 = F ([MPF ] > max)
The variable space associated to φ 1 and ψ 1 is R and corresponds to the sole variable max. The objective is v = 0.3, i.e the target peak value of MPF is 0.3. We have been able to find valid parameter values, denoted k * T yson , satisfying vd(π, φ 1 , ψ 1 ) = 0 where π is the corresponding simulated trace (see Figure 1b). k * T yson is given in Table 3. As the plot shows, for these parameter values essential features of the curve, especially repeated MPF peaks, are conserved although it was not enforced by the specification. In particular, a constantly growing amount of MPF would have also resulted in a null violation degree of this formula.

All computations have been performed on an Intel Core 2 Duo 2GHz with 2GB RAM. Note that as the optimization method CMA-ES uses a probabilistic neighborhood two consecutive runs can yield different results. In this example answers are typically obtained in less than 30 s after around 150 numerical simulations and violation degree computations.

Search problem S2 : amplitude of MPF oscillation (2 parameters unknown)

In this example we refine the previous query by constraining the minimum level of MPF. We search for k4 and k6 values that preserve at least two periods of MPF oscillations having same amplitudes as those observed for k T yson .

In order to specify that the amplitude is at least 0.19, we use the variable space R corresponding to only one variable, amp, with objective value v = 0.19, i.e the target amplitude is 0.19. This value corresponds to the amplitude obtained for k T yson .

φ 2 = F ([MPF ] > max ∧ F ([MPF ] < min ∧ F ([MPF ] > max ∧ F ([MPF ] < min)))) ∧ max -min > 0.19 ψ 2 = F ([MPF ] > max ∧ F ([MPF ] < min ∧ F ([MPF ] > max ∧ F ([MPF ] < min)))) ∧ max -min > amp
Starting from a different value k 2 for k4 and k6, we try to recover the behavior of k T yson . We found such parameters, given in Table 3 and referred to as k * 2 . To illustrate the path followed during the search from k 2 to k * 2 we computed the violation degree landscape in the k4, k6 parameter space. The resulting landscape is displayed in Figure 2. Note that as all constants of the formula have been abstracted by variables, the violation degree can only be finite. In particular when no oscillations are present in the trace amp will be equal to 0, thus leading to a violation degree of 0.19. Regions where the violation degree is 0.19 correspond to regions of steady state behavior whereas regions with a violation degree between 0 and 0.19 correspond to regions of oscillations.

Under mild assumptions Tyson determined linear equations defining a region in the k4, k6 plane where oscillations occurs, also represented in Figure 2. Our results are fully consistent with his analytical analysis, and provide more information on the amplitude of oscillation w.r.t. parameters k4 and k6.

Search problem S3 and S4 : amplitude and period of oscillations (all 8 parameters unknown)

To illustrate the scalability of the method we carry out two parameter searches on all 8 parameters of the model. The first one (problem S3) is the same query as above with formulae φ 2 and ψ 2 but with all parameters unknown.

The second one is a more complex query used to find shorter oscillation periods of Cdc2: to specify that the target period is 20, we use the variable space R corresponding to the variable per with objective value v = 20.

φ 3 = F (d([Cdc2 ])/dt < 0 ∧ X(d([Cdc2 ])/dt > 0 ∧ Time > t1 ∧ X(F (d([Cdc2 ])/dt > 0 ∧ X(d([Cdc2 ])/dt < 0 ∧ T ime < t2)) ∧ t2 -t1 < 20 ψ 3 = F (d([Cdc2 ])/dt < 0 ∧ X(d([Cdc2 ])/dt > 0 ∧ T ime > t1 ∧ X(F (d([Cdc2 ])/dt > 0 ∧ X(d([Cdc2 ])/dt < 0 ∧ T ime < t2)) ∧ t2 -t1 < per
Search problem S3 starts from parameter values k 3 satisfying the constraints on their order of magnitude given in [START_REF] Chen | Kinetic analysis of a molecular model of the budding yeast cell cycle[END_REF]. k 3 does not give rise to oscillations. Search problem S4 starts from k 4 = k T yson . In both cases parameter values are found satisfying the query (in around 350 s for S3 and 30 s for S4). Results are given in Table 3.

Evaluation on MAPK Signal Transduction Model

The MAPK signal transduction model of [START_REF] Levchenko | Scaffold proteins may biphasically affect the levels of mitogen-activated protein kinase signaling and reduce its threshold properties[END_REF] is used to test the scalability of the parameter search method on a larger model, and illustrates the possibility of getting oscillations in a purely directional cascade of reactions [START_REF] Qiao | Bistability and oscillations in the huang-ferrell model of mapk signaling[END_REF]. The MAPK signalling network consists of a cascade of phosphorylation reactions from input membrane receptor activated molecule RAF, to output doubly phosphorylated molecule MAPK p1,p2 which enters the nucleus. The cascade of phosphorylations is depicted in Figure 3.

The model consists of the following reaction rules with mass action (MA) law kinetics:

(MA(k1), MA(k2)) for RAF + RAFK <=> RAF-RAFK.

(MA(k3),MA(k4)) for RAF~{p1} + RAFPH <=> RAF~{p1}-RAFPH.

(MA(k5),MA(k6)) for MEK~$P + RAF~{p1} <=> MEK~$P-RAF~{p1} where p2 not in $P. (MA(k7),MA(k8)) for MEKPH + MEK~{p1}~$P <=> MEK~{p1}~$P-MEKPH. (MA(k9),MA(k10)) for MAPK~$P + MEK~{p1,p2} <=> MAPK~$P-MEK~{p1,p2} where p2 not in $P. (MA(k11),MA(k12)) for MAPKPH + MAPK~{p1}~$P <=> MAPK~{p1}~$P-MAPKPH. MA(k13) for RAF-RAFK => RAFK + RAF~{p1}. MA(k14) for RAF~{p1}-RAFPH => RAF + RAFPH. MA(k15) for MEK~{p1}-RAF~{p1} => MEK~{p1,p2} + RAF~{p1}. 3 Resulting parameter values for search problems S1, S2, S3 and S4. All these search problems stop when the target specification is satisfied, i.e. when the violation degree reaches 0. Computation times are the average values evaluated on 3 runs. Total simulation time is the total time spent for numerical simulations, while vd computation time is the total time spent to compute violation degrees. In this example, we investigate the use of our parameter search method as a curve fitting tool at specific time points, on 22 parameter values. In order to express the classical distance between two curves at time points 30 and 60 for instance, we use the following pattern of formulae :

k 2 k * 2 k 3 k * 3 k 4 k * 4 k1 1.50e-2 1.50e-2 1.50e-2 1.50e-2 1.00e-2 1.14e-2 1.50e-2 2.41e-2 k3 2.00e2 2.00e2 2.00e2 2.00e2 1.00e2 1.13e2 2.00e2 2.83e2 k4p 1.80e-2 1.80e-2 1.80e-2 1.80e-2 1.00e-2 8.77e-3 1.80e-2 2.24e-2 k4 1.80e2 8.99e2 2.00e1 1.94e2 1.00e2 1.82e2 1.80e2 2.28e2 k6
k4 k6 . . k * 2 . k A . k 2 . k B . k C .
ψ 4 = G( T ime = 30 → [MEK-RAF ˜{p1}] = x ∧ T ime = 60 → [MEK-RAF ˜{p1}] = y)
The parameter space of this formula is R 2 and is defined by the two variables x and y. We set the objective v to the target values of [MEK-RAF ˜{p1}] at time 30 and 60, which defines completely φ 4 .

Note that if the time points given in this formula are not present in the numerical simulated trace of the model the formula will always be true, resulting in a null violation degree. In order to effectively use this formula for curve fitting, we need to ensure that numerical simulations contain these specific time points. This can be done by computing concentration values at these time points by linear interpolation of neighboring values. Besides, this formula can be extended to any number of time points and molecules in order to perform a complete curve fitting, if it is relevant.

This Note that a search problem with such dimensionality could not be solved with parameter scanning, which has exponential time complexity in the number of parameters, and which is what is done with a binary valuation of temporal logic specifications [START_REF] Calzone | Machine learning biochemical networks from temporal logic properties[END_REF]. By allowing the use of any continuous optimization algorithm, we provide a parameter search method with respect to temporal specification as scalable as the continuous optimization algorithm used. In our application we use CMA-ES [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] which has been shown to be applicable for searching up to one hundred parameters.

Here, the computation time is more dependent on the type of problem (formula used and initial values of the parameters) and on the landscape of the violation degrees than on the number of parameters.

Search problem S6 : find oscillations (30 kinetic parameters and 7 initial conditions unknown)

In [START_REF] Qiao | Bistability and oscillations in the huang-ferrell model of mapk signaling[END_REF], oscillations have been found in the MAPK cascade model of [START_REF] Levchenko | Scaffold proteins may biphasically affect the levels of mitogen-activated protein kinase signaling and reduce its threshold properties[END_REF] although this model does not contain any negative feedback reaction. This does not contradict Thomas' necessary condition for sustained oscillations as such a purely directional cascade does contain negative feedback in its influence graph as shown in [START_REF] Fages | From reaction models to influence graphs and back: a theorem[END_REF] and analyzed in [START_REF] Ventura | A hidden feedback in signaling cascades is revealed[END_REF]. However, to know whether these negative circuits in the influence graph are functional, one needs to search for kinetic parameter values and initial conditions that exhibit sustained oscillations.

Just by defining the following formulae:

φ 5 = F ([MAPK ˜{p1, p2}] > max ∧ F ([MAPK ˜{p1, p2}] < min)) ∧ max -min > 0.5 ψ 5 = F ([MAPK ˜{p1, p2}] > max ∧ F ([MAPK ˜{p1, p2}] < min)) ∧ max -min > amp
that use the variable space R for the single variable amp < max -min, and ask that the amplitude be at least 0.5, setting an objective value v = 0.5, parameter values leading to sustained oscillations, such as the ones depicted in Figure 5, were found in a few minutes. The notion of satisfaction degree allows us to define in a mathematically precise way a degree of robustness of a system's behavior property described in temporal logic, with respect to a set of system's perturbations, and estimate it computationally. This robustness degree is defined as the average satisfaction degree of the property of interest over all admissible perturbations, possibly weighted by their probabilities: Definition 6.1 Let P be a set of perturbations, prob(p) be the probability of perturbation p, π(p) be the timed trace of the system under perturbation p ∈ P . The robustness degree R φ,ψ,P of properties φ and ψ with respect to perturbation P is the real value R φ,ψ,P = p∈P sd(π(p), φ, ψ) prob(p) dp If the set of perturbations is finite (eg, gene knock outs), the robustness degree is simply the inverse of a finite weighted sum and can be exactly computed. If the set of perturbations is infinite, the robustness degree can be estimated by computing the violation degree between the behavior of the perturbed system π(p) and the specification φ for sufficiently many perturbations.

This definition is an adaptation of the general definition given by Kitano [START_REF] Kitano | Towards a theory of biological robustness[END_REF] to our temporal logic setting. It is described in more detail with applications to synthetic genetic networks in [START_REF] Rizk | A general computational method for robustness analysis with applications to synthetic gene networks[END_REF].

Evaluation on Cell Cycle Model

Using the same cell cycle model as in section 5.2, we compare the robustness of oscillation properties with regard to perturbations of parameter values k4 and k6 for different points in the 4 Total CPU time, simulation time and violation degree computation time (s) for search problems S5 and S6. As the violation degree associated to curve fitting problem S5 is always strictly positive, the search is stopped as soon as vd gets below arbitrary threshold value 0.005. S6 stops when the violation degree gets to 0. parameter space.

We consider that parameter values for k4 and k6 are normally distributed around their reference value with coefficient of variation equal to 0.2. We also enforce that k4 ≥ 0 ∧ k6 ≥ 0. We examine the robustness of the property expressed by φ 2 and ψ 2 , that is, MPF oscillations are of amplitude at least 0.19.

The robustness degree of this property is compared for three different values of k4 and k6. These three points in the parameter space of k4 and k6 are indicated by the three points k A , k B and k C in Figure 2. In all cases, the estimation of the robustness degree is done by computing the mean value of the violation degree for 500 samples.

The estimated degree of robustness for parameters k A , k B and k C are respectively 133, 12.9 and 13.5. This is consistent with the location of points k A , k B and k C . Perturbations around point k A have high probabilities of staying in the region satisfying the specification whereas perturbations around point k B have high probabilities of moving the system to the region with no oscillation. k C is more robust than k B even though, as opposed to k B , its violation degree is non null. This can be explained by the abrupt transition between oscillating and non oscillating regions near k B compared to the smoother transition near k C .

The robustness degree can be estimated for perturbations on any number of parameters. For instance, by computing a robustness estimate for perturbations on all parameters, with coefficient of variation 0.2 for specification φ 2 , ψ 2 and parameter values k T yson and k 3 , the estimated robustness degrees for k T yson and k 3 are 20.7 and 27.1 respectively. This indicates that the oscillations are more robust to variations of the parameters values for k 3 than for the parameters given in the original model of Tyson.

Related Work

Probabilistic temporal logics and probabilistic model checking have been used in systems biology [START_REF] Kwiatkowska | Using probabilistic model checking in systems biology[END_REF], e.g. for an analysis of a probabilistic model of the MAP kinase signaling cascade. However these techniques provide information on the probability that a given property is exactly satisfied. They thus provide no quantitative information on unsatisfied formulae and cannot be compared to the satisfaction degree presented in this paper.

More closely related to our continuous satisfaction degree are the linear metrics for quantitative transition systems defined in [START_REF] De Alfaro | Linear and branching metrics for quantitative transition systems[END_REF]. These metrics apply to traces and can be characterized by quantitative LTL formulae. LTL formulae are interpreted on the [0, 1] interval. However, no implementation is proposed, and the applicability of this approach to solving optimization and robustness problems is not discussed.

Fainekos and Pappas have also proposed a satisfaction degree for temporal logic specifications [START_REF] Fainekos | Robustness of temporal logic specifications[END_REF]. Although the two approaches share some similarities, a significant difference is that in [START_REF] Fainekos | Robustness of temporal logic specifications[END_REF] the satisfaction degree corresponds to a distance between a trace and the set of traces satisfying a formula, whereas in our case the violation degree corresponds to a distance between a formula and the set of formulae satisfied by the trace. One major advantage of our approach is that the dimensionality of the formula space (number of constants abstracted by variables) is generally much lower than the dimensionality of the trace space (trace length). Reasoning in low dimensional space may strongly affect the practical applicability of these methods. Note however both approaches, handling sets of traces [START_REF] Fainekos | Robust sampling for MITL specifications[END_REF][START_REF] Dang | Sensitive state space exploration[END_REF], and our approach, handling sets of formulae, are a priori compatible, and that their combination might combine their benefits.

Another notion of violation degree has been recently proposed by Donaldson and Gilbert [START_REF] Donaldson | A model checking approach to the parameter estimation of biochemical pathways[END_REF], also based on the definition of a validity domain of temporal logic formulae. However, the computation of the (finite) validity domain is made by sampling the formula space rather than by constraint solving.

Concerning robustness, in [START_REF] Dayarian | Shape, size and robustness: Feasible regions in the parameter space of biochemical networks[END_REF], Chaves and colleagues propose a quantitative measure of robustness corresponding to the volume of the set of valid parameters in the parameter space. This measure thus reflects the proportion of parameters that satisfy exactly the property, as opposed to our measure that represents how close to satisfying the property the system is for various parameters. These two measures provide complementary information on robustness. In [START_REF] Batt | Robustness analysis and tuning of synthetic gene networks[END_REF], robustness is similarly defined with respect to temporal logic specifications. However, it has a Boolean interpretation, since a property is defined as robustly satisfied by an ODE system if it is satisfied by the system for all possible perturbations. As stated earlier, obtaining a quantitative measure of robustness is more informative for many practical problems.

Conclusion

We have defined a continuous measure of satisfaction of an LTL(R) formula in a numerical trace and shown that it can be computed using the QFLTL(R) constraint solving algorithm of [START_REF] Fages | On temporal logic constraint solving for the analysis of numerical data time series[END_REF]. This measure is more informative that the Boolean interpretation of the formulae and can be used in many situations in systems biology to reason about numerical traces.

We have shown that this measure can be used as a fitness function in state-of-the-art optimization tools to efficiently guide the search of kinetic parameter values in biochemical reaction models -a central problem of computational systems biology. Several tenths of parameters could be found at the same time to satisfy non-trivial semi-qualitative semi-quantitative properties formalized in QFLTL(R).

We have also shown that this satisfaction measure can similarly be used to estimate the robustness of a model w.r.t. temporal logic specifications, in accordance to Kitano's notion of robustness for systems biology. The generalization of model-checking to temporal logic constraint solving which is at the basis of our computation method thus seems to open new research avenues for the use of temporal logics in systems biology, and more generally for design problems from temporal specifications.

Example 4 . 1

 41 Let us consider the LTL formula φ = F ([A] > 20) specifying it was observed in experiments that after some time the concentration of compound A becomes greater than 20. Take ψ = F ([A] > x) as QFLTL pattern formula and R as variable space indicating that we are interested in looking at the current values reached by [A]. The objective value is v = 20.

Example 4 . 2

 42 Consider the QFLTL formula φ = F ([A] ≥ x) ∧ F ([A] ≤ y) ∧ x -y ≥ 10. The possible values of v range from -∞ to the maximum value M reached by [A] in the trace. Similarly, w ranges from the minimum m of [A] to ∞. Consequently, the quantity x -y ranges from -∞ to M -m, such that the constraint x -y ≥ 10 constraints the maximal amplitude of variations of [A] in the trace.

Example 4 . 3

 43 Consider again the formula φ of example 4.1:

Example 4 . 4

 44 Consider the formulae of example 4.2

Figure 1 .

 1 Figure 1. Dynamical behavior of the cell cycle model. The plots represent total cyclin (YT) and maturation promoting factor (MPF). (a) Oscillatory behavior obtained with parameter values k Tyson . (b) Higher MPF peaks obtained with k * T yson (solution of problem S1). (c) Shorter oscillations period obtained with k * 4 (solution of problem S4).

Figure 2 .

 2 Figure 2. Violation degree landscape of problem S2. This violation degree measures amplitude of oscillations. Non oscillating regions have highest violation degree.

Figure 3 .

 3 Figure 3. Diagram of the MAPK signal transduction cascade.

  pattern of formulae can be used to search the values of all the 22 parameters of the model to fit the concentration [MEK-RAF ˜{p1}] at six time points. The objective values for these time points are the values of the original model, obtained by simulation with the original parameters k MAPK . The initial values for the search are some random altered values k MAPK alt . Numerical simulations obtained with k MAPK , k MAPK alt and the resulting parameter values are given in Figure 4. It took 290 s to obtain the result.

Figure 4 .

 4 Figure 4. Dynamical behavior of the MAPK model. The curves display [MEK-RAF ˜{p1}]. (a) Reference curve obtained with k MAPK (b) Simulated curve obtained with altered parameter values k MAPK alt . Points are the reference values taken from curve (a). (c) Simulated curve obtained after curve fitting (solution of problem S5).

Figure 5 .

 5 Figure 5. Oscillations of MAPK found with CMA-ES in BIOCHAM

  denotes Mass Action law kinetics with parameter k while ~{p1} and ~{p1,p2} denote phosphorylated forms of a molecule. The rate of reaction 4 is described by:

	MA(k1)	for		_ => Cyclin.
	MA(k3)	for Cyclin + Cdc2~{p1} => Cdc2-Cyclin~{p1,p2}
	MA(k4p)	for Cdc2-Cyclin~{p1,p2} => Cdc2-Cyclin~{p1}
	AUTOCAT(k4) for Cdc2-Cyclin~{p1,p2} => Cdc2-Cyclin~{p1}
	MA(k6)	for Cdc2-Cyclin~{p1}	=> Cyclin~{p1} + Cdc2
	MA(k7)	for	Cyclin~{p1} => _
	MA(k8)	for	Cdc2 => Cdc2~{p1}
	MA(k9)	for	Cdc2~{p1} => Cdc2
	MA(k)

Acknowledgements

This work is partly supported by the EU FP6 STREP project TEMPO and by the INRA Agrobi project INSIGHT. We acknowledge fruitful discussions with our partners of these projects.


Thomas Guérin
email: thomas.guerin@univ-lr.fr

Xavier Bertin
email: xavier.bertin@univ-lr.fr

Eric Chaumillon
email: eric.chaumillon@univ-lr.fr

Wave control on the rhythmic development of a wide estuary mouth sandbank : a process-based modelling study

Keywords: Longe de Boyard sandbank, morphodynamic modelling, stratigraphy, wave action, multi-sized sediment, seismic reflector

This study presents the numerical hindcast of the morphological changes and stratigraphy of a wide estuary mouth sandbank located along the French Atlantic coast (Marennes-Oléron Bay), and primarily investigates the respective contribution of tide and waves to its pluri-decadal evolution. Firstly, the morphodynamic hindcast appears strongly improved when tide and waves are considered in the forcing rather than only tide. Secondly, the modelling results suggest that the strong seasonal variation of the wave climate in this area explains the rhythmic lateral accretion of the bank (i.e. normal to the currents main direction) observed in seismic reflection profiles. The grain size variation in the area of sediment accretion is also characterized by a seasonal cycle. Finally, repetitive bathymetric surveys and seismic profiles revealed that the lateral accretion of the bank is associated with about one seismic reflector per year, which is in agreement with our modelling results when suggesting that these seismic reflectors are due to the seasonal variation of the grain size. These new results may be considered for the morphological evolution and the stratigraphy of other tidal sandbanks potentially affected by waves.

Introduction

Sandbanks appear widely in coastal and shelf seas, where sand is abundant enough and hydrodynamic conditions are sufficient to move the sediment material. Understanding the evolution of these bed forms is important for coastal activities such as navigation and dredging, but also to quantify sediment exchange with adjacent shorelines. Due to the different classifications of sandbanks and ridges in the literature (e.g. [START_REF] Pattiaratchi | Mechanisms for linear sandbank formation and maintenance in relation to dynamical oceanographic observations[END_REF]Swift et al., 1991;Dyer and Huntley, 1999, hereinafter DH99), some confusion has been observed in their terminology. In order to unify the approaches of marine geologists and physical oceanographers, DH99 proposed a classification of sandbanks and ridges which will be used in the present study.

Though the formation and maintenance of tide-dominated sandbanks like open shelf or headland associated sandbanks (type 1 and 3 in DH99, respectively) are well explained theoretically (e.g. [START_REF] Huthnance | On one mechanism forming linear sand banks[END_REF]Signell and Harris, 2000), estuary mouth sandbanks (type 2 in DH99), which include sandbanks emplaced in wide estuaries and in tidal inlets (ebb and flood deltas), have for a long time received less attention (DH99). More precisely, while the emergence of morphodynamic modelling provided a substantial improvement in the understanding of the processes responsible for the development and morphological changes of tidal inlets (e.g. [START_REF] Cayocca | Long-term morphological modeling of a tidal inlet : the Arcachon Basin, France[END_REF][START_REF] Dastgheib | Long-term process-based morphological modeling of the Marsdiep Tidal Basin[END_REF][START_REF] Dissanayake | Effect of sea level rise in tidal inlet evolution : A numerical modelling approach[END_REF][START_REF] Nahon | Process-based 2DH morphodynamic modeling of tidal inlets : a comparison with empirical classifications and theories[END_REF], wide estuary mouth sandbanks have received very little attention, potentially because they are less frequent than ebb and flood deltas. Yet, these sandbanks are very interesting because their development in shallow waters and their exposure to waves could suggest that they may not be only tidally controlled.

Numerous studies have provided qualitative information on this type of sandbanks (e.g. [START_REF] Ludwick | Tidal currents and zig-zag sand shoals in a wide estuary entrance[END_REF][START_REF] Harris | Large-scale bedforms as indicators of mutually evasive sand transport and the sequential infilling of wide-mouthed estuaries[END_REF][START_REF] Gómez | Largo Bank : A shoreface-connected linear shoal at the Bahía Blanca Estuary entrance, Argentina[END_REF][START_REF] Chaumillon | Évolution temporelle et architecture interne d'un banc sableux estuarien : la Longe de Boyard (littoral atlantique, France)[END_REF][START_REF] Kapsimalis | Formation and growth of an estuarine sandbank : Saint-Georges Bank, Gironde Estuary (France)[END_REF] but, to our knowledge, no study has managed to quantify the respective contribution of tides and waves to the long-term morphodynamic of these sandbanks. In order to do so, process-based morphodynamic modelling systems appear nowadays as attractive tools since they allow to perform coastal morphodynamic simulations by coupling currents, waves, sediment transport, and bottom change. Moreover, the recent development of multi-class and multi-layer (MCML) methods in sediment transport models allows accounting for the granulometric heterogeneity of the sediment, and even computing numerical stratigraphy. One can cite the work of [START_REF] Geleynse | Modeling of a mixed-load fluvio-deltaic system[END_REF], [START_REF] Geleynse | Controls on river delta formation ; insights from numerical modelling[END_REF] and [START_REF] Viparelli | Comparison between experimental and numerical stratigraphy emplaced by a prograding delta[END_REF] who have successfully simulated the formation of river deltas and their stratigraphy. Nevertheless, such promising coupled approaches still have to be applied to other coastal environments, such as wide estuary mouth sandbanks.

In this study, a horizontally two-dimensional (2DH), unstructured grid, morphodynamicstratigraphic modelling system was used to investigate the respective contribution of tides and waves to the 40-year retrospective evolution (period 1960 to 2000) of a wide estuary mouth sandbank located along the French Atlantic coast in the Marennes-Oléron Bay. To this aim, two morphodynamic simulations starting with the bathymetry of 1960 were considered, the first including tidal forcing only and the second including tidal and waves forcing. The bathymetry of 2000 together with seismic data were then used to assess the morphological and stratigraphic numerical results.

In the next section, the study area is presented along with a summary of previous knowledge available for this sandbank. The morphodynamic modelling system is described in the third section, with emphasis on the MCML method which has been implemented into the sediment transport module. Finally, the confrontation between numerical results and field data are presented and discussed.

Study area

Geomorphic setting

The Longe de Boyard Sandbank is a wide estuary mouth sandbank located along the French Atlantic coast, in the northern entrance of the Marennes-Oléron Bay (Fig. 1). It is about 8 km long and 2 km wide with a tide-dominated morphology (flood lobe), and its bed elevation ranges from about -17 m from mean sea level (MSL) near the flanks up to -3 m MSL at the southern end of the sandbank. The internal architecture of this sandbank was revealed by seismic profiling and three main seismic units were identified : an upper one characterized by fine and medium sand, lying on an old core of coarse sand and gravel, itself lying on a transgressive valley-fill [START_REF] Chaumillon | Évolution temporelle et architecture interne d'un banc sableux estuarien : la Longe de Boyard (littoral atlantique, France)[END_REF][START_REF] Chaumillon | Multi time-scale evolution of a wide estuary linear sandbank, the Longe de Boyard, on the French Atlantic coast[END_REF].

Hydrodynamic setting

This sandbank is subjected to a semi-diurnal, macrotidal regime with tidal range going from less than 2 m during neap tides up to more than 6 m during spring tides. The flood tidal prism is about 10% more than the ebb tidal prism at the northern entrance of the bay due to the connection of the Marennes-Oléron Bay with the ocean at both ends of the Oléron Island. The main incoming freshwater discharge is due to the Charente River, with a mean annual discharge of 70 m 3 /s (Toublanc et al., 2015). The dominant offshore wave conditions are characterized by significant wave heights ranging from 1 to 2 m and potentially exceeding 6 m during storms, peak periods ranging from 8 s to 12 s, and waves coming from West to NorthWest. The Bay of Biscay is subjected to strong seasonal and interannual wave climate variations [START_REF] Dodet | Wave climate variability in the North-East Atlantic Ocean over the last six decades[END_REF], with winter periods potentially characterized by very energetic storm waves [START_REF] Idier | Contributions respectives des courants et de la houle dans la mobilité sédimentaire d'une plate-forme interne estuarienne. Exemple : le seuil interinsulaire, au large du pertuis d'Antioche, France[END_REF] episodically exceeding 10 m of significant wave height in deep water [START_REF] Bertin | The contribution of shortwaves in storm surges : Two case studies in the Bay of Biscay[END_REF].

Previous knowledge on the sandbank

From 1824 to 2000, an anticlockwise rotation of the sandbank was observed, with sediment accretion of its North-West (NW) and South-East (SE) parts and erosion of its NE and SW parts [START_REF] Chaumillon | Évolution temporelle et architecture interne d'un banc sableux estuarien : la Longe de Boyard (littoral atlantique, France)[END_REF]. Seismic reflection exploration of the sediment accretion areas, presented in [START_REF] Chaumillon | Évolution temporelle et architecture interne d'un banc sableux estuarien : la Longe de Boyard (littoral atlantique, France)[END_REF][START_REF] Chaumillon | Multi time-scale evolution of a wide estuary linear sandbank, the Longe de Boyard, on the French Atlantic coast[END_REF], showed inclined reflectors dipping channelward and recording the rhythmic lateral accretion on both sides of the sandbank. Bertin and Chaumillon (2005) explained the observed sediment accretion at the SE part of the sandbank by the crestward convergence of the residual sediment transport. They proposed that this convergence may be due to the crestward deviation of currents in this area, in response to the observed gradient of the free surface elevation above the crest of the bank which is caused by the obliqueness of the crest compared to the axis of the bank.

Since this sandbank presents a tide-dominated morphology and is oriented along the dominant tidal currents direction (i.e. SE-NW axis), one may suggest that tidal processes prevail over wave-induced processes. However, [START_REF] Chaumillon | Multi time-scale evolution of a wide estuary linear sandbank, the Longe de Boyard, on the French Atlantic coast[END_REF] proposed that wave action could explain the erosion of the seaward part of the bank over the period 1824 to 2000, because of both the relative weakness of tidal currents and the local amplification of wave heights due to convex refraction in this area. Indeed, the strong decrease in wave height from the central part of the bank to its southern end and adjacent channels would induce flows in the direction of wave propagation and adjacent channels.

Nevertheless, these previous studies left us with several open scientific questions of first importance : (i ) Can the respective contribution of tide and waves to the sandbank long-term evolution be quantified ? (ii ) How to explain that lateral accretion occurs, perpendicularly to the main current direction ? (iii )

Where does the rhythmic behavior of the sediment accretion come from ? These questions, which could also be relevant for other wave-influenced sandbanks and sandbanks characterized by lateral accretion, are investigated in the present work.

3. The morphodynamic modelling system 3.1. Model description

Hydrodynamic module

The core of the modelling system used in this study is the unstructured-grid, parallelized, Semi Eulerian-Lagrangian Finite Element (SELFE) hydrodynamic model of [START_REF] Zhang | SELFE : A semi-implicit Eulerian-Lagrangian finite-element model for cross-scale ocean circulation[END_REF]. The main feature of this model is the use of the Eulerian-Lagrangian method to treat the advection in the momentum equations in addition to semi-implicit schemes, which relaxes the numerical stability constraints of the model (e.g. the Courant-Friedrichs-Lewy (CFL) condition can be exceeded). The 2DH configuration of SELFE is used in this study, i.e. the Saint-Venant equations are solved at each hydrodynamic time step to provide the evolution of the surface elevation and the depth-averaged currents over the all computational grid.

Wave module

The Wind Wave Model II (WWM-II) of Roland et al. (2012) (third generation, spectral wave model) is coupled to SELFE and simulates gravity waves generation and propagation by solving the wave action equation (WAE) [START_REF] Komen | Dynamics and Modelling of Ocean Waves[END_REF]. WWM-II uses a residual distribution scheme [START_REF] Abgrall | Residual distribution schemes : Current status and future trends[END_REF] to solve the geographical advection in the WAE, which also relaxes CFL constraints and allows using large time step without compromising the numerical stability.

Regarding the spectral advection (i.e. in the frequency and direction space) in the WAE, it is solved with an Ultimate Quickest scheme [START_REF] Leonard | The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection[END_REF].

Sediment transport and bed evolution module

The 2DH sediment transport/bottom evolution module SED2D [START_REF] Dodet | Morphodynamic modelling of a wave-dominated tidal inlet : the Albufeira Lagoon[END_REF], adapted from MORSYS2D model [START_REF] Bertin | Simulating morphodynamics with unstructured grids : Description and validation of a modeling system for coastal applications[END_REF], allows computing the sediment transport using classical semi-empirical formulations based on depth-averaged velocity, water depth, bottom roughness, sediment properties and wave parameters.

The choice of the sediment transport formulation used to run the morphodynamic simulations was done by comparing measured transport rates from the SANDPIT project (van Rijn et al., 2005) with transport rates predicted by the semi-empirical formulations of van Rijn et al. (2004) [hereinafter TR04], van Rijn (2007a,b) [VR07], [START_REF] Larson | A unified sediment transport model for inlet application[END_REF] [CL11], and [START_REF] Wu | Nonuniform sediment transport under non-breaking waves and currents[END_REF] [WL14]. Two series of sediment transport measurements were considered : a first one corresponding to shallow conditions and a second one corresponding to deeper conditions (see Tab. 1).

SANDPIT tests 5 to 9 SANDPIT tests 19 to 24 h 3.9 to 5.3 m 13.9 to 14.5 m is the significant wave height, Tp is the wave peak period, Uc is the depth-averaged current velocity, and d10, d50 and d90 are the grain diameters for which 10%, 50% and 90% of the grains by mass is finer, respectively.

H s 0.
Only the suspended-load was measured in the scope of the SANDPIT project but, since the suspended-load transport (q s ) is the main contribution of the total transport in this type of coastal conditions, these datasets allowed us to perform relevant comparisons between the four transport formulations.

The scatter plots related to these comparisons are shown on Fig. 2. It can be seen on this figure that q s is globally better predicted with the formulations of TR04 and WL14 than with VR07 and CL11 formulations. The scatter index (SI) values are indeed much higher for VR07 and CL11 formulations, because of the significant overestimation of q s for the deeper conditions (blue circles) which contributes more to increase the SI than underestimating q s . Moreover, although the formulations of TR04 and WL14 perform well in terms of SI score (both SI equal about 60%), the TR04 formulation tends more to overestimate q s than the WL14 formulation. Morphodynamic simulations were run with these four formulations and, consistently with the previous comparisons of predicted transport rates, the bathymetry appeared to degenerate significantly with the TR04, VR07 and CL11 formulations as attested by strongly negative Brier Skill Scores (as defined in Fortunato et al. ( 2014)) ranging from -7.51 to -1.08 (not shown). Therefore, we chose to use only the WL14 formulation thereafter. The formulation of WL14 was developed to estimate the transport of a multisized sediment under the combined action of (non-breaking) waves and currents.

For each grain size class i, the bed-load transport q b,i = (q b,i,x , q b,i,y ) is compu-ted as :

q b,i = q b,i,on + q b,i,of f (1)
where q b,i,on and q b,i,of f are respectively the onshore and offshore part of the bed-load transport during the wave period T w . The respective direction of q b,i,on and q b,i,of f is computed according to the sketch of [START_REF] Wu | Nonuniform sediment transport under non-breaking waves and currents[END_REF], and the norm of these two vectors is computed as :

q b,i,on(of f ) = 5.3.10 -3 f i (s -1)gd 3 i T on(of f ) T w τ b,on(of f ) τ cr,i -1 2.2 (2)
where f i is the local available fraction of sediment class i, s is the ratio between sediment density ρ s and water density ρ, d i is the grain size diameter for sediment class i, and T on(of f ) is the onshore (offshore) half wave period.

The onshore and offshore grain shear stresses are computed as :

     τ b,on = 1 2 ρf cw (U 2 c + U 2 wm,on + 2U c U wm,on cos(φ)) τ b,of f = 1 2 ρf cw (U 2 c + U 2 wm,of f + 2U c U wm,of f cos(π -φ)) (3)
where f cw is the grain friction coefficient for combined current and waves, U c is the depth-averaged current velocity, U wm,on(of f ) is the root-mean-square value of the wave orbital velocity over T on(of f ) , and φ is the angle between current and wave directions.

τ cr,i is the critical shear stress for the incipient motion of sediment class i, which is computed as :

τ cr,i = g(ρ s -ρ)d i θ cr p e,i p h,i -0.6 (4)
where θ cr = 0.03 is the critical Shields number for incipient motion of the median diameter (d50) of the bed sediment which was calibrated by [START_REF] Wu | Nonuniform sediment transport under non-breaking waves and currents[END_REF] using laboratory and field data. p e,i and p h,i are respectively the exposed and hidden probabilities of grain size d i in the bed material, computed following [START_REF] Wu | Nonuniform sediment transport in alluvial rivers[END_REF] :

p e,i = N j=1 f j d i d i + d j and p h,i = N j=1 f j d j d i + d j (5)
where N is the number of grain-size classes, and f j is the fraction of sediment size class j in the bed material.

In the absence of waves, the norm of the bed-load transport is computed as :

q b,i = 5.3.10 -3 f i (s -1)gd 3 i τ b,c τ cr,i -1 2.2 (6)
where τ b,c is the grain shear stress due to currents which is computed as :

τ b,c = 1 2 ρf c U 2 c (7)
where f c is the grain friction coefficient due to currents. In this case, the direction of q b,i is the current direction.

The suspended-load transport q s,i is assumed to be along the current direction, and its norm is computed as :

q s,i = 2.62.10 -5 f i (s -1)gd 3 i τ b τ cr,i -1 U c w s,i 1.74 (8)
where w s,i is the settling velocity of sediment class i. The total bed shear stress τ b is computed as :

τ b = τ 2 b,c + τ 2 b,wm + 2τ b,c τ b,wm cos(φ) (9)
where τ b,c = 1 2 ρf c U 2 c and τ b,wm = 1 4 ρf w U 2 w are the bed shear stresses due to current only and due to waves only, with f c and f w being the related friction coefficients, and U w being the amplitude of the wave orbital velocity near the bed.

Finally, the bed evolution is computed at each grid node in SED2D by solving the Exner equation :

∂z b ∂t + 1 1 -λ ∇ x • q = 0 ( 10
)
where x = (x, y), z b (x, t) is the bed level elevation (positive upwards), λ is the sediment porosity, and q = (q x , q y ) = N i=1 (q b,i + q s,i ) is the depth-integrated sediment total transport rate (in m 3 /m/s) computed at element centres. The finite volume WENO-based method of [START_REF] Guérin | A numerical scheme for coastal morphodynamic modelling on unstructured grids[END_REF] is used to solve this equation while preventing the development of numerical oscillations.

Model implementation

Computational grid and external forcings

The unstructured computational grid used in this study extends from the south of the Maumusson inlet to the north of the Pertuis Breton, and the water depth reaches about 60 m at the offshore boundary (Fig. 1 (B)). Its resolution ranges from 2.5 km offshore to 100 m in the sandbank area. The wave forcing is obtained from the hindcast of [START_REF] Bertin | A significant increase in wave height in the North Atlantic Ocean over the 20th century[END_REF] which used the spectral wave model WaveWatch-III (Tolman and others, 2009) forced with NCEP reanalysis wind fields [START_REF] Kalnay | The NCEP/NCAR 40-year reanalysis project[END_REF].

A morphological factor of 10 was used in order to speed up the simulation, meaning that the morphological changes computed during the morphological time step (i.e. the time between two bathymetric updates) was multiplied by 10 while the total simulation time was 4 years instead of 40. To be consistent with the use of a morphological factor, a representative tidal forcing was considered as proposed in the work of [START_REF] Latteux | Techniques for long-term morphological simulation under tidal action[END_REF]. By taking as a reference the residual sediment transport computed over the sandbank area for a one-year simulation with realistic tidal and wave forcings, the representative tide was adjusted such as the corresponding representative simulation gives a residual sediment transport close to the reference result over the whole sandbank area.

The representative tide selected is composed of the quarter-diurnal constituent (M4) and a M2 * constituent which has the same period than the semi-diurnal constituent M2 (T M 2 = 12.42 h) but an amplitude 30% higher than the M2 amplitude. Consequently to the use of a morphological factor, the time corresponding to the realistic wave forcing was squeezed (i.e. 10 minutes became 1 minute).

As explained in [START_REF] Zhang | SELFE : A semi-implicit Eulerian-Lagrangian finite-element model for cross-scale ocean circulation[END_REF], the Eulerian-Lagrangian method employed in SELFE constrains the use of relatively large time steps (i.e. for which the CFL number is close to 1 or even superior to 1) for better accuracy.

In our case, a stable configuration was obtained with a hydrodynamic time step of 120 s, giving a CFL number less than 2 (maximum current velocity in the fine resolution (100 m) area is about 1.5 m/s leading to a CFL number of 1.5*120/100 = 1.8). The morphological time step was set to the period of the main tidal constituent (T M 2 = 12.42 h) in order to speed up the simulation, while remaining sufficiently small compared to the relatively slow morphological changes.

Since both the grid resolution and the modelling system were not configured to handle the diverse physical processes taking place all over the computational grid (e.g. the high fraction of cohesive sediments in the channels along the sandbank sides and in the southern part of the bay, and the wave induced dynamics at the surrounding shoreline which would require a finer grid resolution), a morphological ramp was used in order to update the bathymetry only in the sandbank area and to avoid unrealistic morphological changes elsewhere. To verify that this ramp doesn't impact the morphodynamic hindcast, we performed a simulation where the surrounding bathymetry of 1960 was replaced by the one of 2000. This simulation didn't yield significantly different morphological changes of the sandbank, which would suggest that the use of such a ramp is a reasonable simplification. A negligible impact on the morphodynamic hindcast was also obtained concerning the sea level rise that took place between 1960 to 2000 (about +8.4 cm according to [START_REF] Gouriou | Reconstruction of a twocentury long sea level record for the Pertuis d'Antioche (France)[END_REF]).

The multi-class and multi-layer (MCML) method

Considering a multiple-sized sediment led us to implement a MCML method into our sediment transport and bed evolution module SED2D. The method proposed by [START_REF] Reniers | Observations and modeling of steep-beach grain-size variability[END_REF] was used in this study, allowing us to obtain the temporal evolution of the sediment composition 3D distribution (grain-size and sorting). The initial configuration consisted of 400 layers of 1 cm thickness, along with three grain-size classes of 0.1 mm, 0.2 mm, and 0.35 mm diameter and respective fractions of 10%, 80%, and 10% according to the vibrocore data presented in [START_REF] Chaumillon | Multi time-scale evolution of a wide estuary linear sandbank, the Longe de Boyard, on the French Atlantic coast[END_REF]. At every time step, the median diameter was computed in each layer using a geometric mean (initial value of nearly 0.2 mm).

Modelling results

Water levels

Bottom pressure measurements were performed with a Nortek Vector current meter equipped with a pressure sensor at a station located at the NW side of the sandbank (see station La Brée on Fig. 1 4 (a). Regarding these two stations the related root-mean-square error (RMSE) ranges between 0.14 and 0.21 m for the total water depth (i.e. sum of tide and surge), which yields an error of about 4 to 5% when normalized by the mean tidal range. Such result is satisfactory when compared to previous modelling studies on this area [START_REF] Bertin | Importance of wave age and resonance in storm surges : The case Xynthia, Bay of Biscay[END_REF]Toublanc et al., 2015).

Waves

The computation of wave parameters from the bottom pressure measurements recorded at La Brée station was done by first correcting these measurements from sea level pressure measured at the nearby station of Chassiron.

A bottom pressure spectra was then computed using Fast Fourier Transform over 20 min bursts, and this pressure spectra was converted into elevation spectra using linear wave theory. The significant wave height H m0 was computed as 4

√ m 0 where m 0 is the 0 th moment of the elevation spectra, and the peak period T p as the period where the maximum energy was found in the range Regarding the Antioche station, the instrument was not calibrated to extract wave parameters.

Currents

Vertical profiles of current velocity were measured at the Antioche station (see Fig. 1 (B)) from January to February 2010, and were also measured several times along a same transect located at the south end of the sandbank (see transect A on Fig. 1 (C)) to cover a tidal cycle on November 6 th of 2002.

These two datasets were first used to estimate the local Manning's number (n) which is needed as an input parameter in our hydrodynamic module to handle the bottom friction. According to the method of Soulsby (1997), a linear regression was performed on each scatter plot {U (z); ln(z)}, where z is the vertical distance from the bottom and U (z) is the related horizontal fluid velocity. This linear regression allowed us to compute the friction velocity u * and the bed roughness length z 0 by using the expression U (z) = Aln(z) + B where A = u * /K and B = -Aln(z 0 ), K being the von Karman's constant. The drag coefficient was then computed as C d = [K/(1 + ln(z 0 /h))] 2 , where h is the water depth, and the Manning's number was finally computed as n = C d h 1/3 /g. Since this method was developed for current-only conditions, periods of calm weather (i.e. weak waves and wind) were favored and velocity measurements in the sub-surface layer were removed for this computation. Though this analysis was not sufficient to build a spatial distribution of the Manning's number over the all computational grid, it allowed us to estimate a mean value of the Manning's number which was then used in all our numerical simulations. This value appeared equal to 0.03 and was used for all our numerical simulations.

The depth-averaged current velocity (U c ) obtained from measurements was compared with the simulated one at Antioche station (see Fig. 4 (b)) and along the transect A (not shown). While the model underestimates the current velocity at Antioche station (mean bias of -0.06 m/s), this trend disappears for the results along the transect A (mean bias ranging between about -0.05 and 0.05 m/s). The related normalized root-mean-square error (NRMSE) at the Antioche station is about 22%, while for the transect A it reaches about 15% when U c > 0.5 m/s (transects made during mid-tide) and about 30 to 40% when U c < 0.5 m/s (transects made close to high tide and low tide). Regarding the current direction, the bias between simulated and measured values at the Antioche station is almost constant (about 20 • , see Fig. 4 (c)), but it becomes more variable along the transect A (ranging between about -15 • and 15 • ).

Morphological evolution

Bathymetric survey of the bank was done in 1960 and 2000 (Fig. 1 (C) and (D)) with bathymetric uncertainty of 0.35 m (Bertin et al., 2005), allowing us to study the 40-year morphological evolution of this sandbank. The measured morphological evolution (Fig. 5 (a)) is characterized by lateral sediment accretion of the SE and NW flanks of the bank, and sediment erosion mainly at the south of the bank and along its SW flank. Considering that both tide and wave-related processes may explain the observed evolution of the sandbank, two morphody-

Discussion

Limitations of our modelling approach

The comparison between model results and the available field data revealed firstly that water levels, wave parameters and depth-averaged currents are reproduced with an accuracy corresponding to the state of the art for such an estuarine environment. However, the 2DH configuration of our modelling system is a first limitation since 3D circulations may lead to a possible veering of the currents along the vertical direction, thus influencing the sediment transport. Then, more current velocity measurements would allow to enhance the parametrization of the bottom friction, for instance it would be relevant to use a realistic spatial distribution of the Manning's number over the all computational grid. Regarding the sediment composition, though the integration of a multiple-sized sediment through the MCML method is an important improvement compared to a single-sized sediment, simulations would still benefit from a spatially-varying d50 initialization. Sediment availability could also be improved since the assumption of an unlimited quantity of sand over the full computational grid is considered in our study. An additional simplification is the use of a morphological factor, which prevents us to analyze the impact of the tidal range on the morphological evolution of the bank. Finally, the model misses the role of cohesive sediments, mainly present in the channels adjacent to the sandbank and further south between the sandbank and the Charente Estuary (Fig. 1).

Cohesive sediments can significantly impact the hydrodynamics, for instance the bottom shear stress may decrease compared to the case of non-cohesive sediment. The sediment transport can also be modified in several ways due to the presence of mud, including the process of dynamic bed armouring, where the resuspension of a sediment mixture (sand+mud) is limited or prevented by the present mud fraction (Le [START_REF] Hir | Dynamics of sand and mud mixtures : a multiprocess-based modelling strategy[END_REF].

Despite these limitations, simulating the evolution of a sandbank subjected to tides and waves over 40 years remains very challenging and we propose that our modelling results still provide some relevant leads about the processes that control the long-term evolution of this wide estuary mouth sandbank.

Driving mechanisms for the long-term evolution of the sandbank

The first result of the present study is the significant improvement of the modelled morphological evolution obtained when a realistic wave forcing was added to the tidal forcing. This result confirms the assumption of Chaumillon Since the bed evolution is proportional to the divergence of sediment transport, computing the divergence of the residual sediment transport shows where erosion and accretion occur over a tidal cycle (Fig. 6 : colors). Indeed, at the southern tip of the bank, the decreasing water depths towards the bank crest induces an increase of the residual sediment transport (i.e. minus their divergence is negative ; blue color on Fig. 6), and therefore an erosion of the corresponding area. Whereas along the SE flank of the bank, both the convergence of the residual transport and the increasing water depths (thereby decreasing sediment transport) towards the channel induce an accretion of the area (i.e. minus the divergence of residual transport is positive ; red color on Fig. 6). Interestingly, the erosion and accretion patterns computed over a single tidal cycle match the long term evolution of the bank, which would suggest that, in this area, this long term evolution is well explained by the sediment transport in the presence of energetic waves and averaged over a tidal cycle.

The sediment transport direction appeared also affected by wave action, especially along the SE flank of the sandbank where the main accretion is observed. The channelward component of residual sediment transport in this area appeared clearly visible when waves were added, explaining the lateral sediment accretion observed along this flank. Therefore, the fact that sediment deposition occurred in a direction normal to the dominant tidal currents direction (i.e.

the sandbank main axis), which could be thought rather counterintuitive, is explained by the presence of waves. This result is of great importance as lateral accretion of sandbanks is widely documented (e.g. [START_REF] Houbolt | Recent sediments in the southern bight of the North Sea[END_REF][START_REF] Harris | Large-scale bedforms as indicators of mutually evasive sand transport and the sequential infilling of wide-mouthed estuaries[END_REF]; Berné [START_REF] Berné | Architecture and long term evolution of a tidal sandbank : The Middelkerke Bank (southern North Sea)[END_REF][START_REF] Dalrymple | Estuarine dunes and bars[END_REF][START_REF] Olariu | Anatomy of a laterally migrating tidal bar in front of a delta system : Esdolomada Member, Roda Formation, Tremp-Graus Basin, Spain[END_REF][START_REF] Chaumillon | Tidal and fluvial controls on the internal architecture and sedimentary facies of a lobate estuarine tidal bar (The Plassac Tidal Bar in the Gironde Estuary, France)[END_REF].

These large differences of residual transport patterns between the tide-only and the tide and waves cases can be explained by comparing the instantaneous sediment transport for current alone and for combined current and wave over a wave period (Fig. 7). Though adding the wave orbital velocity to the current velocity does not modify the mean total velocity compared to the case with current alone (i.e. in both cases the mean total velocity equals to the current velocity), the corresponding mean sediment transport for current and wave is about one order of magnitude higher than the one for current alone. This difference is due to the non-linear dependence between the sediment transport and the fluid velocity (i.e. the sum of the current velocity and the wave orbital velocity), and explains the increase of sediment transport rates shown on Fig. 6.

Origin of the rhythmic accretion

Since the wave control on the sandbank evolution was found to be essential, the next step of this study was to analyze how the strong seasonal variation of Moreover, the integration of a MCML approach into SED2D enabled us to study the impact of the wave climate variation on the sediment composition.

Conclusions

This process-based modelling study highlighted the strong control of waves on the morphological evolution of a wide estuary mouth sandbank, in spite of its global tide-dominated morphology. The divergence of residual sediment transport over one single tidal cycle and in the presence of energetic waves is consistent with the long-term evolution of the bank, which would suggest that this long-term evolution is well explained by the short-term sediment dynamics during energetic periods. It was shown that these energetic waves, mostly occuring during winter periods, induce an important channelward sediment transport in the main accretion area of the bank, which explains the observed lateral accretion. This rhythmic accretion consists of an alternation of thicker sand layers, emplaced during periods of higher wave energy, and thinner sand layers, emplaced during periods of lower wave energy. A seasonal cycle is also found in both the temporal and vertical variation of the sediment granulometry. Based on the fact that seismic reflectors can be related to grain size variations, a frequency of about one reflector per year was obtained with the model, which appears supported by recent seismic data. These results may be considered for other estuary mouth sandbanks and other sandbanks that are potentially affected by waves.

Figure 1 :

 1 Figure 1: (A) Location of the study area in the Bay of Biscay. (B) Bathymetric map of the study area with respect to mean sea level, with location of the Longe de Boyard sandbank (black square) in the Marennes-Oléron Bay, near the Charente Estuary. (C) and (D) Bathymetry of the sandbank in 1960 and 2000 respectively, with location of the seismic profile realized in 2014, the transect A corresponding to current velocity measurements, the wave station MO1, and the numerical core station MO2. Coordinates of (B), (C) and (D) are in meters (Lambert-93).

Figure 2 :

 2 Figure 2: Comparison between suspended-load transport rates measured in the scope of the SANDPIT project with predicted values given by the formulations of : (a) van Rijn et al. (2004) [TR04], (b) van Rijn (2007a,b) [VR07], (c) Larson et al. (2011) [CL11], and (d) Wu and Lin (2014) [WL14].

  (B)) from November to December 2009, and with a Nortek AWAC current meter at a second station located in the "Pertuis d'Antioche" (see Antioche station on Fig. 1 (B)) from January to February 2010. The comparison between measured and simulated water levels at these two stations are shown on Fig. 3 (a) and Fig.

[ 0 .

 0 33 -0.04Hz]. Simulated H m0 and T p without bathymetric update (i.e. only SELFE-WWM-II, with the bathymetry of 2000) were compared against the values obtained from measurements (Fig.3 (b)-(c)). RMSE and normalized RMSE are respectively 0.21 m and 26.6% for H s , and 1.27 s and 9.6% for T p , which are reasonably good results given the very coarse spatio-temporal resolution of the atmospheric forcing (about 1.9 • in latitude and longitude, and 6h in time).

Figure 3 :

 3 Figure 3: Comparison of measurements and simulated results for water level, significant wave height (Hs), and peak period (Tp) at La Brée station, from November 18 th to December 12 th 2009.

Figure 4 :

 4 Figure 4: Comparison of measurements and simulated results for water level, and depthaveraged current amplitude and direction at the Antioche station, from January 23 rd to February 8 th 2010.

Figure 5 :

 5 Figure 5: Differential bathymetric maps for the 40-year period (1960 to 2000). (a) Measured map ; (b) Simulated map with tidal forcing only ; (c) Simulated map with tidal and waves forcing.

  et al. (2008) according to which waves might have a significant impact on this sandbank evolution, and more generally on estuary mouth sandbanks, though this type of bank can be characterized by an overall tide-dominated morphology and classified as a tidal sandbank. To further investigate the impact of waves on sediment transport, residual sediment transport maps of the southern part of the bank were computed (Fig.6: vectors), for tidal forcing only and for tide and energetic offshore wave conditions (H s = 4 m, T p = 11 s, Dir p = N 290 • , such wave conditions leading to H s reaching about 1.3 m in the sandbank area).

Figure 6 :

 6 Figure 6: Residual sediment transport (arrows) with related accretion and erosion areas (colors) (computed as minus the divergence of the residual transport), for a tidal cycle and in the south area of the bank : (a) tide-only forcing, (b) tide and waves forcing (offshore wave parameters : Hs = 4 m, Tp = 11 s, Dirp = N290 • ).

Figure 7 :

 7 Figure 7: (a) Instantaneous fluid velocities for current alone (1 m/s), and for current with additional sinusoidal wave (wave height of 0.5 m, period of 12 s, wave direction equals to the current direction, and water depth of 4 m). (b) Related instantaneous sediment transport (same grain-size parametrization than in the present modelling study), with mean transport for the combined current and wave case.

Fig. 8

 8 Fig. 8 (b) compares H s (t) at station MO1 with the absolute rate of change of the d50 computed in the surface layer at station MO2. A correlation coefficient of 0.32 with a 95% confidence interval is obtained when considering monthly averaged signals over the 40-year period, which reflects the tendency of the grain size to vary preferentially during periods of energetic waves (i.e. winter periods).

Figure 8 :Figure 9 :

 89 Figure 8: Temporal variation of significant wave height (Hs) at station MO1 (see Fig. 1 (C)) with absolute rate of change of (a) bed level and (b) d50 of the surface layer at station MO2, between 1970 and 1990.

Figure 10 :

 10 Figure 10: Seismic profile in the SE accretion area of the sandbank recorded in July 2014 (see Fig. 1 (C)), with superposed 1960 and 2000 bathymetric profiles.

Acknowledgements

The AWAC data collected at the Antioche station and the Vector data collected at La Brée station were provided by the French Oceanographic and Hydrographic Institute (SHOM), acquired during the projects MOUTON (funded by DGA PEA 012401) and EPIGRAM (funded by LEFE/IDAO and ANR, agreement ANR-08-BLAN-0330-01).

By considering a "numerical core" in the SE accretion area (see Fig. 9), we can also analyze how the vertical variation of the sediment granulometry is recorded through time. Firstly, it can be observed that the largest accretion occurs during winter periods (more dark blue bars than yellow or light blue bars), i.e. periods of energetic waves, which is in agreement with the important wave control on the morphological evolution previously highlighted (Fig. 8 (a)).

Secondly, a seasonal variation of the median grain size is obtained, which is also in agreement with the previous results shown on Fig. 8 (b). It is finally interesting to see that winter periods can be related to peaks of relatively coarse sediment (like the winters of 1983-84 and 1984-85 for instance), while they can also be related to peaks of relatively fine sediment (like the winter of 1987-88 for instance). This complex behavior reflects probably the complex variation, both in time and space, of the local grain size fractions.

A seismic reflection survey of the sandbank realized in July 2014 allowed us to compare seismic and modelled stratigraphy in the SE accretion area. Inclined and channelward dipping reflectors are observed on the seismic profile (Fig. 10), showing contrast of acoustic impedance that can be interpreted as rhythmic lateral migration of the sandbank. Interestingly, the number of reflectors observed on this seismic profile between the bathymetric profiles of year 2000 and 2014 is of the order of 15, which corresponds to approximately one reflector per year.

The number of reflectors between 1960 and 2000 is less easy to count, but a reflector frequency ranging between 0.75 and 1.25 reflector per year is observed on the same figure. Seismic reflectors observed in this area were already correlated with vertical grain size variations [START_REF] Billeaud | Evidence of a major environmental change recorded in a macrotidal bay (Marennes-Oléron Bay, France) by correlation between VHR seismic profiles and cores[END_REF]. Therefore, the vertical distribution of the d50 obtained numerically (Fig. 9) reinforces this assertion, since this variation is characterized by a seasonal cycle (i.e. a frequency of one reflector per year). This result also confirms the previous interpretations of [START_REF] Berné | Architecture and long term evolution of a tidal sandbank : The Middelkerke Bank (southern North Sea)[END_REF] and [START_REF] Chaumillon | Multi time-scale evolution of a wide estuary linear sandbank, the Longe de Boyard, on the French Atlantic coast[END_REF] who proposed that reflectors recording the lateral accretion of sandbanks could result from wave action.

Note that such numerical investigations would not have been possible if an input reduction technique was considered for the wave forcing (e.g. a randomly van Rijn, L., Soulsby, R., Hoekstra, P., Davies, A., 2005. SANDPIT, Sand Transport and Morphology of Offshore Mining Pits. Aqua Publications.

van Rijn, L., Walstra, D., Ormondt, M.v., 2004