Céline Blanc

Laurent Rouvreau

Mathilde Scamps

Excavated soil reuse -Development of a French management framework -Specific or generic approach?

Excavated soil reuse Development of a French management framework

Specific or generic approach ?

Céline BLANC, Laurent ROUVREAU, Mathilde SCAMPS BRGM, 3, avenue Claude Guillemin, BP 36009, 45060 Orléans Cedex 2 -France After the revision of the Waste Framework Directive and its transposition in the French legal framework, it was necessary to define rules and methods to guide the sustainable reuse of soil ensuring human health and environmental protection. One provision of the European Directive sets that excavated soils not reused on site should be classified as waste and managed as such under national policies. To clarify these rules and provide stakeholders a common operational framework, the French Ministry for Environment asked BRGM and INERIS to draft a methodological guidance document related to the off-site reuse of excavated soil (for road construction, or development projects) and develop the tools necessary for implementation at local / regional levels. After a review of European countries practices, these methodological tools have been developed with the support of working groups composed of a wide range of stakeholders (professional associations, environmental protection association, planners, lawyers ...). They are now public available on the French Portal for Contaminated land management (http://www.developpement-durable.gouv.fr/Guide-de-reutilisationhors-site.html)

At the crossing of two management policies

The French policy on contaminated land focuses on two main concepts:

 Risk analysis and management rather than consideration of an intrinsic level of pollution ;

 Management based on the use of the site.

The presentation focuses on the difficulty to initiate a French management policy on excavated soils in agreement with the French policy on contaminated land (site specific approach) and the Waste policy based on a generic approach with guidelines values.

Excavated soil off-site reuse methodology

The presentation will focus on the main principles of the methodology, especially the conditions of reuse.

After characterization of soils to be reused, the methodology provides two types of reuse on a receiver site :

-For road construction, -As part of a development project (for which a building permit or an EIS is issued) and if the three following criteria are respected:

-Criterion 1: soil quality of the receiving site is maintained; -Criterion 2: water resources and ecosystems are preserved beneath the receiver site, -Criterion 3: Excavated soils are compatible with the future use of the receiver site (only in the context of development projects).

Tools developed as part of the french soil management framework

After presenting shortly the French Soil management framework, its aims and the general scheme, this second part of the presentation will focus on tools developed to support the implementation program. Two main tools were developed in order to support the methodology.

The first tool, HydroTex, concerns groundwater risk assessment. It was developed to verify that the offsite reuse of excavated soil will not affect the groundwater resource quality. The Hydrotex spreadsheet provides a specific result for each off-site receiver site and for each substance of concern (on single substance basis). The tool proposes harmonized entrance parameters for the groundwater modeling model. This assessment determines the maximum soil concentration for the relevant substance to protect the groundwater resource at the target endpoint.

To keep track of excavated soil quality, an information traceability system is implemented through the TERRASS database, an interactive tool for excavated soil banking, available online. These traceability measures allow monitoring the excavated soil volumes from a producer site to a receiver site. This receiver site can be a transit hub / cluster, a treatment center or directly a receiver site (for a development project or road construction). The database will also provide functioning indicators.


Assyr Abdulle
email: assyr.abdulle@epfl.ch

Gilles Vilmart
email: gilles.vilmart@bretagne.ens-cachan.fr

Konstantinos C Zygalakis
email: k.zygalakis@soton.ac.uk

Mean-square A-stable diagonally drift-implicit integrators of weak second order for stiff Itô stochastic differential equations

Keywords: 4), 2009] Stiff SDEs, drift-implicit stochastic methods, mean-square stability AMS subject classification (2010): 65C30, 60H35, 65L20

published or not. The documents may come     Mean-square A-stable diagonally drift-implicit integrators of weak second order for stiff Itô stochastic differential equations

Introduction

Consider the system of Itô stochastic differential equations dX(t) = f (X(t))dt + m r=1 g r (X(t))dW r (t), [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF] where X(t) is a random variable with values in R d , f : R d → R d is the drift term, g r : R d → R d , r = 1, . . . , m are the diffusion terms, and W r (t), r = 1, . . . , m are independent one-dimensional Wiener processes. The drift and diffusion functions are assumed smooth enough, Lipschitz continuous and to satisfy a growth bound in order to ensure a unique (mean-square bounded) solution of (1) [START_REF] Arnold | Stochastic differential equations: theory and applications[END_REF]. Except for some very special cases, (1) cannot be solved analytically and numerical methods are needed. The accuracy of a numerical solution is usually measured in terms of strong error (the rate at which the mean of the error norm decays) and weak error (the rate at which the difference between the mean of a smooth functional of the exact and the numerical solutions decays) [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF]. Besides the strong and weak error, the stability of a numerical integrator is an essential issue for many problems. In the case of numerical methods for ordinary differential equations (ODEs) this is a very well studied problem and one desirable property is the so-called A-stability, especially when dealing with stiff problems [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF]. More precisely, considering the linear test problem dX = λXdt, λ ∈ C, whose solution is stable if and only if lim t→+∞ X(t) = 0 ⇐⇒ λ ∈ S ODE := λ ∈ C; ℜ(λ) < 0 and applying a Runge-Kutta method to it leads to the one step difference equation X n+1 = R(p)X n , where the stability function R(p) is a rational function of p = λh. The numerical method is then stable for this problem if and only if lim n→∞ X n = 0 ⇐⇒ p ∈ S num := p ∈ C; |R(p)| < 1 , and the method is called A-stable if and only if S ODE ⊆ S num . For SDEs, different measures of stability are of interest and in this paper we focus on meansquare stability [START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF] and mean-square A-stability [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF] (the generalization of the A-stability for ODEs to SDEs). One considers the following test problem [START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF][START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF][START_REF] Burrage | Numerical methods for strong solutions of stochastic differential equations: an overview[END_REF][START_REF] Tocino | Mean-square stability of second-order Runge-Kutta methods for stochastic differential equations[END_REF]

dX = λXdt + µXdW (t), X(0) = 1, (2)
in dimensions d = m = 1, with fixed complex scalar parameters λ, µ. The exact solution of (2), given by X(t) = exp((λ +1 2 µ 2 )t + µW (t)), is mean-square stable if and only if lim t→∞ E |X(t)| 2 = 0 and mean-square stability can be characterized as the set of (λ, µ) ∈ C 2 such that ℜ(λ)+ 1 2 |µ| 2 < 0, that will be called S M S SDE [START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF][START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF]. Another measure of stability that will be briefly mentioned in this paper is that of asymptotic stability. In particular, the solution of ( 2) is said to be stochastically asymptotically stable if and only if lim t→∞ |X(t)| = 0, with probability 1. Asymptotic stability can be characterized as the the set of (λ, µ) ∈ C 2 such that ℜ λ -1 2 µ 2 < 0. Applying a numerical method to the test SDE (2) usually yields the following one step difference equation [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF] X

n+1 = R(p, q, ξ n )X n , (3)
where p = λh, q = µ √ h, and ξ n is a random variable. We can then characterize the meansquare stability domain of the method as

lim n→∞ E(|X n | 2 ) = 0 ⇐⇒ S M S num := (p, q) ∈ C 2 ; E(|R(p, q, ξ)| 2 ) < 1 . (4)
A characterization of the numerical asymptotic stability domain can also be derived, assuming R(p, q, ξ) = 0 with probability 1 1 and E((log |R(p, q, ξ)|) 2 ) < ∞, as [13, Lemma 5.1] the set of (p, q) ∈ C 2 such that E(log |R(p, q, ξ)|) < 0. Finally, a numerical integrator is called

• mean-square A-stable if S M S SDE ⊆ S M S num ;
• mean-square L-stable, if it is mean-square A-stable and if the limit

E(|R(p k , q k , ξ)| 2 ) → 0 holds for all sequences (p k , q k ) ∈ S M S SDE with ℜ(p k ) → -∞.
If we restrict (p, q) ∈ R 2 then the domains of mean-square or asymptotic stability are called regions of stability. Mean-square A-stable numerical methods are necessarily drift-implicit and it is shown in [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF] that the stochastic θ-methods, which have strong order 1/2 and weak order 1 for general SDEs of the type (1) are mean square A-stable for θ ≥ 1/2. In [START_REF] Higham | A-stability and stochastic mean-square stability[END_REF] it is shown that A-stable methods which have strong and weak order 1 can be built using the θ-method, with meansquare A-stability achieved for θ ≥ 3/2 (notice that such methods might have large error constants and are usually not used in the deterministic case). We mention that a class of strong order one implicit schemes for stiff SDEs, based on the so-called Balanced method, were recently proposed in [START_REF] Alcock | Stable strong order 1.0 schemes for solving stochastic ordinary differential equations[END_REF] with the aim of achieving large asymptotic stability regions. High order strong methods for SDEs are usually difficult to implement due to the need of computing numerically involved stochastic integrals. In contrast higher order weak methods are easier to simulate as the stochastic integrals can in this case be replaced by discrete random variables. However, constructing mean-square A-stable higher order integrators is a non trivial task. In [START_REF] Komori | Stable ROW-type weak scheme for stochastic differential equations[END_REF] a method of ROW type [START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems[END_REF] of weak second order is proposed for Itô SDEs that is mean-square stable under the assumption of real diffusion coefficients. Recently, a class of singly diagonally drift-implicit Runge-Kutta methods of weak second order was proposed in [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF]. These methods, called SDIRK in the numerical ODE literature [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF]Chap. IV.6] are of interest because they are cheaper to implement than fully drift-implicit methods (see also the discussion in Section 2). However, none of the weak second order diagonally implicit Runge-Kutta methods proposed in [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF] are mean square A-stable. Moreover except for the variation of the θ-Milstein method for Itô SDEs that was proposed in [START_REF] Abdulle | High order weak methods for stochastic differential equations based on modified equations[END_REF], we are not aware of any other weak second order mean-square A-stable integrator.

In this paper we derive a class of singly diagonally drift-implicit integrators of weak second order (indexed by parameter γ), that we call S-SDIRK methods, for multidimensional SDEs with non-commutative noise. These methods have the same computational cost as the methods derived in [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF], but much better stability properties. More precisely, for a particular choice of γ, the mean-square A-stability for general parameters (p, q) ∈ C 2 can be proved. For another choice of γ that leads to a third order method for deterministic problems, for which the mean-square A-stability can be checked numerically. Comparison with the methods derived in [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF] is discussed and numerical experiments on a nonlinear test problem with non-commutative noise corroborate the weak second order of convergence predicted by our analysis. Finally, a new stabilization procedure introduced in this paper also allows to improve the stability of the strong and weak order 1 methods introduced in [START_REF] Higham | A-stability and stochastic mean-square stability[END_REF] based on the θ-method. In particular, mean-square A-stable methods for any value θ ≥ 1/2 are constructed.

Mean-square A-stable diagonally drift-implicit integrators

Instead of considering the general framework of stochastic Runge-Kutta methods [START_REF] Komori | Weak second-order stochastic Runge-Kutta methods for non-commutative stochastic differential equations[END_REF][START_REF] Rößler | Second order Runge-Kutta methods for Itô stochastic differential equations[END_REF] we derive our S-SDIRK methods by stabilizing the simplest Taylor based method of weak second order, namely the Milstein-Talay method [START_REF] Talay | Efficient numerical schemes for the approximation of expectations of functionals of the solution of a SDE and applications[END_REF] following the methodology developed in [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF] for explicit stabilized stochastic methods. Consider

K1 = X 0 + hf (X 0 ), K2 = K1 + √ h m r=1 g r (X 0 )ξ r , K3 = X 0 + K1 2 , X1 = X 0 + h 2 f (X 0 ) + f ( K2 ) + 1 2 m r=1 g r X 0 + m q=1 g q (X 0 )J q,r -g r X 0 - m q=1 g q (X 0 )J q,r + √ h 2 m r=1 g r ( K3 + h 2 m q=1 g q (X 0 )χ q ) + g( K3 - h 2 m q=1 g q (X 0 )χ q ) ξ r , (5)
where

J q,r =    h(ξ r ξ r -1)/2, if q = r, h(ξ q ξ r -χ q )/2, if r < q, h(ξ q ξ r + χ r )/2, if r > q, (6)
and χ l , ξ l , l = 1 . . . m are independent discrete random variables satisfying respectively

P(χ l = ±1) = 1/2, P(ξ l = ± √ 3) = 1/6, P(ξ l = 0) = 2/3. ( 7
)
The method ( 5) is obtained from a second weak order Taylor method [START_REF] Talay | Efficient numerical schemes for the approximation of expectations of functionals of the solution of a SDE and applications[END_REF] replacing stochastic integrals with discrete random increments and derivatives by finite differences. Notice that each step of the above derivative-free scheme only involves five evaluations of the functions g r , r = 1, . . . , m, independently of the dimension m of the Wiener processes, thanks to suitable finite differences involving also noise terms, as first proposed in [START_REF] Rößler | Second order Runge-Kutta methods for Itô stochastic differential equations[END_REF]. A direct proof of the weak second order of the method (5) can easily be established and we refer to [4, Lemma 3.1] for details. However, as it can be seen in Figure 1, this method (5) has a very restricted meansquare stability domain, which results in a stepsize restriction in the case of stiff problems. In order to relax this restriction, one needs to stabilize the method [START_REF] Alcock | Stable strong order 1.0 schemes for solving stochastic ordinary differential equations[END_REF]. There are two alternatives ways for doing this: use stabilized explicit integrators [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF][START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF], or use implicit integrators, and in particular diagonally drift-implicit ones, which, as mentioned above, is the focus of this paper.

Diagonally implicit Runge-Kutta methods In the case of numerical methods for ODEs, diagonally implicit Runge-Kutta methods are integrators of the form

K i = Y 0 + i j=1 a ij hf (K j ), Y 1 = Y 0 + s i=1 b i hf (K i ), ( 8
)
where s is the number of internal stages, and a ij , b i are the coefficients of the method. The advantage of diagonally implicit Runge-Kutta methods over fully implicit Runge-Kutta methods is that one can treat one internal stage after the other (nonlinear systems of size d × d) instead of solving the full nonlinear system of size (d • s) × (d • s). An additional advantage is that the values for the internal stages already computed can be used to find a good starting value for the next implicit stage that needs to be computed [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF]. Moreover, choosing identical diagonal coefficients a ii = γ permits to use at each step a single LU-factorization for all quasi-Newton iterations in all internal stages of the method. In this case the corresponding methods are called singly diagonally implicit Runge-Kutta methods (SDIRK).

New weak order two diagonally implicit A-stable integrators We introduce the following stabilized integrator of weak order two for the integration of (1).

K 1 = X 0 + γhf (K 1 ), K 2 = X 0 + (1 -2γ)hf (K 1 ) + γhf (K 2 ), K * 1 = X 0 + β 1 γhf (K 1 ) + β 2 γhf (K 2 ), K * 2 = X 0 + γhf (K 1 ) + D -1 (K * 1 -X 0 ), K * 3 = K * 1 + β 3 hf (K * 2 ), X 1 = X 0 + h 2 f (K 1 ) + h 2 f K 2 + √ h m r=1 g r (K * 2 )ξ r + 1 2 m r=1 g r K * 2 + m q=1 g q (K * 2 )J q,r -g r K * 2 - m q=1 g q (K * 2 )J q,r + √ h 2 m r=1 g r K * 3 + h 2 m q=1 g(K * 2 )χ q + g r K * 3 - h 2 m q=1 g(K * 2 )χ q ξ r , (9)
where

β 1 = 2-5γ 1-2γ , β 2 = γ 1-2γ , β 3 = 1 2 -2γ
, and ξ r , χ r , J q,r satisfy ( 7) and ( 6) respectively. For γ = 0, we have K * 1 = K * 2 = X 0 and K * 3 = X 0 + (h/2)f (X 0 ) and we recover the explicit Milstein-Talay method (5). For the stage K * 2 , we use D -1 to stabilize K * 1 -X 0 , where D = Iγhf ′ (X 0 ). This stabilization procedure (used in numerical ODEs to stabilize the error estimator of an integrator) is well-known in ODEs and has been introduced by Shampine [11, Sect. IV.8], its use for SDEs is motivated in Remark 3.1 below. We emphasize that it does not represent a computational overhead as the LU -factorization of D needed to compute D -1 (K * 1 -X 0 ) is already available from the solution of the nonlinear system for the stages (K 1 , K 2 ) (see Remark 2.1).

We shall consider two choices for γ that yield mean-square A-stable integrators:

• the S-SDIRK(2,2) method for the value γ = 1 -

√ 2
2 which gives a weak order 2 A-stable method with deterministic order 2;

• the S-SDIRK(2,3) method for the value γ = 1 2 + √ 3

6 which gives a weak order 2 A-stable method with deterministic order 3.

We notice that the value of γ for the S-DIRK(2,3) yields in the deterministic case a method of order 3 which is strongly A-stable, i.e. |R(∞)| < 1, while the value of γ for the S-DIRK(2,2) yields a method of order 2 which is L-stable, i.e. it is A-stable and R(∞) = 0. L-stability is desirable in the case of very stiff deterministic problems as it permits to damp the very high frequencies.

Complexity In addition to the solution of the deterministic two stage SDIRK method (which yields the stages (K 1 , K 2 )) one step of the scheme (9) costs one evaluation of the drift function f , and 5 evaluations of each diffusion functions g r , and the generation of 2m random variables. The cost is similar to the diagonally implicit methods proposed in [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF] (in particular the number of evaluation of the diffusion functions g r , r = 1, . . . , m is independent of the number of Wiener processes m).

Remark 2.1. We emphasize that the computation of D -1 (K * 1 -X 0 ) in the scheme (9) does not represent any computational overhead. Indeed, as for any deterministic or stochastic diagonally implicit method [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF][START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF], the usual procedure for evaluating K 1 , K 2 is to compute the LU -factorization of D = Iγhf ′ (X 0 ) (f ′ (X 0 ) is usually further approximated by finite differences) and make the quasi-Newton iterations 2

LU (K k+1 i -K k i ) = -K k i + X 0 + δ 2i (1 -2γ)hf (K 1 ) + γhf (K k i ), i = 1, 2, ( 10
)
where δ 2i is the Kronecker delta function. The same LU -factorization is then used to compute

D -1 (K * 1 -X 0 ) by solving LU Y = K * 1 -X 0 ,
whose cost in negligible: the cost of evaluating K * 2 together with K * 3 is the same as one iteration of (10).

Remark 2.2. The stabilization procedure used in the method (9) can also be used to improve the stability of the strong order one methods studied in [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF] by considering the following variant of the stochastic θ-method

K 1 = X 0 + (1 -θ)hf (X 0 ) + θhf (K 1 ), K 2 = X 0 + θ(2 -θ)(1 -θ)hf (X 0 ) + θ 2 (2 -θ)hf (K 1 ) + θ 2 D -1 hf (K 1 ), (11)
X 1 = K 1 + √ h m r=1 g r (K 2 )ξ r + 1 2 m r=1 g r K 2 + m q=1 g q (K 2 )I q,r -g r K 2 - m q=1 g q (K 2 )I q,r ,
where ξ r ∼ N (0, 1) are independent random variables, and I q,r are the multiple stochastic integrals I q,r = h 0 s 0 dW q (t) dW r (s). Here, D = Iθhf ′ (X 0 ). In the case of commutative noise, these multiple stochastic integrals do not need to be simulated and can be simply replaced by h(ξ q ξ rδ qr )/2, where δ qr is the Kronecker delta function. The advantage of the integrator [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF] is that it can be shown to be mean-square A-stable for all θ ≥ 1/2 (see Remark 3.5). In contrast, the strong order one θ-Milstein method, whose stability has been analyzed in [START_REF] Higham | Mean-square and asymptotic stability of the stochastic theta method[END_REF], is mean-square A-stable only for θ ≥ 3/2 and these values of θ yield large error constants.

We next show that the integrator (9) has weak second order. Theorem 2.3. Consider the SDE (1) with f, g r ∈ C 6 P (R d , R d ), Lipschitz continuous. Then, for all fixed γ = 1/2, the integrator (9) satisfies

|E(φ(X(nh))) -E(φ( Xn ))| ≤ Ch 2 , 0 ≤ nh ≤ T, for all φ ∈ C 6 P (R d , R)
, where C is independent of n, h. 2 In the implementation, we use the initializations K 0 1 = X0 and K 0 2 = X0 + (1γ)hf (K1) and we consider the stopping criteria (

K k+1 i -K k i = 0 or K k+1 i -K k i ≥ K k i -K k-1 i
) which guaranties a convergence up to machine precision for the iterations [START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems[END_REF]. Other stopping criteria, such as K k+1 i -K k i < T ol where T ol is a prescribed tolerance could also be considered.

Proof. We base our proof on a well-known theorem by Milstein [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF], which states that under our smoothness assumptions a local weak error estimate of order r + 1 guarantees a weak order of convergence r. Since the derivative free Milstein-Talay method ( 5) is already of weak second order (see [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF]Lemma 3.1] for a short direct proof), it is sufficient to show that

|E(φ(X 1 )) -E(φ( X1 ))| ≤ Ch 3 , (12)
where X 1 , X1 are the one step numerical approximations given by ( 9) and ( 5), respectively. A Taylor expansion argument shows

K * 2 = X 0 + 3γhf (X 0 ) + O(h 2 ), K * 1 + 1 -4γ 2 hf (K * 2 ) = X 0 + h 2 f (X 0 ) + O(h 2 ),
and

K 1 = X 0 + γhf (X 0 ) + O(h 2 ), K 2 = X 0 + (1 -γ)hf (X 0 ) + O(h 2 ),
from which we deduce

h 2 f (K 1 ) + h 2 f (K 2 + √ h m r=1 g r (K * 2 )ξ r ) = h 2 f (X 0 ) + f ( K2 ) + h 5/2 R 1 + O(h 3 ), ( 13
)
where E(R 1 ) = 0. Furthermore, we notice that the last two lines of ( 5) and ( 9) are identical, with the exception that X 0 is replaced by K * 2 and K3 is replaced by K * 3 . This induces a perturbation of the form

h 2 R 2 + h 5/2 R 3 + O(h 3 ) where E(R 2 ) = E(R 3 ) = 0 and E(R 2 ξ r ) = 0
for all r (a consequence of E(J q,r ξ j ) = 0 for all indices q, r, j). We deduce

X 1 -X1 = h 2 R 2 + h 5/2 (R 1 + R 3 ) + O(h 3 ). Using X1 = X 0 + √ h m r=1 g r (X 0 )ξ r + O(h), we obtain φ(X 1 ) -φ( X1 ) = φ ′ ( X1 )(X 1 -X1 ) + O(h 3 ), = φ ′ (X 0 )(h 2 R 2 + h 5/2 (R 1 + R 3 )) + h 5/2 m r=1 φ ′′ (X 0 )(g r (X 0 )ξ r , R 2 ) + O(h 3 ).
We deduce that the local error bound ( 12) holds. To conclude the proof of the global error, it remains to check that for all r ∈ N all moments E(|X n | 2r ) are bounded uniformly for all 0 ≤ nh ≤ T, with h small enough. These estimates follow from standard arguments [18, Lemma 2.2, p. 102] using the linear growth of f, g r , assumed globally Lipschitz.

We now illustrate Theorem 2.3 numerically. In particular, we consider the following non stiff nonlinear test SDE from [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF] with a non-commutative noise with 10 independent driving Wiener processes,

dX(t) = X(t)dt + 10 j=1 a -1 j X(t) + b -1 j dW j (t), X(0) = 1, ( 14
)
where the values of the constants a j , j = 1, . . . , 10 are respectively [START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems[END_REF][START_REF] Komori | Weak second-order stochastic Runge-Kutta methods for non-commutative stochastic differential equations[END_REF][START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF]25,40,25,[START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF][START_REF] Komori | Weak second-order stochastic Runge-Kutta methods for non-commutative stochastic differential equations[END_REF][START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF]25, and the values of b j , j = 1, . . . , 10 are respectively 2, 4, 5, 10, 20, 2, 4, 5, 10, 20. For this problem, by applying Itô's formula to φ(x) = x 2 , taking expectations and using the fact that E(X(t)) = e t , one calculates E(X 2 (t)) = (-68013 -458120e t + 14926133e 2t )/14400000.We apply the S-SDIRK methods to the problem [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF] and approximate E(X 2 (T )) up to the final time T = 1 for different step sizes h. In Figure 2, we plot the relative errors using 10 9 realisations for the new integrators S-SDIRK(2,2) (solid line) and S-SDIRK(2,3) (dashed line). For comparison, we also include the results of the derivative-free Milstein-Talay method (5) (dashed-dotted line) and the DDIRDI5 method (dashed-dotted-dotted line) from [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF] (with

c 1 = c 2 = 1/2 + √ 3/6,c 3 = c 4 = 1 in [8]
). We note here that the same set of random numbers is used for all four integrators. We observe the expected line of slope 2 for S-SDIRK(2,2) (compare with the reference slope in dotted lines) which confirms the weak order two of the methods predicted by Theorem 2.3. In addition, we also observe that the methods S-SDIRK(2,3) and DDIRDI5, which have deterministic order 3, are about one magnitude more accurate than the Milstein-Talay method for steps of size h ≃ 10 -1 . As the noise is small in this example, the third order deterministic accuracy can improve the convergence of S-SDIRK(2,3) or DDIRDI5. We emphasize that this behavior does not hold in general when the diffusion function is not small compared to the drift function as can be seen in Figure 5 and its related problem. Figure 2: Weak convergence plots for the nonlinear problem ( 14) for the Milstein-Talay method (5) (dashed-dotted line), S-SDIRK(2,2) (solid line), S-SDIRK(2,3) (dashed line), DDIRDI5 (dashed-dotted-dotted line) [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF]. Second moment error at final time T = 1 versus the stepsize h, where 1/h = 1, 2, 3, 4, 6, 8, 11, 16, 23, 32. Averages over 10 9 samples.

Mean-square A-stability

In this section, we study the mean-square stability of the integrators S-SDIRK(2,2) and S-SDIRK [START_REF] Abdulle | High order weak methods for stochastic differential equations based on modified equations[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF]. The stability functions (3) have the form where P(ξ = ± √ 3) = 1/6, P(ξ = 0) = 2/3, and

R(p, q, ξ) = A(p) + B(p)qξ + C(p)q 2 ξ 2 -1 2 , DDIRDI5 [8]
A(p) = 1 + (1 -2γ)p + (γ 2 -2γ + 1/2)p 2 (1 -γp) 2 , B(p) = 1 + (1 -3γ)p (1 -γp) 3 , C(p) = 1 (1 -γp) 3 . (15)
We deduce for all p, q ∈ C,

E(|R(p, q, ξ)| 2 ) = |A(p)| 2 + |B(p)| 2 |q| 2 + |C(p)| 2 |q| 4 2 =: S(p, q). ( 16
)
Remark 3.1. We observe that by removing the term involving D -1 in the S-SDIRK methods (9), the denominators of B(p) and C(p) would scale at best as (1-γp) 2 . The resulting methods would no longer be mean-square A-stable.

In Figure 3 we visualize the mean-square and asymptotic stability regions for S-SDIRK(2,2) and compare them with the ones of the diagonally-implicit method DDIRDI5 introduced in [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF] [START_REF] Arnold | Stochastic differential equations: theory and applications[END_REF] ). The dashed lines {|p| = q 2 /2} indicate for p < 0 and p > 0 respectively the boundaries of the mean-square and asymptotic stability regions for the exact solution. As we can observe, the mean-square stability regions (dark gray) of S-SDIRK(2,2) are much bigger than the ones of DDIRDI5, and include the ones of the exact solution. This relates with the fact that DDIRDI5 and the class of weak second order methods introduced in [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF] are not mean-square A-stable as shown in [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF], whereas S-SDIRK(2,2) is a mean-square A-stable integrator, as proved in the next theorem.

] (c 1 = c 2 = 1 2 + √ 3
Theorem 3.2. The integrator S-SDIRK(2,2) is mean-square A-stable.

Proof. Since S(p, q) given by ( 16) is an increasing function of |q| 2 , we set for z ∈ C, s(z) := S(z, √ -2ℜz) and the method is mean-square A-stable if and only if

sup ℜz<0 s(z) = sup ℜz<0 |A(z)| 2 + 2|B(z)| 2 (-ℜz) + 2|C(z)| 2 (ℜz) 2 ≤ 1. ( 17
)
Using γ = 1 -√ 2/2 and putting z = x + iy, a calculation yields

∂s(z) ∂y = y a 1 (x)y 4 + a 2 (x)y 2 + a 3 (x) |1 -γz| 8 , (18)
where

a 1 (x) = 35 √ 2 -99/2, a 2 (x) = (-99 + 70 √ 2)x 2 + (-198 √ 2 + 280)x -17 + 12 √ 2, a 3 (x) = (35 √ 2 -99/2)x 4 + (-198 √ 2 + 280)x 3 + (64 √ 2 -91)x 2 + (6 √ 2 -8)x.
It can be checked that a i (x) < 0 for all x < 0 and all i = 1, 2, 3. For a given x < 0, we consider [START_REF] Milstein | Stochastic numerics for Mathematical Physics. Scientific Computing[END_REF] as a function of y, say g(y) and we observe that g(0) = 0 and g(y) > 0, for all y < 0, g(y) < 0, for all y > 0. Thus, since s(x + iy) is a smooth function of y that tends to zero for y → ∞, we deduce sup ℜz=x s(z) = s(x), for all x < 0.

Finally, an elementary study of the quantity s(x) as a function of the real parameter x < 0 yields s(x) ≤ 1. This implies the bound [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] and concludes the proof.

Remark 3.3. We observe that in the proof of Theorem 3.2 we also have lim x→-∞ s(x) = 0.

Hence the S-SDIRK(2,2) method is also mean-square L-stable.

Remark 3.4. Using [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF], it can be checked numerically that the integrator S-SDIRK(2,3) is mean-square A-stable (see an illustration for real (p, q) in Figure 4). A rigorous proof is however more tedious to derive because [START_REF] Rößler | Second order Runge-Kutta methods for Itô stochastic differential equations[END_REF] does not hold for this integrator (notice that the scheme is not L-stable for deterministic problems).

Remark 3.5. For the strong order one θ-methods (11) a simple calculation gives for the stability function [START_REF] Komori | Stable ROW-type weak scheme for stochastic differential equations[END_REF],

S(p, q) = 1 + (1 -θ)p 1 -θp 2 + q 2 |1 -θp| 4 + q 4 2|1 -θp| 4 .
Similarly as in Theorem 3.2 the mean-square A-stability of (11) can be proved for θ ≥ 1/2. We now exhibit the advantage of our method over the Milstein-Talay method (5) and the weak second order drift-implicit methods considered in [START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF]. In particular, we consider the linear test problem (2) and compare the behaviour of the three different methods for a range of parameters λ, µ for which the solution of ( 2) is mean-square stable. As we can see in Table 1, even for a moderate stiff problem (-λ = µ 2 = 5) in contrast to the S-SDIRK methods introduced here, there is quite a big stepsize restriction in order for the numerical solution to be mean-square stable for the Milstein-Talay and the DDIRDI5 methods. Furthermore, as expected we observe that the stepsize restriction for the other two methods becomes more severe as we increase the stiffness of the problem.

S-SDIRK(2,3)

We finally compare the performance of the introduced stochastic integrators on a nonlinear stiff system of SDEs with a one-dimensional noise (d = 2, m = 1),

dX(t) = α(Y (t) -1) -λ 1 X(t)(1 -X(t)) dt -µ 1 X(t)(1 -X(t))dW (t), dY (t) = -λ 2 Y (t)(1 -Y (t))dt -µ 2 Y (t)(1 -Y (t))dW (t), (20)
which is inspired from a one-dimensional population dynamics model [START_REF] Gard | Introduction to stochastic differential equations[END_REF]Chap. 6.2]. Notice that if we linearise [START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF] around the stationary solution (X, Y ) = (1, 1), we recover for α = 0 the linear test problem [START_REF] Abdulle | High order weak methods for stochastic differential equations based on modified equations[END_REF]. We take the initial conditions X(0) = Y (0) = 0.95 close to this steady state and use the parameters λ 2 = -4, µ 2 = 1, α = 1. We take for the deterministic part of the problem the stiff parameter λ 1 = -500 and we shall consider for the noise parameter µ 1 either the stiff value µ 1 = √ 500 or the nonstiff value µ 1 = 1. We plot in Figure 5 the errors for E(X(T ) 2 ) at the final time T = 1 versus stepsizes h for the integrators DDIRDI5, S-SDIRK(2,2), S-SDIRK(2,3) taking the averages over 10 8 samples. Reference solutions where computed using the Milstein-Talay method [START_REF] Alcock | Stable strong order 1.0 schemes for solving stochastic ordinary differential equations[END_REF] with stepsize h = 10 -4 . We consider the two cases of a non-stiff noise (µ 1 = 1) and a stiff noise (µ 1 = √ 500). In the non-stiff noise case (left picture), the results of S-SDIRK(2,3) are nearly identical to those of DDIRDI5 with hardly distinguishable curves, while in the stiff noise case (right picture), the results for DDIRDI5 are not included because this method is unstable for the considered stepsizes, as predicted by the linear stability analysis (see the stepsize restrictions in Table 1). It is remarkable in both cases that S-SDIRK(2,2) is more than four magnitudes more accurate than S-SDIRK(2,3) for steps with size ∼ 10 -1 , a regime for which curves with slope two can be observed. We believe that the mean-square L-stability of the S-SDIRK(2,2) method is responsible for this behavior.

MilsteinFigure 1 :

 1 Figure 1: Mean-square stability region (dark gray) and asymptotic stability region (dark and light grays).

Figure 3 :

 3 Figure 3: Mean-square stability region (dark gray) and asymptotic stability region (dark and light grays) for DDIRDI5 [8] (left pictures) and S-SDIRK(2,2) (right pictures).

Figure 4 :

 4 Figure 4: Mean-square stability region (dark gray) and asymptotic stability region (dark and light grays) of S-SDIRK(2,3).

Figure 5 :

 5 Figure 5: Weak convergence plots for the nonlinear stiff problem[START_REF] Saito | Stability analysis of numerical schemes for stochastic differential equations[END_REF] for S-SDIRK(2,2) (solid line), S-SDIRK(2,3) (dashed line), DDIRDI5 (dashed-dotted-dotted line)[START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF]. Error for E(X(1) 2 ) versus the stepsize h, where 1/h = 1, 2,[START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF][START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF][START_REF] Arnold | Stochastic differential equations: theory and applications[END_REF][START_REF] Debrabant | Diagonally drift-implicit runge-kutta methods of weak order one and two for itô sdes and stability analysis[END_REF][START_REF] Hairer | Solving ordinary differential equations II. Stiff and differentialalgebraic problems[END_REF][START_REF] Komori | Stable ROW-type weak scheme for stochastic differential equations[END_REF] 23, 32. Averages over 10 8 samples.

Table 1 :

 1 ). Comparison of mean-square stability constraints.

	SDE test problem: dX(t) = λX(t)dt + µX(t)dW (t)
		mean-square	stepsize restriction for mean-square stability
	method Milstein-Talay (5) DDIRDI5 [8] S-SDIRK (2,2) or (2,3)	A-stability no no yes	-λ = µ 2 = 5 -λ = µ 2 = 50 -λ = µ 2 = 500 h ≤ 0.236 h ≤ 0.0236 h ≤ 0.00236 h ≤ 0.246 h ≤ 0.0246 h ≤ 0.00246 no restriction no restriction no restriction

Notice that if R(p, q, ξ) = 0 with a non-zero probability, then (3) is clearly numerically asymptotically stable.

Acknowledgements. The research of A. A. is partially supported under Swiss National Foundation Grant 200021 140692.


Nadia Haddadou

Abderrezak Rachedi

Yacine Ghamri-Doudane

Trust and Exclusion in Vehicular Ad Hoc Networks: An Economic Incentive Model based Approach

In this body of work, we are interested in road safety applications such as advanced driver assistance systems, based on Vehicular Ad Hoc Networks (VANETs). One of the particular characteristics of this kind of networks is the continuous sharing of safety information by its nodes. Since this kind of information is time sensitive, a node cannot spend much time to verify its validity with an authority. However, the presence of malicious and selfish nodes in VANETs corrupts exchanged data, and lowers the overall data reception ratio in the network. To tackle this, we propose a new incentive model with exclusion for malicious nodes called VIME. VIME is inspired from the signaling theory from economics. It is based on managing a credit count that each node receives at the initialization of the application. Straightforwardly, VIME is based on two pillars. On the one hand, a node pays an appropriate cost for each sent message, which is seen by the receivers as a guarantee from the source about the truthfulness of the information. On the other hand, nodes get rewarded for cooperating in the network. The proposed economic model allows computing the amounts to be paid and those to be awarded in order to fight selfish and malicious nodes. We validate our approach via simulations. We show that VIME is able to detect and evict gradually all malicious nodes in the network, and decreases the ratio of corrupted and false sent data until reaching zero. Moreover, it has a positive impact on the participation of selfish nodes, as our approach increases the average ratio of sent data as to equal the ideal case's percentage, when no selfish node is present.

I. INTRODUCTION

Road traffic has greatly benefited from the recent network expansion. As a result, several applications have emerged having as support Vehicular Ad Hoc Networks (VANETs). Communication between vehicles is ad-hoc, using dedicated short-range communications (DSRC). This type of networks has several proprieties to take into consideration, such as high mobility, extended area, and frequent topology changes.

In this paper, we are interested in safety applications, such as advanced driver assistance systems, built on top of VANETs, where data related to road traffic is broadcasted. This kind of applications is usually time-constrained, and carries important contents whose transmission has to be completely and trusty achieved. Some solutions like ADCD [START_REF] Haddadou | Modeling and performance evaluation of advanced diffusion with classified data in vehicular sensor networks[END_REF] and EMPR-PD [2] proposed to increase the data reception ratio and to reduce the overhead as well as long latency. Nevertheless, there are still remaining security and truthfulness issues. A number of attacks are caused by the diffusion of false information or the removal of an accurate one. There are two kinds of nodes that can affect the good functioning of such systems using VANETs. The first kind is the malicious node. For example, a malicious node can constantly alternate between bad and good behavior by sharing true and false information (i.e. to avoid being detected). The second one is qualified as selfish but rational. It refuses to cooperate by relaying messages, in order to save its own resources. A large number of these nodes cause a decrease in the cooperation in a network, which weakens its connectivity and reduces its offered services.

The usual solution for establishing a trust metric between members of a network is a reputation model [START_REF] Govindan | Trust computations and trust dynamics in mobile ad hoc networks: A survey[END_REF], where members are assigned a reputation within a numeric interval. This allows nodes to have a local view of the network members, and to share it with their neighbors. A node's reputation is negatively affected among its neighbors if a bad behavior (e.g. sending false information) is detected. However, if the deployed reputation system is distributed, due to the lack of infrastructure in the roads, a node can easily be unrecognized in another area and recover a neutral reputation among its new neighbors, by simply changing directions. Due to VANET's characteristics, such as high mobility, extended deployment area, and large number of members, the lack of a global view of the network and its members makes a reputation model unsuccessful and insufficient for safety applications in VANETs.

In order to obtain a larger view and make safer exchanges, we have to impinge resources to nodes, as it is done in creditbased models [START_REF] Buttyán | Enforcing service availability in mobile adhoc wans[END_REF], [START_REF]Stimulating cooperation in self-organizing mobile ad hoc networks[END_REF] and [START_REF] Zhong | Sprite: A simple, cheat-proof, credit-based system for mobile ad-hoc networks[END_REF]. Unlike a reputation value, a credit amount does not change because of the mobility of a particular node, which allows to quickly detect and evict a malicious node based on this value. These creditbased solutions focus only on improving nodes' cooperation. To tackle these limitations, we propose VIME, a VANET Incentive Model with Exclusion for malicious nodes, to deal with both selfish nodes to improve their cooperation, and with malicious nodes to detect and exclude them from the network. VIME implements cost and reward functions to handle the account of each member according to its behavior, thus inciting nodes to act well. Each node is credited at its first connection to the network with a fixed amount of credits. For each sent message, the source pays some cost depending on different factors, such as the relevance of the information and, the node's reputation among its direct neighbors. This represents a guarantee of trustfulness for the receivers, and an investment for the source since a corresponding reward will be given to the node if the shared data is considered as valid. Therefore, a node can only receive and send messages if it has credits left, or else it is detected and excluded from the network.

The system described here is similar to the signaling theory in economics [START_REF] Spence | Signaling in retrospect and the informational structure of markets[END_REF], whose main concept is to deal with the information asymmetry in a market, between customers and sellers, regarding the quality of the offer. As in a market, the information asymmetry about the nodes' behavior in a network, between message senders and recipients, is caused by the lack of overview among them, thus inducing the emergence of malicious nodes. To detect them and to incite selfish nodes to cooperate, we propose using a signal, which is an amount with a corresponding cost for the source node, used for each sent message, to influence the neighbors' decision about the message's validity. To reward good members, some credits are given to them from their neighbors. This reward depends on the chosen signal value, and represents an incentive to cooperate and to continually behave well, as the credit is used to have access to many privileges in the network such as receiving neighbors' messages. If a node runs out of credits, then it is considered as malicious, and is evicted from the network.

The remainder of this paper is organized as follows. In Section II we present some related work in incentive solutions and their applicability to VANETs. In Section III we present VIME, our proposed approach. A simulation-based performance study is presented in Section IV. Finally, Section V concludes this paper.

II. RELATED WORK

Existing solutions about cooperation propose to motivate nodes to forward other's packets in return for a reward, in order to minimize the number of selfish nodes. This is the general concept of incentive methods using cost/reward. The approach presented in [START_REF] Buttyán | Enforcing service availability in mobile adhoc wans[END_REF] provides a solution using nuggets as method of payment for the incentive part. A node loads into its packet a number of nuggets corresponding to an estimate of the reward for the required intermediate nodes in order to reach the destination node. The main limitation of this method is the reward estimation since, in both cases of underestimation and overestimation, a part of the nuggets is lost. A second scheme presented in [START_REF] Buttyán | Enforcing service availability in mobile adhoc wans[END_REF] is the packet Trade Model. It uses a negotiation between intermediate nodes, such that an intermediate node buys a packet and sells it to the next node in such a way that it will make a benefit. But the destination node, being the last receiver, supports all the forwarding costs.

In [START_REF]Stimulating cooperation in self-organizing mobile ad hoc networks[END_REF], authors propose rules for the nodes to choose their degree of cooperation within the network according to the node's amount of nuggets. This method requires nodes to cooperate quantitatively, which does not always ensure a constant cooperation over time. Our system differs in the effort required by a node for cooperation. Reputations are frequently updated and their value have a great impact on the costs and the rewards for sending messages. This incites nodes to cooperate well and regularly, especially since a good reputation greatly reduces the cost of sending messages.

Another incentive solution proposed in [START_REF] Li | Frame: An innovative incentive scheme in vehicular networks[END_REF], introduces a sweepstakes component to enhance more nodes' collaboration. The authors incite intermediate nodes by proposing them weighted rewards, in addition to a fixed reward to one of them, chosen probabilistically. This solution depends on the existence of infrastructures, responsible of the rewarding and the reception of nodes' receipts about their cooperation actions. Nevertheless, this solution does not deal with the presence of malicious nodes.

Another approach presented in [START_REF] Zhong | Sprite: A simple, cheat-proof, credit-based system for mobile ad-hoc networks[END_REF] proposes an incentive model for mobile nodes using infrastructures as a Credit Clearance Service (CCS). Each node has to keep a receipt of its actions in the network. Hence, upon verification by the CCS, it receives a reward. The main drawback of this solution is the multiple verifications to ensure a correct functioning of the whole chain. It introduces a lot of latency and overhead which do not guarantee the scalability in VANET. In [START_REF] Sun | An identity-based security system for user privacy in vehicular ad hoc networks[END_REF], authors propose to use an identity-based cryptosystem, to allow network authorities to handle misbehaved nodes. This solution needs traceability information about nodes, obtained with the deployment of roadside infrastructures and law enforcement authorities. However, we argue that solutions are more scalable and heterogeneous if their operations do not depend on an infrastructure, which is the case of our approach, VIME.

Most existing incentive models do not deal with message diffusion. Moreover, regarding unicast, some solutions rely on estimation of rewards, which is not a useful parameter in VANETs as routes change very often; or they use a constant price for the reward, which does not always represent an incentive for the nodes. Others aim at solving either the selfish nodes problem or the malicious node problem. In our case, we believe that both should be handled at the same time. Considering these limitations, our work uses different parameters to calculate a fitted price so that each node is a winner while applying our mechanism. Moreover, it links an incentive and an exclusion module to eradicate malicious nodes and motivate selfish ones.

III. VIME: A VANET INCENTIVE MODEL WITH EXCLUSION FOR MALICIOUS NODES

VIME is a novel approach to ensure the truthfulness of shared information and handle high mobility in VANETs. It uses a credit management system to eradicate dishonest nodes, and to increase the participation of selfish nodes. VIME implies that each node has a unique identification in the network, and a tamper-proof credit count [START_REF]Trusted computing group: Tpm main specification[END_REF], which is credited initially with a fixed value. A node uses its credit for two cases. The first one is for sending a message and the second one, is for rewarding a neighbor about a shared information if it considers the content as well as the paid cost of sending it as valid. The enrichment or impoverishment of a node are related to its behavior and its degree of cooperation in the network.

A. Credit Management Functions 1) Sending Cost: According to the market strategy, a source node takes on the role of a service provider that must offer a guarantee about its products to its customers. This guarantee depends on many parameters such as the source node's reputation, a standard cost set by the application, and the importance of the data. A guarantee cost for a message is individual and differs for the same information, even between the source node and the relayed node for this message. In our environment, this guarantee is translated into a cost called Cs t i (N (i)) used to diffuse a message from node i to its neighbors N (i). This cost has to be substantial for the source node according to its initial amount of credits given at its first connection to the network, init_credit. It also has to be consistent with respect to the importance of the shared information, i_msg, set on a range of 0 to 1 as discussed in [START_REF] Haddadou | Modeling and performance evaluation of advanced diffusion with classified data in vehicular sensor networks[END_REF]. Finally, it has to meet at the same time all of the neighbors' expectations in terms of guarantee, so that the receivers consider and accept the message. The source node pays a cost as computed in equation ( 1) by decrementing its credit count via a tamperproof counter.

Cs t i (N (i)) = init_credit G × i_msg × (2 -R t i (i)) (1)
A parameter G is used for the cost calculation. It divides the initial credit that a node receives, in order to establish a reference value for message costs. G allows to increase and to decrease the number of sent messages whose cost is paid only with the initial credit. Using a small value for G, allows quickly depleting a malicious node's credit count (e.g. in case of no rewards), thus making it unable to participate in the network. Using a large value for G allows the unknown nodes, for example, to have more chances of sharing their messages in order to increase their reputation among their neighbors, especially in the case of low density and high mobility. The second parameter of the equation, i_msg, is about the diffused data relevance. Indeed, the more the data is important, the more there is consequences if it is corrupted. Therefore there is a correlation between it and the cost. The third parameter is R i (i) (i.e. a reputation about the node i held by itself at time t), based on the received reputations from its neighbors about it. The reputation values are computed as described in [START_REF] Zouridaki | E-hermes: A robust cooperative trust establishment scheme for mobile ad hoc networks[END_REF], these values are continuous and belong to the [-1, 1] range. R i (i) represents a personal and momentary estimation of the node's reputation among its current neighbors. Hence, the cost is also dependent on the reputation of the source among its neighbors, which are the receivers and the decision-makers about the truthfulness of the data. The receivers check if the paid cost corresponds to the message and their estimations of it, and then check the validity. A reputation makes the cost more or less expensive in order to offer a greater guarantee, especially when the reputation is bad. To meet the receivers' expected guarantees, and to receive rewards from them, a source node has to be able to know or estimate its reputation among them.

2) Reputation Relevance for the Reward: VIME uses a function for the reputation relevance, W(R t i (j)). It assigns a weight to R t i (j), the reputation value of a sender/forwarder node j computed by a receiver node i at time t. This weight is used for the reward estimation. It allows higher rewards for high reputation nodes compared to lower reputation nodes for the same action, thus inciting nodes to maintain good reputations.

W(R t i (j)) is computed with equation ( 2), where a hyperbolic tangent is used because of its characteristics, namely its behavior that closely resembles that of an exponential function for both positive and negative values. The tanh result is multiplied with a factor in order to obtain results in [-1, 1] range.

W (R t i (j)) = tanh(R t i (j)) × e 2 + 1 e 2 -1 (2)
3) Reward Value Computation: At the reception of a message, a node checks the cost paid by the source. To do that, the receiver calculates a range of tolerated paid cost for the sent message, [Cs t i (j) min , Cs t i (j) max ] as described in equation ( 3). This is performed in order to cope with information asymmetry in the network, i.e. in case of divergence about the reputations of each node, which may be caused by the node's mobility and topology changes.

Cs t i (i)min= init_credit G ×i_msg×(2-max(R t N (j) (i))), Cs t i (i)max= init_credit G ×i_msg×(2-min(R t N (j) (i))), (3)
The minimum boundary of this range is set to maintain a mandatory guarantee for the message. The maximum boundary is to prevent any abuse from malicious nodes. For example, they could use a greater cost for the only trusted information that they send in order to obtain a greater reward corresponding to their paid cost, thus sharply increasing their credit count. The minimum cost calculated by a receiver node j is based on the maximum held reputation of the source node attributed by it or received from its neighbors N (j), and vice versa for the maximum cost. The reward's aim is to make sure that the system remains incentive for the good nodes and dissuading for the malicious one. Moreover, it represents a cost to receive information. This decreases nodes' credit count when rewarding source nodes. Furthermore, in order to earn credits, and reward received messages, selfish nodes are encouraged to cooperate.

When the cost paid by a source node for its message lives up to the receiver's expected guarantee, the receiver examines the message and decides about its confidence on the data. If the receiver considers the message as valid, it rewards the source node by sending to it some credit via its tamper-proof credit count. The amount depends on the cost paid by the source, the density around it, which correspond to the potential rewarder nodes, and the held reputation value about the source node. The reward Rew t j (i) from a node j to a node i at time t, is then calculated as follows in equation ( 4). A reward, Rew, is estimated according to a weight , W , given to the source node's reputation, R t i (j), in order to provide higher remunerations to the nodes with greater reputation values, and to delay the credit count increase of the less reputed ones. For nodes with average reputation, their reward only recovers their paid cost, Cs t i , for the message, so there is no loss and no profit for them. However, their reputation can increase which promises them a weighty gain later. To maintain consistency, the reward value depends also on the cost paid by the source, and on the number of neighbors, N (i), which are supposed to reward the source. In order for a reward to be divided among them, even if the estimation of each one of them differs because of asynchronous information, an approximate reward can be reached, ensuring a certain level of coherence for the reward of the source node.

Rew t j (i) = (W (R t j (i)) + 1) × Cs t i (N (i)) |N (j)| (4)

B. Tamper-Proof Credit Count

To guarantee a smooth functioning for VIME, we assume that each node has a Trusted Platform Module [START_REF]Trusted computing group: Tpm main specification[END_REF], TPM, which is a secure piece of hardware with cryptographic capabilities, able to implement elliptic curves in order to generate keys according to the IEEE 1609.2 Elliptic Curve Digital Signature Algorithm [12], and to store data in shielded locations. To guarantee its security and integrity, the TPM stores a fingerprint of the application, which leads to detect any change by an attacker.

TPM is frequently used to build secure architectures for Vehicular Ad Hoc Networks [START_REF] Guette | Using tpms to secure vehicular ad-hoc networks (vanets)[END_REF]. In our case, the TPM is used to manage a tamper-proof credit 1 . The security requirements of our solution are about trustworthiness of paid costs, which aims at deducing the paid sum from the node's account, and updating its credit amount, while supporting large-scale deployment. To do this, TPM encrypts the signaling cost and reward, and stores the credit account into its shielded storage. The use of a TPM allows VIME to be distributed, and not based on the presence of any infrastructure. Fig. 1 illustrates the message exchange for either new captured information or forwarded messages, between nodes A and B using their 1 The TPM is not exclusive to VIME. It can also be used by other applications and for other purposes Fig. 2. VIME process TPM (T P M A and T P M B , respectively). Firstly, node A updates its own reputation, R T A (A), using its neighbor's shared impressions about it. This value is added to the data to share, and is sent to T P M A as M sg A , in order to calculate the corresponding cost, Cs t A , using equation [START_REF] Haddadou | Modeling and performance evaluation of advanced diffusion with classified data in vehicular sensor networks[END_REF]. Then, T P M A returns a signed message M A to node A, which is the cost added to the data encryption including a timestamp, by using its private key pk. A will then broadcast M A to its neighbors.

In our example, node B is a neighbor of A and receives M A . Foremost, node B asks its TPM about the used cost by A and compares it to the reputation that it holds about it. If the two values are coherent, according to equation ( 3), and if node B holds no information that would cast doubt on the veracity of the received data, (e.g. a too bad reputation of the source, or an incoherence in the timestamps and the localization of the event as discussed in [START_REF] Lo | Illusion attack on vanet applications -a message plausibility problem[END_REF]), then node B accepts the information, and therefore has to reward node A. Then, T P M B decrypts the data and returns it to node B, and subtracts a corresponding reward for node A according to its deployed cost and its reputation, using equation ( 4). Finally, it returns the new message, containing the signed reward value to node B, which is sent subsequently to node A, in order that its TPM increases its credit count. Thereafter, the reputation of the two nodes can be updated. A's reputation is updated because of the validity of its sent data, and B's reputation is updated for the cooperation by sending it a reward.

To avoid frauds, the calculations of the cost and the reward are made and signed by the TPM, which adds or subtracts immediately this value from the credit account.

C. Example of VIME's Operations

Fig. 2 illustrates the workings of VIME. Vehicles are able to send messages after paying a cost, calculated using equation [START_REF] Haddadou | Modeling and performance evaluation of advanced diffusion with classified data in vehicular sensor networks[END_REF]. According to their reputation, vehicle A should pay 6, and vehicle C 8, 5 when init_credit is 100, G parameter is 20, and i_msg is 1. The transaction is secured using the TPM of the sending node, which authenticates the subtraction from its credit count. Vehicle E pays more than vehicle C, which in turn pays more than vehicle A because of their different reputations. When credits are insufficient to pay the corresponding cost, as in the case of vehicle D, a node cannot send its message, and when its credit is exhausted, it is evicted from the application and considered as malicious.

As sent messages are encrypted, a receiver has to choose to accept it or to refuse it beforehand. Therefore, at each message reception, nodes verify if the cost paid by the source is coherent with equation (3). Once this verification is done, a node informs its TPM about its decision, and if it considers the data as valid, then the TPM decrypts it and delivers it to the node, while computing a reward for the source node, using equation ( 4) (e.g. vehicle B rewards A with 2.24, E with 1.62, and C with 2.37). Finally, the TPM deducts the reward from the credit count of the receiver. As for the sending node, the reward process is secured by the TPM of the receiving nodes. By performing this entire procedure, VIME encourages cooperation among nodes, and aims at automatically selecting good members in the network by proposing them interesting rewards and low costs (e.g. vehicle A earns after deducting its sending cost 6.7 in the case where the five neighbors participate in the reward process, while vehicle C earns 3.35, and E loses 2.9. Credit counts of the other vehicles diminish during the reward process, because of the reward they pay to remain informed. Conversely, they will earn credits when sending correct data. Indeed, to be able to be informed at all times, nodes have to always posses enough credits to decrypt the received messages. This is possible provided that a node cooperates and behaves well, which incites selfish nodes to cooperate more.

IV. PERFORMANCE EVALUATION

We evaluate the performance of our model by analyzing its effectiveness on detecting malicious nodes and inciting selfish nodes to cooperate, in order to improve the received ratio of truthful data and decrease that of corrupted data. We conducted a set of simulations using NS2 [START_REF]The NS-2 website[END_REF] for 100 nodes moving in a 5 kilometer highway, making round trips and some stops at both ends, according to a mobility scenario with a velocity in the range of 90 -160 km/h during 3600 seconds, such as an event is detected and transmitted each one second. We used VanetMobiSim [START_REF] Fiore | Vehicular mobility simulation for vanets[END_REF] as vehicular mobility model. The simulated radio transmission range is 250 meters, using IEEE 802.11 as the MAC layer protocol. We use ADCD [START_REF] Haddadou | Modeling and performance evaluation of advanced diffusion with classified data in vehicular sensor networks[END_REF] to optimize the diffusion of data and to minimize the communication overhead.

During the simulation, we introduce different ratios of malicious and selfish nodes to evaluate their effects on the network, and how much our solution remedies to that. For our evaluation, we selected the three following performance metrics:

• Percentage of detected nodes: Either malicious nodes or false positives (non-malicious nodes detected as malicious).

• Average ratio of corrupted data: Measures the impact of malicious nodes on the network.

• Average ratio of received data: Measures the impact of selfish nodes on the network. In our scenario, a malicious node always sends false data and, during the forwarding of a received message, it corrupts the data before retransmitting it. On the other hand, a selfish node is rational. It cooperates and forwards others' messages when it requires credits to provide for its own needs. For these simulations, we set the threshold of needed credits for a selfish node to the half of the initial credits received in the beginning, fixed to 1150 during the simulation.

First, we studied the detection percentage of malicious nodes, and the percentage of false positives in a network composed of 16% and 25% of malicious nodes. Results are presented in Fig 3 . We can evidence that the percentage of malicious node detection gradually increases, which occurs when they run out of credit. By the end of our simulated time, 100% of the total number of malicious nodes in the network have been detected and then removed from the network application. Moreover, the percentage of false positives, which is the percentage of good nodes detected as malicious, does not exceed 4% in both cases.

By analyzing in detail the simulation results, we found out that nodes that have been mistakingly evicted are those who send/forward very few packets (between 10 and 20 packets during the simulation time), compared to the number of packets they receive. In comparison, other nodes send/forward five to ten times more packets (around 100 packets). These "false positive" nodes are forwarding less packets, not because they behave selfishly, but because of the use of ADCD [START_REF] Haddadou | Modeling and performance evaluation of advanced diffusion with classified data in vehicular sensor networks[END_REF]. Indeed, due to ADCD's optimizations with the identification of target areas to the broadcasting strategy, according to the data importance in order to avoid overhead, some nodes that are in the boundary of the diffusion zone often consume packets (i.e. send rewards), but do not forward them (i.e. do not receive rewards). One possible solution for reducing the number of false positive exclusions is to have a threshold under which a node in the boundary of the diffusion zone (i.e. the last one to receive the packet) will not accept any packet as it will not be able to forward them and thus get rewarded. Even though this solution will decrease the performance of the diffusion approach, it allows avoiding the false positive exclusions. This feature can be tied to message priority. High priority messages will always be accepted, while low priority ones can be rejected.

The second studied metric is the average ratio of corrupted data sent through the network during the simulation time. Since a new event is detected every second and sent using the target diffusion protocol, ADCD, the aim is to only deliver the data to the concerned vehicles. In the case where a malicious node sends a false data or a corrupted one about a truthful previous one, a receiver well-behaving node may not detect this immediately and may thus retransmit the false data, which increases the amount of false data in the network. Therefore, it is important to rapidly detect and remove the malicious nodes from the network. Fig. 4 shows the average ratio of corrupted and false sent data in a network with different percentages of malicious nodes: 16% and 25%. We compare the performance of VIME with a network not using any solution, and network using MEB_Trust, a majority based and experience based trust algorithm adapted from [START_REF] Minhas | Towards expended trust management for agents in vehicular ad-hoc networks[END_REF], where a node always asks its neighbors about the validity of a received message before accepting it. Feature of our solution decreases this ratio until reaching 0 after 3600s in both studied cases. The decrease of the received ratio of corrupted data is due to the detection and eradication of the malicious nodes. Unlike MEB_Trust, which does not eradicate malicious nodes from the network. The third metric, represented in Fig. 5, shows the effectiveness of inciting selfish nodes to cooperate in the network. Different compositions of networks are made, with 0%, 16% and 25% of selfish nodes. The case where there are 0% of selfish nodes represents the ideal case. Note that after 2000s of simulation time, the ratio of received data in networks with 16% and 25% of selfish nodes using our incentive model is close to the ideal case, which proves the effectiveness of our approach in creating the need for selfish nodes of having credit, and thus increasing their cooperation in the network.

Furthermore, the generated messages using our solution do not interfere on the scalability of the network. Indeed, as our solution aims primarily at chiefly ensuring safety information, which generally occurs at sparse times and areas, the generated traffic is not continuous.

V. CONCLUSION

In this work, we propose VIME, an incentive model with exclusion for malicious nodes, inspired from signaling theory from economics. It is based on the management of an initial credit, that each node receives at the initialization of the application. VIME is an effective solution to cope with both malicious and selfish nodes in VANETs, without requiring the deployment of any infrastructure.It detects and evicts malicious nodes from the network when they have no more credits left. Furthermore, it incites selfish nodes to cooperate more by proposing interesting rewards. We showed via simulations that VIME detects gradually all malicious nodes in the network, while at the same time decreasing the ratio of corrupted and false sent data. Our solution also increases the participation of selfish nodes in a network, as to equal the percentage of the ideal case, with no selfish nodes in the network.

Fig. 1 .

 1 Fig. 1. Message Exchanges Using TPM

Fig. 3 .

 3 Fig. 3. Percentage of detected nodes, malicious and false positives for a network with 16% and 25% of malicious node.

Fig. 4 .

 4 Fig. 4. Average ratio of corrupted data for a network with 16% and 25% of malicious node.

Fig. 5 .

 5 Fig. 5. Average ratio of received data for a network with 0%, 16% and 25% of selfish nodes.