FRANC ¸OIS BOUCHUT Vladimir Zeitlin
email: zeitlin@lmd.ens.fr

Franc ¸ois Bouchut
email: francois.bouchut@ens.fr

Dma

A ROBUST WELL-BALANCED SCHEME FOR MULTI-LAYER SHALLOW WATER EQUATIONS

Keywords: multi-layer shallow water, nonconservative system, nonhyperbolic system, well-balanced scheme, splitting scheme, hydrostatic reconstruction, discrete entropy inequality

The numerical resolution of the multi-layer shallow water system encounters two additional difficulties with respect to the one-layer system. The first is that the system involves nonconservative terms, and the second is that it is not always hyperbolic. A splitting scheme has been proposed by Bouchut and Morales, that enables to ensure a discrete entropy inequality and the wellbalanced property, without any theoretical difficulty related to the loss of hyperbolicity. However, this scheme has been shown to often give wrong solutions. We introduce here a variant of the splitting scheme, that has the overall property of being conservative in the total momentum. It is based on a source-centered hydrostatic scheme for the one-layer shallow water system, a variant of the hydrostatic scheme. The final method enables to treat an arbitrary number m of layers, with arbitrary densities ρ 1 , . . . , ρ m , and arbitrary topography. It has no restriction concerning complex eigenvalues, it is well-balanced and it is able to treat vacuum, it satisfies a semi-discrete entropy inequality. The scheme is fast to execute, as is the one-layer hydrostatic method.

Introduction

The multi-layer shallow water model is used to describe incompressible flows in the shallow water regime, in the situation where several layers with different densities can be identified. With rotating Coriolis force, the model can be used to describe large scale atmosphere or ocean flows. In the setting of this model, the following assumptions are needed: shallowness of the layers with respect to horizontal scales, small viscosity, hydrostatics, small slope for the topography, almost uniform velocities across each layer, uniform densities for each layer.

Suitable numerical methods for solving one-layer shallow water models have been developed in the last fifteen years, and in particular for getting the so called well-balanced property, that enables to exactly resolve steady states at rest. Well-recognized methods are those of [9] and [4] (that are indeed identical at first order), in the context of the Roe solver, and the hydrostatic reconstruction method of [3].

The two-layer shallow water model have been studied more recently, and involves much more difficulties, because of nonconservative terms and nonhyperbolicity. The only works that deal with it are [9], [2], [START_REF] Abgrall | Two-layer shallow water system: a relaxation approach[END_REF], [6].

Here we introduce a variant of the splitting method of [6], that enables to correct the problem of wrong solutions with unexpected discontinuities found in [6]. The scheme involves a source-centered hydrostatic reconstruction scheme for the one-layer system, a variant of the hydrostatic reconstruction scheme. In contrast to the previous works, this new scheme allows to treat nonhyperbolicity, arbitrary number of layers, and arbitrary densities in a very robust manner, it is wellbalanced and satisfies a semi-discrete entropy inequality.

One-layer shallow water

Before writing the multi-layer model, let us first recall the main features of the one-layer shallow water model.

2.1. The system. The equations of the one-dimensional one-layer shallow water problem are:

∂ t h + ∂ x (hu) = 0, ∂ t (hu) + ∂ x (hu 2 + gh 2 /2) + hg∂ x z = 0, (1)
where h(t, x) ≥ 0 is the water height, u(t, x) ∈ R is the velocity, z(x) is the topography, and g > 0 (see Figure 1). The main properties are z(x) : topography h(t, x) : water level steady state at rest Figure 1. One layer that

• The system has a convex entropy η = hu 2 /2 + gh 2 /2 + hgz with entropy flux G = ( η + gh 2 /2)u.

• The steady states are characterized by hu = cst, u 2 /2 + g(h + z) = cst.

(2)

The steady states at rest are those for which u = 0 and h + z = cst. • The system is hyperbolic, except at resonant points (crossing of eigenvalues), defined by u = ± √ gh.

2.2. Numerical difficulties. The desirable properties of a numerical method for solving the one-layer shallow water system are of:

• keeping the water height h nonnegative,

• being able to compute dry areas where h = 0,

• preserving the total mass,

• maintaining steady states at rest (well-balanced property),

• satisfying a discrete entropy inequality,

• producing stable computations (no oscillations) for all data, including transcritical cases.

2.3. Hydrostatic reconstruction scheme. The hydrostatic reconstruction scheme proposed in [3] satisfies all of the above properties (except the entropy inequality which is only semi-discrete), and is cheap computationally. With U = (h, hu) it can be written with classical notations (see [7])

U n+1 i -U i + ∆t ∆x i F i+1/2--F i-1/2+ = 0, (3)
F i+1/2-= F l (U i , U i+1 , z i , z i+1 ), F i+1/2+ = F r (U i , U i+1 , z i , z i+1 ), (4)
F l (U l , U r , z l , z r ) = F (U * l , U * r ) + 0 p(h l ) -p(h * l ) , F r (U l , U r , z l , z r ) = F (U * l , U * r ) + 0 p(h r ) -p(h * r ) , (5)
where p(h) = gh 2 /2, F (U l , U r ) is a consistent numerical flux for the shallow water problem without source (z = cst), and the reconstructed states U * l , U * r are defined by

U * l = (h * l , h * l u l ), U * r = (h * r , h * r u r ), h * l = max(0, h l + z l -z * ), h * r = max(0, h r + z r -z * ), z * = max(z l , z r ). (6)

Multi-layer shallow water system

The one-dimensional multi-layer shallow water system can be written

     ∂ t h j + ∂ x (h j u j ) = 0, ∂ t (h j u j ) + ∂ x (h j u 2 j + gh 2 j /2) + gh j ∂ x z + k>j h k + k<j ρ k ρ j h k = 0, (7)
where h j ≥ 0, j = 1, . . . , m are the fluid depths, u j are the velocities, and z(x) is the topography. The constants g,

0 < ρ 1 ≤ • • • ≤ ρ m (8)
are respectively the gravity and the densities of the fluids. The fluids 1, . . . , m are labeled from top to bottom. The situation is represented on Figure 2 for two layers. This system admits a convex entropy, and

z(x) : topography h 2 (t, x) h 1 (t, x)
steady state at rest

∂ t j ρ j h j u 2 j /2 + gh 2 j /2 + h j gz + g j,k,k<j ρ k h k h j +∂ x j ρ j u j h j u 2 j /2 + gh 2 j + h j g z + k>j h k + k<j ρ k ρ j h k ≤ 0.
(9) The system has the steady states at rest

u j = 0, ∂ x h j + z + k>j h k + k<j ρ k ρ j h k = 0, for j = 1, . . . , m. (10) Notice that if ρ 1 < • • • < ρ m , this reduces to u j = 0, ∂ x (z + h m ) = 0, ∂ x h j = 0 for j < m, (11)
while if ρ 1 = • • • = ρ m , (10) reduces to u j = 0, ∂ x (z + h 1 + • • • + h m ) = 0. ( 12
)
As for the one-layer shallow water system (m = 1), the numerical difficulties related to this system are positivity of the depths h j , with the possibility of treating out-or in-cropping ("drying" of one ore more layers), the exact preservation of the steady states at rest (wellbalanced property), and the property to have a discrete entropy inequality. Overall, the multi-layer system has extra difficulties which are the non-conservativity of the system (even for smooth topography z) that implies the lack of appropriate Rankine-Hugoniot relations, and the possibility of having complex eigenvalues (the system is not everywhere hyperbolic). This last property raises the question of the mathematical well-posedness of the Cauchy problem. The answer is not clear, but some hope is possible, see in particular [START_REF] Pego | Instabilities in Glimm's scheme for two systems of mixed type[END_REF] for a study on a related system. In physical terms the loss of hyperbolicity corresponds to the onset of the well-known Kelvin -Helmholtz (KH) instability, i.e. instability due to the strong velocity shear between adjacent fluid layers. It was shown in [START_REF] Lesommer | Singularity formation during relaxation of jets and fronts toward the state of geostrophic equilibrium[END_REF] that criteria of hyperbolicity loss and of KH instability coincide in two-layer shallow water with a rigid lid. This argument can be easily extended to the (present) free-surface case. The same line of argument is applicaple to the sheared continously stratified hydrostatic flows, as was shown recently in [10], where the so-called Howard -Miles criterion of instability coincides with that of hyperbolicity loss. The shear instabilities persist even in the presence of dissipation, i.e. friction between the layers. In the multi-layer model, the mixing is replaced by numerical dissipation mechanisms, however of unknown nature since the system is not hyperbolic. One can hope that such dissipative modelling gives a rough description of the missing mixing, in the same spirit as the idea that a shock replaces multivalued solutions in the case of wave-breaking in fluid dynamics. Thus, the mathematical difficulty of having a mixed system corresponds to well-established physical phenomenon, which reinforces the hope of properly dealing with the Cauchy problem.

Several attempts have been made in order to solve the multi-layer shallow water system. In [9] and subsequent papers of the Malaga-Sevilla school, the two-layer case is treated by a Roe type method. A special treatment is done in order to recover positivity, and a special treatment is performed for complex eigenvalues, making the scheme inconsistent in this case. A relaxation method is proposed in [START_REF] Abgrall | Two-layer shallow water system: a relaxation approach[END_REF], with similar properties, but which is not able to treat outcropping, nor complex eigenvalues. In [2], the m-layer system is treated in the case ρ 1 = . . . ρ m , without restrictions on the eigenvalues and including drying. However, topography is not included, and the scheme is not extensible to arbitrary densities. None of these solvers is able to treat the general case with complex eigenvalues in a consistent way.

Splitting method

In [6] we introduced the operator splitting method for the multilayer shallow water system, that leads to solving successively (7) for each j, with time independent h k for k = j. This means that we solve successively m one-layer shallow water systems for U j = (h j , h j u j ), with topography

z j = z + k>j h k + k<j ρ k ρ j h k . ( 13
)
Applying a one-layer solver with numerical fluxes F l , F r gives

U j,n+1 i -U j i + ∆t ∆x i F l (U j i , U j i+1 , z j i , z j i+1 ) -F r (U j i-1 , U j i , z j i-1 , z j i ) = 0. ( 14
)
The advantages of the method are that

• The scheme is well-balanced if the resolution of each subsystem is performed with a well-balanced scheme. This is true even in the case of constant densities ρ 1 = • • • = ρ m . • The scheme is nonnegative in heights and entropy satisfying if the resolution of each subsystem is performed with an entropy satisfying scheme, because the entropy of the system is an entropy for each of the subsystems, as explained in [6]. • No estimate of the (possibly complex) eigenvalues of the full system is necessary. In particular the scheme is always consistent, and one can use the CFL condition associated to each subsystem.

4.1. Sum scheme. A variant of the splitting scheme is the "sum scheme", that can be defined as follows for two operators. Consider an ordinary differential equation

dU dt + A(U) + B(U) = 0. ( 15
)
Assume that we have two schemes for solving dU/dt + A(U) = 0 and dU/dt + B(U) = 0 respectively, given by

U n+1 -U n + ∆t A(U n ) = 0, U n+1 -U n + ∆t B(U n ) = 0. ( 16
)
Then the splitting scheme can be written

U n+1/2 -U n + ∆t A(U n ) = 0, U n+1 -U n+1/2 + ∆t B(U n+1/2 ) = 0.
(17) while the sum scheme is simply

U n+1 -U n + ∆t (A(U n ) + B(U n )) = 0. (18)
In [6], the splitting scheme (17) was used. The application of the sum scheme gives the same properties (entropy inequality, well-balanced property, nonnegativity...), and give very similar results. Note that in our context, applying the splitting scheme means that we solve (14) successively for each j, while applying the sum scheme means that we solve ( 14) simultaneously for all j.

4.2. The problem of momentum conservation. The application of the splitting method of [6], using the one-layer hydrostatic solver, leads to wrong solutions, as seen on Figure 4. It happens also with the sum scheme. In [6] we thought the problem was due to only a semidiscrete entropy inequality. Here we show that the problem is in fact due to non-conservativity of total momentum.

By summing up the momenta equations in (7) with weights ρ j , we deduce the total momentum balance

∂ t j ρ j h j u j + ∂ x j ρ j (h j u 2 j + gh 2 j /2) + j,k,j<k gρ j h j h k + j ρ j gh j ∂ x z = 0, ( 19
)
that follows from the identity

j ρ j h j ∂ x z + k>j h k + k<j ρ k ρ j h k = ∂ x j,k,j<k ρ j h j h k + j ρ j h j ∂ x z. (20)
The equation ( 19) is conservative as soon as z is smooth. Since this conservative equation is physically meaningful, one would like to use a scheme that is also conservative in total momentum. At the discrete level, this means that one would like a discrete analogue of (20), telling that the discrete nonconservative terms add up to give a conservative contribution. A natural way of achieving this is by relying on the quadratic identity

h 1 l + h 1 r 2 (h 2 r -h 2 l ) + h 2 l + h 2 r 2 (h 1 r -h 1 l ) = h 1 r h 2 r -h 1 l h 2 l , (21)
where indices l and r refer to left and right values with respect to a given interface between two cells, and h 1 and h 2 denote two different layers. Applying this several times leads to

j ρ j h j l + h j r 2 ∆z + k>j (h k r -h k l ) + k<j ρ k ρ j (h k r -h k l ) = ∆z j ρ j h j l + h j r 2 + j,k,j<k ρ j (h j r h k r -h j l h k l ), (22)
which is the desired discrete analogue of (20) (we denote ∆z = z rz l ). Since the right-hand side of ( 22) is in conservative form (in the limit ∆z → 0, which corresponds to the assumption of continuous z(x)), it is enough in order to get the discrete conservation that the nonconservative term in the momentum equation ( 7) is discretized for each j as

g h j l + h j r 2 ∆z + k>j (h k r -h k l ) + k<j ρ k ρ j (h k r -h k l ) = g h j l + h j r 2 ∆z j , (23)
where ∆z j = z j rz j l and

z j = z + k>j h k + k<j ρ k ρ j h k (24)
is the topography seen by the j's layer in the splitting scheme. However, in order not to have different times involved, such as t n+1/2 in (17) (that would break the conservation), we need to apply the sum scheme (18), that involves only the data at time n. One can observe that the one-layer hydrostatic reconstruction method (3)-(6) gives a nonconservative contribution F l -F r = (0, g h l +hr 2 ∆z), explaining the failure of conservation of total momentum.

Note that the centered discretization ( 23) is consistent with the following Rankine-Hugoniot condition for the j's momentum

h j u 2 j + gh 2 j /2 + g h j l + h j r 2 z j -s h j u j = 0, ( 25
)
where as usual [. . . ] denotes the jump of a quantity through a spacetime discontinuity, s denotes the speed of the discontinuity, and indices l and r refer to the values on each side. These Rankine-Hugoniot conditions are also those corresponding to the scheme of [9]. Thus we expect to get the same solution, even if no rigorous proof of convergence to generalized Rankine-Hugoniot conditions exists for nonconservative systems, see [8].

Source-centered hydrostatic reconstruction

As explained in the previous section, the sum scheme for the multilayer shallow water system satisfies the discrete conservation of total momentum whenever the one-layer solver used for U = (h, hu)

U n+1 i -U i + ∆t ∆x i F i+1/2--F i-1/2+ = 0, ( 26
)
F i+1/2-= F l (U i , U i+1 , z i , z i+1 ), F i+1/2+ = F r (U i , U i+1 , z i , z i+1 ), (27)
with left/right numerical fluxes F l , F r , verifies the source-centered identity

F r (U l , U r , z l , z r ) -F l (U l , U r , z l , z r ) = 0, -g h l + h r 2 (z r -z l ) . ( 28
)
We can observe that this property (28) is true for Roe type solvers like those of [4], [9], because in this context the source discretization can be chosen arbitrarily. However this property is not true for the hydrostatic reconstruction solver, that we would like to keep because of its robustness. We provide below a correction to it in order to enforce the source-centered property, but without breaking the semi-discrete entropy inequality.

5.1. Dispersive correction. In order to provide some intuition, at the continuous level, the idea is to find a dispersive correction to the shallow water system, that is neutral at the level of entropy. A parametrization of possible modified systems is as follows,

       ∂ t h + ∂ x (hu) -∂ x µu∂ x w = 0, ∂ t (hu) + ∂ x (hu 2 + gh 2 /2) + hg∂ x z -µ(∂ x w)g∂ x (h + z) -∂ x µu 2 ∂ x w = 0, (29)
with ∂ t z = 0, and µ(t, x), w(t, x) arbitrary functions. Then the entropy η = hu 2 /2 + gh 2 /2 + hgz satisfies

∂ t η + ∂ x ( η + gh 2 /2)u -∂ x g(h + z) + u 2 2 µu∂ x w = 0, ( 30
)
showing that the correction is conservative in entropy. It also keeps steady states at rest.

The aim is to find a correction µ∂ x w of order ∆x (it will be indeed of order ∆z), in such a way that the nonconservative part of the momentum equation -µ(∂ x w)g∂ x (h+ z) exactly balances the noncentered part of the numerical source coming from the hydrostatic reconstruction scheme.

At the discrete level, omitting the standard shallow water terms, a dispersive correction can be performed with U = (h, hu) as

U n+1 i -U i + ∆t ∆x i J i+1/2--J i-1/2+ = 0, (31)
J i+1/2-= J l (U i , U i+1 , z i , z i+1 ), J i+1/2+ = J r (U i , U i+1 , z i , z i+1 ). ( 32
) Denoting J i+1/2± = (J 0 i+1/2± , J 1 i+1/2± ), we require conservativity in h J 0 i+1/2-= J 0 i+1/2+ ≡ J 0 i+1/2 , (33)
and we take the momentum formulas

J 1 i+1/2-= δ i+1/2-+ J 1 i+1/2 , J 1 i+1/2+ = -δ i+1/2+ + J 1 i+1/2 , (34)
where

J 1 i+1/2 = u i (J 0 i+1/2 ) + + u i+1 (J 0 i+1/2 ) - (35)
and X + ≡ max(0, X), X -≡ min(0, X). The scalars δ i+1/2-, δ i+1/2+ parametrize the nonconservative part δ i+1/2-+ δ i+1/2+ of the momentum correction. Writing (31) by interface,

U n+1 i = 1 2 U i -2 ∆t ∆x i J i+1/2-+ 1 2 U i + 2 ∆t ∆x i J i-1/2+ , (36)
we deduce sufficient nonnegativity conditions by interface for h,

2 ∆t ∆x i J 0 i+1/2 ≤ h i , -2 ∆t ∆x i+1 J 0 i+1/2 ≤ h i+1 . (37)
Let us now write down discrete entropy inequalities. We decompose the entropy η = hu 2 /2 + gh 2 /2 + hgz as η(U, z) = K(U) + L(U, z) with

K(U) = 1 2 hu 2 , L(U, z) = 1 2 gh 2 + hgz. ( 38
)
Let us define J i+1/2 = (J 0 i+1/2 , J 1 i+1/2 ), and

U ′ i = U i -2 ∆t ∆x i J i+1/2 , U ′′ i = U i + 2 ∆t ∆x i J i-1/2 . ( 39
)
Then (37) say that U ′ i and U ′′ i have nonnegative first component. Writing

h ′ i = h i -2 ∆t ∆x i J 0 i+1/2 , h ′ i u ′ i = h i u i -2 ∆t ∆x i J 1 i+1/2 = u i h i -2 ∆t ∆x i (J 0 i+1/2 ) + -2u i+1 ∆t ∆x i (J 0 i+1/2 ) -, (40)
we deduce that

u ′ i = (1 -θ)u i + θu i+1 , for some 0 ≤ θ ≤ 1. Thus (u ′ i ) 2 ≤ (1 -θ)u 2 i + θu 2 i+1
, or in other words

h ′ i (u ′ i ) 2 ≤ u 2 i h i -2 ∆t ∆x i (J 0 i+1/2 ) + -2u 2 i+1 ∆t ∆x i (J 0 i+1/2 ) -. (41)
A similar computation leads to

h ′′ i+1 (u ′′ i+1 ) 2 ≤ u 2 i+1 h i+1 + 2 ∆t ∆x i+1 (J 0 i+1/2 ) -+ 2u 2 i ∆t ∆x i+1 (J 0 i+1/2 ) + . (42)
Therefore, one gets

K(U ′ i ) -K(U i ) + 2 ∆t ∆x i ϑ K i+1/2 ≤ 0, K(U ′′ i+1 ) -K(U i+1 ) -2 ∆t ∆x i+1 ϑ K i+1/2 ≤ 0, (43)
with

ϑ K i+1/2 = u 2 i 2 J 0 i+1/2 + + u 2 i+1 2 J 0 i+1/2 -. (44)
Using (34) and the fact that K is quadratic in hu at fixed h, (43) yields

K U i -2 ∆t ∆x i J i+1/2--K(U i ) + 2 ∆t ∆x i ϑ K i+1/2 ≤ 1 2h ′ i -4 ∆t ∆x i (h i u i -2 ∆t ∆x i J 1 i+1/2 )δ i+1/2-+ 4( ∆t ∆x i δ i+1/2-) 2 , K U i+1 + 2 ∆t ∆x i+1 J i+1/2+ -K(U i+1 ) -2 ∆t ∆x i+1 ϑ K i+1/2 ≤ 1 2h ′′ i+1 -4 ∆t ∆x i+1 (h i+1 u i+1 + 2 ∆t ∆x i+1 J 1 i+1/2 )δ i+1/2+ +4( ∆t ∆x i+1 δ i+1/2+ ) 2 . ( 45
)
Next, we treat L(U, z), which is easier since it does only depend on h and z. We have

L h i -2 ∆t ∆x i J 0 i+1/2 , z i -L(h i , z i ) + 2 ∆t ∆x i ϑ L i+1/2 = -2 ∆t ∆x i gh i J 0 i+1/2 + 2g( ∆t ∆x i J 0 i+1/2 ) 2 -2 ∆t ∆x i gz i J 0 i+1/2 + 2 ∆t ∆x i ϑ L i+1/2 , L h i+1 + 2 ∆t ∆x i+1 J 0 i+1/2 , z i+1 -L(h i+1 , z i+1 ) -2 ∆t ∆x i+1 ϑ L i+1/2 = 2 ∆t ∆x i+1 gh i+1 J 0 i+1/2 + 2g( ∆t ∆x i+1 J 0 i+1/2 ) 2 +2 ∆t ∆x i+1 gz i+1 J 0 i+1/2 -2 ∆t ∆x i+1 ϑ L i+1/2 , (46)
for some ϑ L i+1/2 to be determined. Adding up (45) and ( 46) gives

η U i -2 ∆t ∆x i J i+1/2-, z i -η(U i , z i ) + 2 ∆t ∆x i (ϑ K i+1/2 + ϑ L i+1/2 ) ≤ error, η U i+1 + 2 ∆t ∆x i+1 J i+1/2+ , z i+1 -η(U i+1 , z i+1 ) -2 ∆t ∆x i+1 (ϑ K i+1/2 + ϑ L i+1/2 ) ≤ error.
(47) Now, assume for a moment that the errors in (47) are both nonpositive. Then, from (36) we get

η(U n+1 i , z i ) ≤ 1 2 η U i -2 ∆t ∆x i J i+1/2-, z i + 1 2 η U i + 2 ∆t ∆x i J i-1/2+ , z i ≤ η(U i , z i ) - ∆t ∆x i ϑ K i+1/2 + ϑ L i+1/2 -ϑ K i-1/2 -ϑ L i-1/2 , (48)
showing that a discrete entropy inequality is satisfied with numerical entropy flux

ϑ i+1/2 = ϑ K i+1/2 + ϑ L i+1/2 . ( 49
)
But since requiring nonpositive errors in (47) is too restrictive, we shall only look for a semi-discrete entropy inequality, which is the limit of (48) divided by ∆t when ∆t → 0. For this to hold it is enough that the errors in (47) divided by ∆t give nonpositive contributions in the limit ∆t → 0. These errors coming from (45) and (46), this yields the sufficient conditions

-u i δ i+1/2--g(h i + z i )J 0 i+1/2 + ϑ L i+1/2 ≤ 0, -u i+1 δ i+1/2+ + g(h i+1 + z i+1 )J 0 i+1/2 -ϑ L i+1/2 ≤ 0, (50)
or in other words

g(h i+1 + z i+1 )J 0 i+1/2 -u i+1 δ i+1/2+ ≤ ϑ L i+1/2 ≤ g(h i + z i )J 0 i+1/2 + u i δ i+1/2-.
(51) In order to be able to find some ϑ L i+1/2 satisfying (51), it is necessary and sufficient that the lower bound is less than the upper bound, which is

g(h i+1 + z i+1 -h i -z i )J 0 i+1/2 ≤ u i δ i+1/2-+ u i+1 δ i+1/2+ . (52)
Then one can choose for ϑ L i+1/2 the half sum of the lower and upper bounds in (51),

ϑ L i+1/2 = 1 2 g(h i +z i +h i+1 +z i+1 )J 0 i+1/2 + 1 2 (u i δ i+1/2--u i+1 δ i+1/2+ ). ( 53
)
Noticing that the CFL conditions for nonnegativity (37) are satisfied in the semi-discrete limit, we conclude the following.

Proposition 1. The scheme (31)-(32) with numerical fluxes J l/r (U l , U r , z l , z r ) satisfying J l/r = (J 0 , J 1 l/r ),

J 1 l = δ l + J 1 , J 1 r = -δ r + J 1 , (54)
and

J 1 = u l (J 0 ) + + u r (J 0 ) -, (55)
is semi-discrete entropy satisfying for the entropy η = hu 2 /2 + gh 2 /2 + hgz as soon as

g(h r + z r -h l -z l )J 0 ≤ u l δ l + u r δ r . ( 56
)
Then the entropy flux can be taken as

ϑ = u 2 l 2 (J 0 ) + + u 2 r 2 (J 0 ) -+ 1 2 g(h l +z l +h r +z r )J 0 + 1 2 (u l δ l -u r δ r ). (57)
5.2. SCHR scheme. Let us now explain how to build the sourcecentered hydrostatic reconstruction (SCHR) scheme. Denoting by F HR l/r the numerical fluxes of the original hydrostatic reconstruction scheme (5)-( 6), we take

F l = F HR l + J l , F r = F HR r + J r , (58)
with J l/r = (J 0 , J 1 l/r ) defined by ( 54), (55). We have to take δ l , δ r so that J 1 r -J 1 l = -(δ l + δ r ) will kill the source error between the hydrostatic source and the desired centered source,

δ l + δ r = (F HR r -F HR l ) 1 + g h l + h r 2 (z r -z l ), (59)
that gives (28). Then we define J 0 in order to have equality in (56). Notice that even if the construction (58)-(59) looks possible for any scheme instead of F HR l/r , the issue is in the possibility to compute J 0 from (56), which somehow involves the property that h r + z rh lz l is in factor in δ l and δ r , and thus must be factorizable from the right-hand side of (59). Indeed, the well-balanced property ensures that whenever u l = u r = 0 and h r +z r -h l -z l = 0, the right-hand side of (59) vanishes. But this does not imply that h r + z rh lz l can be factorized, unless (F HR r -F HR l ) 1 is independent of u l , u r . This is wrong in general, but this is true for the hydrostatic solver. Indeed one can check that

(F HR r -F HR l ) 1 + g h l + h r 2 (z r -z l ) = g(h r -h l + ∆z) κ, (60)
with

κ = 1 2 |∆z| +                1 2 (∆z -h l )h r h r + ∆z -h l if ∆z > h l , 0 if -h r ≤ ∆z ≤ h l , 1 2 -(h r + ∆z)h l h l -(h r + ∆z) if ∆z < -h r , (61) ∆z
= min(∆z, h l ) if ∆z ≥ 0, max(∆z, -h r ) if ∆z ≤ 0. ( 62
)
One can observe that

0 ≤ κ ≤ min 1 2 |∆z|, h l + h r 2 . ( 63
)
Thus one would like to take δ l = 1 2 (1 + θ)g(h rh l + ∆z) κ, δ r = 1 2 (1θ)g(h rh l + ∆z) κ, for some -1 ≤ θ ≤ 1, and J 0 = 1 2 (1 + θ)u l + (1θ)u r κ. However, in order for the CFL condition (37) to be finite, we need that |J 0 | ≤ C min(h l , h r ). But one can check that this is not true, and only (63) is valid.

Therefore, our choice is the following. We define

κ = min κ, 5 2 min(h l , h r ) , (64)
and take

δ l = 1 2 (1 + θ)g(h r -h l + ∆z)κ, δ r = 1 2 (1 -θ)g(h r -h l + ∆z)κ, J 0 = 1 2 (1 + θ)u l + (1 -θ)u r κ, (65)
so that the semi-discrete entropy condition (56) is satisfied. The value of θ is not very important, and we have chosen an upwind formula,

θ = min 1, (u l ) + √ gh l -min 1, (-u r ) + √ gh r . ( 66
)
Now, denote by a HR (U l , U r , ∆z) the speed involved in the CFL condition for nonnegativity of the hydrostatic reconstruction scheme, and by a J (U l , U r , ∆z) the one associated to the correction scheme (31)-(32). Indeed, according to (37), one can take

a J (U l , U r , ∆z) = 2 (J 0 ) + h l + (-J 0 ) + h r . ( 67
)
We recall that according to [3], a HR (U l , U r , ∆z) = a F (U * l , U * r ) where a F stands for the speed associated to the homogeneous solver F .

Then, the CFl condition for nonnegativity for the scheme (58

) is ∆t min(∆x i , ∆x i+1 ) a SCHR (U i , U i+1 , z i+1 -z i ) ≤ 1, (68)
with

a SCHR (U l , U r , ∆z) = a HR (U l , U r , ∆z) + a J (U l , U r , ∆z). (69)
Theorem 1. The source-centered hydrostatic reconstruction scheme (SCHR) ( 26)-( 27) with numerical fluxes (58), ( 54), ( 55), ( 65), ( 66), (64), (61), ( 62) is (i) consistent with the one-layer shallow water system (1), (ii) semi-discrete entropy satisfying, (iii) well-balanced, (iv) conservative in h, and nonnegative under the CFL condition (68), (69), (67), (v) source-centered away from sharp variations of h close to vacuum.

Proof. Recall that the hydrostatic reconstruction scheme satisfies (i)-(iv). Then, according to (63), (64), (65) one has

δ l = O g∆z(h r -h l + ∆z) , δ r = O g∆z(h r -h l + ∆z) , J 0 = O (|u l | + |u r |)∆z , (70)
which yields (i) and (iii). Property (ii) follows from the above construction (by applying Proposition 1), as well as property (iv). Finally, for (v), we observe that δ l + δ r = g(h rh l + ∆z)κ. Thus, whenever κ = κ, the identity (60) gives (59), proving (28). Therefore, the only restriction is that κ = κ, which according to (64), (63), holds as soon as

1 2 (h l + h r ) ≤ 5 2 min(h l , h r ).
This inequality is what we mean by "away from sharp variations of h close to vacuum". In practice, this scheme gives numerical results for the one-layer system that are almost identical to the usual hydrostatic reconstruction scheme, except that it needs slightly smaller timesteps because of the restricted CFL condition. This can be improved by diminishing a bit the constant 5/2 in (64) (take for example 3/2), at the cost of losing a bit of source-centering for sharp variations of h close to vacuum.

In order to treat the one-layer shallow water system with transverse velocity v solving ∂ t (hv) + ∂ x (huv) = 0, (71) one has to apply also the correction J . This means that we have to complete (31) with the passive advection scheme

h n+1 i v n+1 i -h i v i + ∆t ∆x i J 2 i+1/2 -J 2 i-1/2 = 0, ( 72
) with J 2 i+1/2 = v i (J 0 i+1/2 ) + + v i+1 (J 0 i+1/2 ) -. (73)
Then the discrete entropy inequality associated to the entropy hu 2 /2 + hv 2 /2 + gh 2 /2 + hgz holds with an additional numerical entropy flux

ϑ 2 i+1/2 = v 2 i 2 (J 0 i+1/2 ) + + v 2 i+1 2 (J 0 i+1/2 ) -. (74)
Equation (71), completed with (1), corresponds to the one and half dimensional one-layer shallow water model. For treating further zeroorder right-hand sides in the momentum equation (such as Coriolis force) in a well-balanced setting, one can apply the apparent topography method of [5], see [7] for details.

5.3. Sum scheme for the multi-layer system. From Theorem 1 we get directly Theorem 2. For the multi-layer shallow water system, the sum scheme (14), ( 13) using the source-centered hydrostatic reconstruction scheme on each layer is (i) consistent for smooth solutions, (ii) semi-discrete entropy satisfying, (iii) well-balanced, (iv) conservative in h j for each layer j, and nonnegative under the CFL conditions (68), (69), (67) for each j, (v) conservative in total momentum away from sharp variations of h j close to vacuum.

It is noticeable to remark that for a sum scheme

U n+1 -U n + ∆t A 1 (U n ) + • • • + A m (U n ) = 0, (75)
if each scheme A j is stable (in the sense of invariant domain or fully discrete entropy inequality) under a CFL condition ∆t a j (U n ) ≤ 1, then the sum scheme is stable under the CFL condition

∆t m j=1 a j (U n ) ≤ 1. ( 76
)
However, since in our case the height h j evolves only through the term A j (and not through A k for k = j), the nonnegativity is true under the weak CFL condition ∆t max j=1,...,m

a j (U n ) ≤ 1, (77)
as stated in (iv), while by nature the semi-discrete entropy inequality does not need any CFL condition (because it involves only the limit ∆t → 0). Therefore, even if one would expect that for a really stable sum scheme, the CFL condition (76) is needed, in practice all our computations with a large number of layers show that (77) is sufficient, at the cost of putting a factor 1/2 in the right-hand side, or a factor 1/4 in two dimensions. Note that the one and half dimensional multi-layer model can be treated similarly by applying the sum scheme, leading to additional resolutions of (71) for each layer. Then, fully two-dimensional finite volume simulations can be performed as usual by applying interface one and half dimensional solvers. 5.4. Upwinding dispersive correction. The SCHR scheme obtained in the previous paragraph for the multi-layer shallow water system shows good numerical results, except for "upwind data" for which all eigenvalues have the same sign. For these we observe some important oscillations, as shown on Figure 3. These oscillations do not increase when refining, thus are not instabilities, but rather dispersive errors. An explanation of this phenomenon is that in the upwind situation, the scheme differs from the standard upwind scheme evaluating all the data on the upwind side, because the apparent topography z j = z + k>j h k + k<j ρ k ρ j h k is treated as a steady variable (instead of upwind).

J 1 = κu 2 l . Thus (J 0 , J 1 ) = κ(u l , u 2 l ). In order for this to compensate (83), we need to have κ = (∆z) + , which justifies (79). The case of negative eigenvalues is similar.

The obtained upwinding dispersive correction hydrostatic reconstruction (UDCHR) scheme, defined by (58), ( 54), ( 55), ( 65), (66), ( 78), (79), ( 61), (62), satisfies the points (i)-(iv) of Theorem 1, and is in some sense "weakly source-centered". The use of this scheme in the context of Theorem 2 has given the best results, even if property (v) is satisfied in a very weak sense.

Numerical tests

We only intend to show here that our approach works on the basic tests proposed in [6]. More involved tests (non-hyperbolic regime, more than two layers, second-order, two dimensions, Coriolis force) will be reported in other places.

In all our tests we take g = 9.81 with two layers, the one-layer UDCHR scheme of Subsection 5.4 being used into the sum scheme. The homogeneous solver F of (5) is taken as the Suliciu relaxation solver (the HLL solver gives similar results). A second-order two-step scheme in time is applied (Heun scheme). The weak CFL condition (77) is implemented with right-hand side 1/2. We use the first-order scheme in space, but nevertheless for Tests 1, 2, 4, 5, a reference solution computed at second-order with a very fine grid is plotted. Except for Test 1, the SCHR scheme gives similar results as the UDCHR scheme. 6.1. Test 1. It is an upwind Riemann problem with ρ 1 /ρ 2 = 0.98, z = 0, u 1 = u 2 = 2.5,

h 1 l = 0.5, h 1 r = 0.55, h 2 l = 0.5, h 2 r = 0.45, (84)
x ∈ (0, 1) and the discontinuity is taken at x = 0.5. The final time is t = 0.05 and we use 100 points. Figure 3 shows the interface z +h 2 . For comparison the SCHR scheme is plotted also (dispersive oscillations).

6.2. Test 2. It is a centered Riemann problem, ρ 1 /ρ 2 = 0.7, z = 0, u 1 = u 2 = 0, h 1 l = 1.8, h 1 r = 0.2, h 2 l = 0.2, h 2 r = 1.8, (85)
x ∈ (0, 10) and the discontinuity is taken at x = 5. The final time is t = 1 and we use 500 points. Figure 4 shows the interface. The result with the HR scheme is plotted also (unphysical discontinuity). The UDCHR scheme gives the same solution as the Roe method, indeed this test was proposed in [9]. Small perturbation of a constant state with complex eigenvalues. We take ρ 1 /ρ 2 = 0.98, z = 0, initial data are u 1 = 0.6, u 2 = -0.6,

h 2 (x) = 0.5 + 0.01 1 + cos (x-5)π 0.1 if |x -5| ≤ 0.1, 0.5 if |x -5| ≥ 0.1, (86)
h 1 + h 2 = 1,
x ∈ (0, 10), the final time is t = 1 and we use 1000 points. Figures 5 and6 show the benefit of using the UDCHR scheme.

= u 2 = 0, h 1 l = 0, h 1 r = 1, h 2 l = 1, h 2 r = 0, (87)
x ∈ (0, 1) and the discontinuity is taken at x = 0.5. We use Neumann boundary conditions, the final time is t = 0.32 and we use 100 points. The interface is plotted on Figure 7.

z(x) = 0 if x ≤ 0.5, 4(x -1/2) if x ≥ 0.5, (88)
with initial data u 1 = u 2 = 0,

h 2 (x) = 0.5 if x < 0.25, 0 if x > 0.25, h 1 (x) = (1 -h 2 (x) -z(x)) + , (89)
x ∈ (0, 1) with 100 points. Wall boundary conditions are used. Figure 9 shows the result for t = 0.5, while Figure 10 is for t = 50.

Conclusion

We have proposed a simple uncoupled method to solve the multilayer shallow water system. Its features are:

• It keeps the water heights nonnegative, and it is able to compute with dry areas, • It is well-balanced, • It satisfies a semi-discrete entropy inequality, • It works for an arbitrary number of layers, arbitrary densities, and arbitrary topography, • It can deal with complex eigenvalues consistently, indeed the eigenvalues of the system are never computed, • It runs very fast since it needs only to run a one-layer solver for each layer. We have performed simple tests in order to validate the method, that gave good results, showing the exceptional robustness of the scheme. More involved tests will be presented in other publications, including [START_REF] Scherer | Symmetric instability in twolayer rotating shallow water and its nonlinear evolution[END_REF].

Figure 2 .

 2 Figure 2. Two layers

Figure 4 .

 4 Figure 3. Test 1: upwind flow

Figure 5 .Figure 6 .

 56 Figure 5. Test 3: complex eigenvalues, HR scheme

Figure 7 .

 7 Figure 7. Test 4: mixing of two pure layers

Figure 9 .

 9 Figure 8. Test 5: initial

Figure 10 .

 10 Figure 10. Test 5: long time (t = 50)

A correction can be done at the level of the one-layer solver, by changing the definition (64) of κ into

with

The idea is that for Froude number F r = u/ √ gh small, κ = κ and the scheme is as previously, while for large Froude number, κ = (∆z) + (F r positive, i.e. positive eigenvalues), or κ = (-∆z) + (F r negative, i.e. negative eigenvalues).

The reason for putting these values is as follows. Consider the case of positive eigenvalues. Then the HR scheme (5)-(6) gives the numerical fluxes, if ∆z ≥ 0,

while if ∆z ≤ 0,

These values have to be compared to the "ideal" upwind values

Taking into account (60), this yields the errors

In (83) we can neglect g(h rh l + ∆z) κ since it is of higher order (recall that 0 ≤ κ ≤ |∆z|/2). Now it remains to take into account the contribution of the correction fluxes (54)-(55) into (58). Since in (65) δ l and δ r can also be neglected, it remains J 0 = κu l , and (55) gives


Tejeddine Mouelhi
email: tejeddine.mouelhi@telecom-bretagne.eu

Yves Le Traon
email: yves.letraon@telecom-bretagne.eu

Benoit Baudry
email: bbaudry@irisa.fr

Transforming and Selecting Functional Test Cases for Security Policy Testing

In this paper, we consider typical applications in which the business logic is separated from the access control logic, implemented in an independent component, called the Policy Decision Point (PDP). The execution of functions in the business logic should thus include calls to the PDP, which grants or denies the access to the protected resources/functionalities of the system, depending on the way the PDP has been configured.

The task of testing the correctness of the implementation of the security policy is tedious and costly. In this paper, we propose a new approach to reuse and automatically transform existing functional test cases for specifically testing the security mechanisms. The method includes a three-step technique based on mutation applied to security policies (RBAC, XACML, Or-BAC) and AOP for transforming automatically functional test cases into security policy test cases. The method is applied to Java programs and provides tools for performing the steps from the dynamic analyses of impacted test cases to their transformation. Three empirical case studies provide fruitful results and a first proof of concepts for this approach, e.g. by comparing its efficiency to an error-prone manual adaptation task.

Introduction

Automatically transforming functional test cases in order to test another concern seems an open issue, while it would drastically reduce the testing effort: no more effort would be spent other than the one already paid when testing system functions. In this paper, we explore this issue in the case of security policy testing (the new concern/property to be tested). The intuition which makes the proposed approach feasible is based on two facts:

-All security mechanisms are necessarily exercised at least once by functional test cases if they cover 100% of the executable statements: the functional test sequences can thus be reused. Only the oracle function has to be transformed to test that the tested security mechanism behaves as expected, -

The typical access control logic is separated from the application code (business logic): this access control logic can thus be manipulated or modified. Indeed, in a typical application which embeds access control mechanisms, the business logic implementation is separated from the access control logic [START_REF] Martinelli | A Runtime Monitoring Environment for Mobile Java[END_REF]. In such systems, the recommended architecture consists of designing a security component, called the Policy Decision Point (PDP), which can be configured independently from the rest of the application. The PDP is configured with respect to a security policy, modeled using an access control modeling language such as OrBAC or RBAC [START_REF] El Kalam | Organization Based Access Control[END_REF][START_REF] Ferraiolo | Proposed NIST standard for role-based access control[END_REF] A security policy is composed of a set of access control rules, which specify the conditions for granting or denying access to protected elements. The execution of functions in the business logic should thus include calls to the PDP, which grants or denies the access to the protected resources/functionalities of the system.

From a security testing point of view, the issue is to determine the correctness of the interactions between the business logic and the PDP. The generation of test cases targeting security thus consists in generating a test sequence ending with a given access control call, and also relies on the creation of the oracle function. Observing that the execution of a given test sequence effectively leads to an access granted does not mean that a security mechanism has correctly checked the access rights for this specific call. To be fully efficient, the oracle should be built by observing the piece of code which corresponds to the call to the PDP.

The problem we are tackling in this paper is related to the definition of an efficient oracle function for security test cases. In particular, we propose an approach to adapt existing functional test cases that already define relevant test scenarios, by updating their oracle to take security concerns into account.

In [START_REF] Le Traon | Testing security policies : going beyond functional testing[END_REF], we studied the overlap between functional test cases (generated for testing the functions of the business logic) and the test cases needed to test all the security mechanisms. To adapt functional test cases, in order to take security into account, we need to identify test cases that trigger security mechanisms. Then, their oracle should be enhanced to also check that the triggered security mechanisms work properly. This paper provides an automated solution for selecting functional test cases which are qualified to test security mechanisms. A dynamic analysis based on mutation of access control rules helps pinpointing which access control rules are already exercised by each test case. These rules need to be tested specifically by the security tests as outlined our previous work [START_REF] Le Traon | Testing security policies : going beyond functional testing[END_REF] Detecting the test cases impacted by a security rule (triggering at least one security mechanism) is also useful to perform regression testing on the parts of the policy that have not changed. It allows testing that the evolution has not introduced unexpected changes.

In this paper, we propose a method and tools (1) to detect the test cases which are impacted by a security policy and (2) transform existing functional test cases in security policy test cases.

To reach this objective, several problems are treated incrementally: 1. Select tests that are impacted by the security policy: this selection step aims at reducing the step 2 execution times, we perform a conservative impact analysis 2. for each selected test case, determine exactly which rules impact the test case, 3. for each test case, modify and adapt the oracle to observe the security mechanism. At the end of this three-step process, the existing functional test cases that trigger security mechanisms are transformed into security policy test cases. The presented approach has been applied to Java programs.

To validate the approach, we apply it to the case study presented in [START_REF] Le Traon | Testing security policies : going beyond functional testing[END_REF] and augmented with two other case studies. We study the feasibility of the approach and show the obtained results. In addition, we study the efficiency of our approach compared to a manual one. Also, we present to what extent the first selection step improves the performances of the approach regarding the execution times of the second step.

Section 2 introduces the context of this work and defines important concepts that are used in this study. Section presents the overall process for selecting functional test cases that trigger security mechanisms and adapt their oracle functions to take security concerns into account. Section 4 and 5 detail the two steps analyses performed to identify the relationships between functional test cases and security rules. Section 6 details how these test cases are adapted to consider security and section 7 discusses the results obtained from the three case studies.

Definitions and Context

Before detailing the proposed approach, we need to provide a few fundamental concepts. In this section, we present the context and definitions of the important concepts.

Context

When building a software system that has specific requirements, the business logic and the security policy are usually modeled, analyzed and developed separately. The security policy is specified as a set of rules, according to one of the various access control languages that are available (RBAC and variants such as XACML, OrBAC…). This abstract specification of the security policy allows reasoning about the rules and performing specific types of analyses on these rules. Eventually, the global system has to take into account both the business and the security concerns into a single implementation. The security policy can be implemented as a separate component called the Policy Decision Point (PDP) that interacts with the business logic through Policy Enforcement Points (PEP). Although the PDP can be automatically generated from the security policy, the location of PEPs in the business logic has to be decided in an ad-hoc way to take into account specific design and implementation decisions.

Figure 1 -Context for testing a secured system

Our goal is then to validate that PEPs are correct in the sense that they implement the desired interactions between the business logic and the PDP in order to enforce the security policy. In this work, we validate the integration of the security policy with the business logic through testing. More precisely we want to adapt existing functional test cases so that they can take security into account when they check the correctness of the system. Our concern in this paper is then to identify the test cases which cover PEPs, and thus trigger security mechanisms

In [START_REF] Mouelhi | A model-based framework for security policy specification, deployment and testing[END_REF], we propose a model-based approach for specifying, deploying and testing an access control policy. The specification of the policy is performed according to a generic meta-model, which allows this policy to be expressed in any access control language (RBAC, OrBAC, DAC, MAC, XACML). The metamodel includes generic mutations operators that define the semantic of mutation operators at the generic model. This framework allows us to produce mutants for any access control language that conforms to the generic model, for instance XACML. The generation of the PDP and the deployment is performed using several model-based tools (AOP, model transformation). The missing component for this fully automated framework is the generation of test cases that exercise the application code to guarantee that the access control policy is fully respected and implemented in the application internal code (no backdoors or hidden mechanisms [START_REF] Le Traon | Test-Driven Assessment of Access Control in Legacy Applications[END_REF]).

Definitions

We need to define or recall some definitions: Business Logic -The business logic is the part of the system which implements the system functionalities.

A security policy defines a set of security rules that specify rights and restrictions of actors on parts and resources of the system. A rule can be a permission or a prohibition. Each rule consists of five parameters (called entities): a status flag S indicating permission or prohibition, a role, an activity, a view, and a context. Test cases aim at detecting security flaws in the systems under test. In this paper, the faults (or "flaws") can be classified into two main categories.

Security faults -interaction faults and hidden security mechanisms. Security interaction faults occur when the interactions between the business logic implementation and the PDP is erroneous. Interaction faults thus correspond to a fault in the security mechanism itself (and maybe its absence at a given point in the software). It may be caused by manual modifications in the business logic code (e.g. adding manually an unexpected call to the PDP, erroneous parameters) or to error in the use of AOP if this technique is used to insert the PEP code. Such faults may generate unexpected interactions between the business logic and the PDP. Hidden security mechanisms correspond to design constraints or pieces of code which may bypass the PEP which is expected to control the access for a specific execution.

SP oracle function:

The oracle function for a SP test case is a specific assertion which interrogates the security mechanism. There are two different oracle functions: -For permission, the oracle function checks that the service is activated (access granted). -For a prohibition, the oracle checks that the service is not activated (access denied). The intent of the functional tests is not to observe that a security mechanism is executed correctly. For instance, for an actor of the system who is allowed to access a given service, the functional test intent consists of making this actor execute this service. Indirectly, the permission check mechanism has been executed, but a specific oracle function must be added to transform this functional test into a SP test.

Adaptation of a functional test case -The adaptation of a functional test case consists of (1) modifying its intent, (2) identifying the security rules the test sequence triggers and (3) adding the SP oracle function.

Point 2 is obtained by determining which test cases are impacted by the security policy, and more precisely by a given security rule.

Impacted Test case -A test case is said to be impacted by the security policy, if its execution triggers at least one of the PEPs in the business code. A test case is said to be impacted by a security rule R of the security policy if it triggers at least one of the PEPs with the parameters corresponding to the security rule R. A test case which is impacted by a security rule is obviously impacted by the security policy.

Overview of the approach

Figure 2 presents the methodology. First we identify the test cases which are not impacted by the security policy and remove them from the initial set of test cases. Second, we analyze the impacted test cases in order to precisely associate the security rules impacting each test case. Then we adapt each functional test case by adding a new oracle function, specific to the security mechanism under test. This new oracle function checks that the correct PEP is executed and that the access is granted/denied consistently with the security policy. In both steps, the analysis is performed dynamically. First, a single execution of all the test cases allows identifying which test case impacts at least an access control mechanism. It thus selects the test cases which are qualified for testing security, in order to reduce the effort needed to perform the second step analysis. This step optimizes the number of test cases to be run in the second step (especially when the test set is large).

In the second step, we perform a mutation analysis which requires executing the test cases for each mu-tant. We systematically inject errors in every security rule in order to create mutant PDPs. This analysis is a kind of sensitivity analysis [START_REF] Voas | PIE : A Dynamic Failure-Based Technique[END_REF], and allows locating precisely the relationships linking a test case to security policy rules in the PDP. If the behavior is the same when executing the test case on the initial and modified PDP, then it means it is not impacted. We execute every test case with every mutant. When a test case kills a mutant m, it means that it is impacted by the rule that has been modified to generate m. This mutation analysis thus allows us to associate each impacted test case with a set of security rules that are covered by this test case.

Step 1 Step 2

Mapping rules to test cases security test suite (TS3)

Test1 + security oracle We assume that the PDP is controllable. Controllability means that the rules in the PDP can be easily modified. In this paper, we consider that the PDP is a separate component. This property is needed to modify the rules and check the sensitivity of the test cases execution to changes in the PDP. The PDP security rules can be implemented in several ways, e.g. in a database or in XACML. In this paper we use the environment proposed in [START_REF] Mouelhi | A Generic Metamodel For Security Policies Mutation[END_REF] to modify the security rules in PDP implemented in several languages (OrBAC, RBAC).

Based on our previous work [START_REF] Mouelhi | A model-based framework for security policy specification, deployment and testing[END_REF], we are able to claim that this assumption is realistic. By using our generic approach, we are now able to modify the rules of policies expressed in different language (RBAC, OrBAC or XACML).

Selecting test cases that trigger security rules

In the first step of our process, we perform a dynamic analysis to select the subset of test cases that are impacted by the SP. The objective of the first dynamic analysis is to filter the test cases which are not impacted by the security policy. One test case is considered to be impacted by the SP if it triggers a PEP during its execution. It is an optimization step which is used to make this second dynamic analysis less time and resource consuming. Experiments of Section 7 will show the benefits of using this first step (test cases selection) compared to a direct use of the second step (test cases impact analysis).

Figure 3 -First step: test cases selection

As shown in Figure 2, the first test selection involves two steps.

We start by producing a mutant policy which replaces each rule from the initial policy with its opposite rule. We replace permissions with prohibitions and prohibitions with permissions, and we replace the default policy decision by its opposite as well. We end up with a mutated access control policy that contains only faulty rules. In the security mechanism, we use this policy instead of the initial one. Then tests are run on the system containing this new mutated policy in its security mechanism. Tests that exercise functionalities protected by the security mechanism are expected to fail due to the several errors injected in the policy. These tests are selected as impacted by the security policy because they fail if the policy is erroneous. This is a conservative approach, since we are sure that the test cases that kill mutants are impacted. We cannot assess that test cases that pass are not impacted by the security policy mechanism.

Relating test cases and security rules

To perform this analysis, we applied mutation adapted to access control policies, we proposed in a previous work [START_REF] Mouelhi | Mutation analysis for security tests qualification[END_REF]. The mutation approach is recalled, then the second step dynamic analysis (test cases impact on security policy rules) is detailed.

Mutation applied to security policies

Mutation analysis is a technique for evaluating the quality of test cases [START_REF] Demillo | Hints on Test Data Selection : Help For The Prac-ticing Programmer[END_REF]. We proposed several mutation operators that are applied to the security policy model and are thus independent from implementationspecific details. This approach has the advantage of defining faults that are actually related to the definition of security rules (prohibition instead of permission, wrong role, etc.). We assume that there is a direct link between these SP rules (specification level) and the corresponding PDP (implementation level). This is why the PDP has to be controllable. In [START_REF] Mouelhi | A Generic Metamodel For Security Policies Mutation[END_REF], a solution to derive directly the PDP from the security policy language has been presented based on a common metamodel to the languages and the mutation operators.

Another

The mutation operators inject mutation faults to the SP. The approach was firstly applied to OrBAC (Organization Based access control [START_REF] El Kalam | Organization Based Access Control[END_REF]), and extended to RBAC and XACML like code [START_REF] Mouelhi | A Generic Metamodel For Security Policies Mutation[END_REF]. Four types of mutation operators were proposed:

• Type changing operators: Turn one permission rule to a prohibition (PPR) or one prohibition to a permission (PRP).

• Parameter changing operators: Take one rule, and replace one of its parameters with a different one.

• Hierarchy changing operators: Replace a rule with one of its descendants.

• Adding rules operator: Add a new rule (ANR), using a combination of parameters and status that is not is the defined rules.

Detecting test cases impacted by a rule

In this second step, we map test cases to the security rules they trigger. The second step dynamic analysis is used to build these relationships. While it could be done with other techniques (e.g. by comparing the execution traces, regression testing techniques), we choose to use the mutation approach we proposed previously [START_REF] Mouelhi | Mutation analysis for security tests qualification[END_REF] to make this analysis possible. As shown in Figure 4, it has the advantage to be non intrusive in the business logic. No assumptions are needed con-cerning the PEP 'observability' and the way the security mechanisms are coded in the business logic. In fact, the business logic is treated as a black box which interacts with the PDP. The idea is to perform some "sensitivity analysis" [START_REF] Voas | PIE : A Dynamic Failure-Based Technique[END_REF] by modifying the PDP security rules and check whether the execution of a test case is impacted by this change. Thus, the mutation analysis is used to create a variant version of the PDP: it could be done using Xie's mutation technique on XACML code [START_REF] Martin | A Fault Model and Mutation Testing of Access Control Policies[END_REF] or even manually if the number of rules is not too large.

The executions (on the initial reference version and on the mutated ones) of a test case i are compared: a difference means that the change in the PDP has impacted the execution of test case i. In other words, it means that the initial execution triggers the PEP related to the security rule (recall that a PEP is a security mechanism which explicitly calls the PDP). At this stage, we don't need to know whether the initial execution result is correct or not from the security point of view (is the correct rule called? Is the result of the PDP taken into account by the business logic code?). It is just necessary to determine the test cases exercising a PEP related to a given security rule.

Figure 4 -Second step: dynamic analysis of test cases impact on security mechanisms using mutation

When a test case TC is executed, the PEP sends a request to the PDP that evaluates a security rule Rj which impacts the behavior of the test case: using mutation, we can determine that rule Rj is evaluated by the PDP when the test case is executed.

For detecting the impacted test cases by a security rule, the first operators are sufficient. The two other categories of operators could be used for a more precise analysis. The ANR operator is not used either since existing functional test cases are unlikely testing the robustness of the security policy in terms of de-fault/unspecified behavior (which is the reason why the ANR operator is used). Since this operator is costly, we do not use it, even if it may be used to identify the test case which may be transformed in advanced security test cases (see [START_REF] Le Traon | Testing security policies : going beyond functional testing[END_REF]) for testing the robustness of the security policy. We generate a set of mutant policies in the PDP by: -Applying PPR (prohibition to permission) and PRP (permission to prohibition), -Executing the test cases on each mutant. If a test case fails with a particular mutant, this reveals that the test case is actually impacted by the rule that has been mutated in this mutant.

Oracle modification in the functional test case

We propose three different levels of quality to implement the oracle function. As stated in Table 1, the first level of oracle just checks that no obvious inconsistency exists between the current security policy and the implementation of the system. Level 1 oracle function is "black-box" in the sense no information is needed to build it except the status of the rule(s) impacted by the test case. If the tested rule corresponds to a prohibition, the oracle must check that an exception or a specific message is raised. For a permission, the oracle checks that no such exception is raised. Level 2 oracle extends the basic level 1 oracle function by observing the PDP logs. The oracle then relates the oracle 1 results with the execution of the expected rule in the PDP. This allows detecting hidden mechanism. For example if the PDP is bypassed due to a hidden mechanism, this is the only way to detect that a permission has been effectively granted by the PDP. Level 2 oracle is especially useful to detect hidden security mechanisms.

Level 3 oracle extends level 1 and check that no unexpected interactions occur between the PEP which is expected to send a request to the PDP and the PDP. Level 3 oracle is thus a good way to detect interaction faults.

In this paper, we automated the weaving of level 2 oracle functions in the impacted test cases using As-pectJ. To make this process possible, the security policy is taken as an input to determine whether the expected result is 'permission' or 'prohibition' (level 1 oracle). The automation is possible only because we force the PEP, which has been woven using AOP, to raise a specific exception when the access is prohibited. With such a constraint, capturing the access status is feasible as well as comparing it with the rule status of the security policy (specification). Figure 2 presents a test case after the modification. The PDP is observable and logs coming from the business logic are observed. The security oracle retrieves the PEP log and compares it with an expected log according to the rules that are impacting the test (obtained by the second step impact analysis).

Experiments and results

The objectives of the experiments are to: -Study the information provided by the two steps.

In particular, we will study the interest of the first step of test selection in the methodology.

-compare various solutions to obtain security policy test cases: o Manual creation of SP test cases o Manual adaptation of functional test cases o Automated detection of impacted test cases and manual modification of the oracle o Fully automated approach.

Presentation of case studies

We applied our approach on three Java case studies applications previously developed by students during group projects and used in a computer science course in the university of Rennes 1:

LMS: The library management system (LMS) offers services to manage books in a public library (see [START_REF] Le Traon | Testing security policies : going beyond functional testing[END_REF] for details).

VMS: The virtual meeting system offers simplified web conference services. The virtual meeting server allows the organization of work meetings on a distributed platform.

ASMS (Auction Sale Management System): The ASMS allows users to buy or sell items online. A seller can start an auction by submitting a description of the item he wants to sell and a minimum price (with a start date and an ending date for the auction). Then usual bidding process can apply and people can bid on this auction. One of the specificities of this system is that a buyer must have enough money in his account before bidding.

The two last studies have been presented in [START_REF] Le Traon | Test-Driven Assessment of Access Control in Legacy Applications[END_REF] for studying the issue of hidden security mechanisms. Table 2 gives some information about the size of the three applications (the number of classes, methods and lines of code).

Examples

Figure 6 presents an example of the approach. An access control rule Ri mutated to R'i by switching the rule status into "permission". The mutated policy makes the test case 'testBorrowBookInMaintenanceDay' fail because the applied rule is 'permission' while the rule expected to be applied is "prohibition". In other words, the case expected an exception or an error message for borrowing a book during a maintenance day. It does not get this behaviour due to the mutated rule that was applied instead of the correct one.

Results

In this section, we show the various steps of the process. The functional test cases have been produced either manually or using the use case system testing approach of [START_REF] Nebut | Automatic Test Generation: A Use Case Driven Approach[END_REF]. The generation process stopped when 100% of code coverage was reached.

a

First step impact analysis results Table 3 shows the number/percentage of test cases which can be filtered thanks to the first step of test filtering. The execution times are less than one second for the three applications (between 0.5s and 1s) and thus can be considered as negligible in our studies. Around 20 to 30 % of the test cases can be identified not impacted by the access control policy. The first step allows saving some execution time when using the second step dynamic analysis. We measured the execution time with and without the first step of test selection. Without the first step, the impact analysis is performed on all tests cases, and then only on test cases selected by the first step. As shown in Table 4, the relative saved time is significant (and reflects the percentage of filtered test cases). This first step reduces the execution time of the second one. The execution time of the first step (test selection) is neglected. For instance, on a real enterprise information system, the execution time of a test case may be impor-tant, and a relative 20-35 % of saved execution time becomes interesting. 5 shows examples of the results that we get by this analysis. For each test case, the analysis pinpoints the exact rules that are impacting it. In addition to the previous table, it is possible to find out for each rule the corresponding impacted tests.

Furthermore, we can analyze more the results and observe the number of tests that are impacted by one rule or two or three etc. In Table 6, we show the result for the VMS application. For instance, there are 9 tests that are impacted by only one rule (each test is impacted by only one rule).

Comparing several degrees of automation

We conducted an experiment to compare between manual and automated approaches and to estimate the efforts in terms of time spent in obtaining the final test suite. We considered 4 scenarios: (1) Manually creating, from scratch, all the sec. tests, (2) Adapting manually existing tests (without step 1), (3) Adapting manually only selected tests (using the first step of the approach), (4) fully automated approach (steps 1, 2 and 3).

Two graduate students were in charge of performing these scenarios. Not all the test cases have been created manually since more than 50% of functional test cases were generated using the use case driven approach of [START_REF] Nebut | Automatic Test Generation: A Use Case Driven Approach[END_REF]. However, the students reported the time spent in the remaining manual tasks. We count in hours of intensive work, which is a lower bound of the real effort that would be needed. Table 7 displays the results showing the interest of the automation. In the three cases studies,, the creation of test cases dedicated to security is very important since it includes the identification of what should be tested and the elaboration of the test sequence. The cost of adapting the test cases selected using the steps 1 and 2 dynamic analyses is still important. Even if allows to significantly reduce the work of analysis which is needed to determine manually which functional test is impacted by a security rule, the remaining task consists in modifying systematically the oracle in the test cases. The fully automated approach is very efficient but can only be applied in a process where the PEP sends specific messages when an access is denied. The install time of the various tools and environments for performing the automated treatments has not been measured but is done once. Even taking into account the bias due to the expertise degree of the students (and all other subjective factors), the benefit of automation seems high.

Related works

As far as we know, no previous work studies how to automatically adapt and reuse functional tests for SP testing. However, several works proposed techniques and tools for automatically testing the PDP implementation for security policies written in XACML [START_REF] Martin | Automated Test Generation for Access Control Policies via Change-Impact Analysis[END_REF][START_REF] Vincent | Conformance Checking of Access Control Policies Specified in XACML[END_REF] or RBAC [START_REF] Masood | Technical report: Scalable and Effective Test Generation for Access Control Systems that Employ RBAC Policies[END_REF]. Fisler et al. proposed Magrave a tool for analyzing XACML policies and performing change-impact analysis [START_REF] Fisler | Verification and change-impact analysis of access-control policies[END_REF]. The tool can be used for regression testing to identify the differences between two versions of the policies and test the PDP. In [START_REF] Martin | Automated Test Generation for Access Control Policies via Change-Impact Analysis[END_REF], Xie et al. proposed a new tool Cirg that automatically generates test for XACML policies using Change-Impact Analysis.

The main difference between their work and ours is that they focus on the PDP alone and testing consists of sending request and getting responses that are compared to the expected ones. We consider the whole system (PDP+PEP) and adapt functional tests to validate the implementation of the security policy. Both works are thus complementary. In fact, it is necessary to begin with making sure that the implementation of the PDP is correct and that the PDP responds as expected according to the specification, before considering the whole system. However, it is obvious that testing the PDP in a separate way does not guarantee that the application behind it actually implements correctly the access control policy. We cannot assume that the application behaves as expected by the access control policy by only checking the front door (the PDP+PEP). For instance, the internal application code (the business code) may contain backdoors that bypass the PDP. Testing the PDP independently does not allow these backdoors to be detected and removed. In addition, as we have previously pointed out [START_REF] Mouelhi | A model-based framework for security policy specification, deployment and testing[END_REF], applications may contain hidden doors (that we called hidden mechanisms) which implement access control rule in a nondocumented ways. These hidden doors can be detected only by testing the whole system w.r.t. the access control policy. For all these reasons, it is absolutely important to test the whole system.

Even if the idea of estimating the impact of test cases is not new (see [START_REF] Al-Hayek | Impact of system partitioning on test cost[END_REF] for design choices), no previous work specifically studies how to automatically adapt and reuse functional tests for SP testing.

Conclusion and Future Work

In this paper, we presented a new automated approach for selecting and adapting the functional tests in the context of SP testing. The approach includes a three steps process. The first step uses a special mutant of the policy having all its rules mutated. This step helps selecting the subset of tests impacted by the access control policy. The second step uses single mutations to provide a mapping between tests and the triggered SP rules. According to the mapping, AOP is used to transform them into security test cases by weaving the security policy oracle function. The approach was successfully implemented and applied to three case studies. The experimental results show the effectiveness of the approach when compared to a manual one.

In further work, we will study regression testing in the case of a security policy evolution: the second step (impact analysis) can be used for selecting and adapt-ing security tests with respect to the security policy modification. In addition, one of the main issues is access control tests automation. This issue can be addressed using combinatorial [START_REF] Pretschner | Model-Based Tests for Access Control Policies[END_REF] or computational intelligence algorithms [START_REF] Baudry | Automatic Test Cases Optimization using a Bacteriological Adaptation Model: Application to[END_REF][START_REF] Baudry | From Genetic to Bacteriological Algorithms for Mutation-Based Testing[END_REF] and will be studies in future work.

Figure 2 -

 2 Figure 2 -Overview of the three steps

  i is impacted by Rule j

Figure 5 -

 5 Figure 5 -Woven security oracle function

Figure 6 -

 6 Figure 6 -Example of rule and impacted test case

  PEP sends a request to the PDP to get the suitable response for the requested service by the user and in the current context. Based on the response of the PDP, if the access is granted the service executes, and if the access is denied the PEP forbids the execution of the service.PDP -ThePolicy Decision Point is the point where policy decisions are made. It encapsulates the Access Control Policy and implements a mechanism to process requests coming from the PEP and return a response which can be deny or permit. System/functional testing -the activity which consists of generating and executing test cases which are produced based on the uses cases and the business models (e.g. analysis class diagram and dynamic views) of the system. By opposition with security tests, we call these tests functional.

	Our domain consists of role names RN, activity names
	PN, view names VN, and context names CN.
	Security Policy (SP): A security policy SP is thus set of
	rules defined by	SP	⊆	S	×	RN	×	PN	×	VN	×	CN	. The
	signature of an access rule is thus:
	Status(Role, Activity, View, Context)
	Examples:
	Permission(Borrower,Borrow,Book,WordingDay)
	Prohibition(Borrower,Consult,PersonnelAccount,
	Default)
	Security mechanism -It denotes any piece of code or
	constraint internal to the business logic which re-
	stricts (or relaxes) the access to some protected re-
	sources/functions of the system. A specific case of
	security mechanism is called PEP, which explicitly
	calls the external PDP component (see the follow-
	ing definitions).
	PEP -The Policy Enforcement Point is the point in the
	business logic where the policy decisions are en-

forced. It is a security mechanism, which has been intentionally inserted in the business logic code. On call of a service that is regulated by the security policy, the Security policy testing (SP testing): it denotes the activity of generating and executing test cases that are derived specifically from a SP. The objective of SP testing is to reveal as many security flaws as possible.

Test case: In the paper, we define a test case as a triplet: intent, input test sequence, oracle function.

Intent of a test case:

The intent of a test case is the reason why an input test sequence and an oracle function are associated to test a specific aspect of a system. It includes at least the following information: (functional, names of the tested functions) for functional test cases or (SP, names of the tested security rules) for SP ones.

Table 1 . The three quality level of oracle functions

 1

	Oracle	Assumptions	Check that..
	level
			the	access	is
	1	no	granted/denied w.r.t. the
		(black box)	rules (1)
			-(1)
		Observable	-the PDP is correctly
	2	PDP	called (2)
		("Business
		logic + PEPs"
		is a black box)
			-(1)
		Observable	-the PDP and the ex-
	3	PDP and	pected PEP are correctly
		Observable	interacting (the right PEP
		PEP	calls the PDP with the
		(glass-box)	expected parameters).

/ test if data was well stored in DB bookReturned

	/**
	* Test for method borrow
	*
	*/
	public void testBorrow() {
	try {
	// init test data
	Book book = new Book("book title"); Teacher teacher = new Teacher("teacher1");	Test sequence
	// teacher borrow book
	bookService.borrowBook(teacher, book);
	// functional oracle
	// test if borrowed and no longer reserved
	assertTrue(teacher.getBorrowed().contains(book));
	/= bookDAO.loadBook("book title");		Existing
	asserEqualsTrue(bookReturned.getCurrentStateString() ,Book.BORROWED));	functional oracle
	} catch (BSException e) {
	fail(e.getMessage());
	} catch (SecuritPolicyViolationException e) {
	fail(e.getMessage());
	}
	// Security oracle
	// get the rules impacting this test
	String rules = getRulesImpactedByTests("testBorrow"); // get both expected and return log	Woven level 2
	String peplogExpected = getExpectedPEPLog(rules); String pepLogReturned = readPEPLog();	security oracle
	// security oracle
	assertEquals(pepLogReturned,peplogExpected);
	}

Table 2 . The size of the three applications

 2

		# classes	# methods	LOC
				(exec. statements)
	LMS	62	335	3204
	VMS	134	581	6077
	ASMS	122	797	10703

Table 3 . First step analysis results

 3

	System	Impacted	Other Tests	All
		Tests
	VMS	41	11 (21%)	52
	ASMS	85	38 (31%)	123
	LMS	19	7 (27%)	26

Table 4 . Execution time with/without the first analysis

 4

		Application	Without 1 st	With 1 st	% Saved
			step	step	time
		VMS	30.5 s	24,3 s	20.4 %
		SMS	114.1 s	74.6 s	34.6 %
		LMS	19.3	16.4	15 %
	b	Second step: test impact analysis results
		Table

Table 5 . Examples of detailed results

 5

	Tests	# rules	Rules
			1. Permission(Teacher, GiveBack,
			Book, WorkingDays)
			2. Permission (Teacher, Borrow,
	Test for the LMS: testBorrow2Borrower()	4	Book, WorkingDays) 3. Permission (Student, GiveBack,
			Book, WorkingDays)
			4. Permission (Student, Borrow,
			Book, WorkingDays)
			1. Permission (Webmaster, Update,
			Account, MaintenanceDay)
	Test for VMS: testUpdateUserAccount()	3	2. Permission (Admin, Update, Account, MaintenanceDay)
			3. Permission (Personnel, Update,
			Account, MaintenanceDay)
			1. Permission (Seller, CreateSale,
	Test for ASMS: testOpenSale_0()	2	Bid, default) Permission (Seller, Update, Bid,
			default)

Table 6 . Tests and impacting rule for VMS application

 6

	# Tests	# Impacting Rules
	9	1
	12	2
	2	3
	9	4
	1	5
	5	6
	1	7
	1	8

Table 7 . Automated and manual approaches

 7

		LMS	VMS	ASMS
	(1)Creating	32 hours	48 hours	64 hours
	all tests
	(2)Adapting	8 hours	16 hours	24 hours
	manually all
	existing tests
	(3)Adapting	1 hour	3 hours	4.5 hours
	selected tests
	(4)Fully	5 min + install	10 min +	13 min +
	automated	time	install time	install time

Acknowledgments: This work is supported by "Région Bretagne" (Britanny Council) through a contribution to a student grant.