Eitan Altman
email: eitan.altman@inria.fr

Aniruddha Singhal

Corinne Touati
email: corinne.touati@inria.fr

Jie Li
email: lijie@cs.tsukuba.ac.jp

Resilience of Routing in Parallel Link Networks

We revisit in this paper the resilience problem of routing traffic in a parallel link network model with a malicious player using a game theoretic framework. Consider that there are two players in the network: the first player wishes to split its traffic so as to minimize its average delay, which the second player, i.e., the malicious player, tries to maximize. The first player has a demand constraint on the total traffic it routes. The second player controls the link capacities: it can decrease by some amount the capacity of each link under a constraint on the sum of capacity degradation. We first show that the average delay function is convex both in traffic and in capacity degradation over the parallel links and thus does not have a saddle point. We identify best responses strategies of each player and compute both the max-min and the minmax values of the game. We are especially interested in the min max strategy as it guarantees the best performance under worst possible link capacity degradation. It thus allows to obtain routing strategies that are resilient and robust. We compare the results of the min-max to those obtained under the max-min strategies. We provide stable algorithms for computing both max-min and min-max strategies as well as for best responses.

Introduction

The current computer networks such as Internet architecture remain remarkably vulnerable to different security attacks and failures which may cause system unaivabilities or performance degradation. It is a great challenge to provide services under such security attacks and failures in computer networks. Resiliency is the ability to provide and maintain an acceptable level of service in the face of faults and challenges to normal operation [START_REF] Smith | Network resilience: a systematic approach[END_REF].

In this paper, we study the resilience problem of routing traffic in a parallel link network model with a malicious player using game theory. Although the network model looks simple, it could be taken as a typical one for a computer network with general network configuration in which there are many paths between a source node and a destination node and a path consists of several communications.

Although our network is a simple one, the network resilience problem in the network model is not a trivial one. We study the resilience problem of routing traffic in a parallel link network model with a malicious player using a game theoretic framework. Consider that there are two players in the network: the first player wishes to split its traffic so as to minimize its average delay, which the second player, i.e., the malicious player, tries to maximize. The first player has a demand constraint on the total traffic it routes. The second player controls the link capacities: it can decrease some amount of the capacity of each link under a constraint on the sum of capacity degradation. We first show that the average delay function is convex both in traffic and the capacity degradation over the parallel links and thus does not have a saddle point. We identify best responses strategies of each player and compute both the max-min and the minmax value of the game. We are especially interested in the min-max strategy as it guarantees the best performance under worst possible unknown link capacity degradation. It thus allows to obtain routing strategies that are resilient and robust. We compare the results of min-max to those obtained at max-min. We provide numerical algorithms for computing both max-min and min-max values and strategies as well as for best responses.

Related work

We restrict in this paper our analysis to the framework of routing in a "parallel link" network. This topology has long been a basic framework for the study of routing, as it is a natural generic framework of load balancing among servers in a network. The study of competitive routing in networks with parallel links using game theory goes back to [START_REF] Orda | Competitive routing in multiuser communication networks[END_REF]. They were further studied in [START_REF] Harks | Stackelberg strategies and collusion in network games with splittable flow[END_REF], [START_REF] Koutsoupias | Worst-case equilibria[END_REF] and many others. The only reference we know that studied adversarial behavior in routing in a model similar to the max-min scenario is [START_REF] Blocq | Worst-case coalitions in routing games[END_REF] but they do not propose an algorithmic solution as we do here. On the other hand, to the best of our knowledge, the min-max setting has not been studied before. While the maxmin problem has a water-filling structure, we show that the min-max policy has a form which extends the water filling policy and we call it the "water distribution" policy. We provide an algorithm for computing it.

System Model and Problem Formulation

Consider a set L = {1, ..., L} of parallel links between a common source s and destination d as shown in Fig. 1.

Let the delay density over link ∈ L of capacity C be given by the following function of the link flow x l :

D(x , C )    1 C -x if x < C , +∞ otherwise. (1)
Let x be the flow vector, x = (x , 1 ≤ ≤ L). Define the system delay as the average delay experienced by the flow on the different links:

D(x, C) ∈L x D (x , C ).
(2)

2 1 L s d • • •
Fig. 1: A system of parallel links Such delay system model was already widely used to describe delay in telecommunication network (see, e.g. [START_REF] Orda | Competitive routing in multiuser communication networks[END_REF]). In this paper, we address resilience in routing for such networks.

The vector x is controlled so as to minimize the system delay under the demand constraint: ∈L x = X. Meanwhile, suppose that the capacity C of link is decreased to C -δ where δ ∈ R + . In this case, the delay of link becomes

D(x , C -δ ) =    1 C -δ -x if x < C -δ , +∞ otherwise.
Let δ be the degradation vector, δ = (δ , 1 ≤ ≤ L). The average system delay is therefore given by

D(x, C -δ) ∈L x D (x , C l -δ ), (3)
and the worst degradation is therefore max δ D(x, C-δ) subject to the constraint

∈L δ = ∆.
In this paper, we study a routing that would be robust under the worst possible impact of removing an amount ∆ of link capacities. Our objectives are as follows:

Objective 1 For a given load vector x, identify the vector δ which is the most harmful. This is addressed in Section 3.

Objective 2 Conversely, for a given attack δ, identify the best possible response x to capacity degradation. This is addressed in Section 4.

Objective 3 Determine the worst possible capacity degradation to arbitrary flow vector. This is addressed in Section 5. More precisely, our problem can be interpreted as a zero-sum game with x playing the role of the minimization player's action and δ the maximization player's. If player Π δ is playing first, then it will aim at finding the attack which reduces the system capacity most:

max δ min x D(x, C -δ) subject to        ∈L x = X and ∈L δ = ∆. ( 4
)
Objective 4 Determine the flow vector x * which guaranties the best possible performance under any possible capacity degradation response. This is addressed in Section 6. That is, if player Π x is playing first, it will aim at choosing the flow vector x * that guarantees the best possible performance under the worst possible reaction of attack δ of player Π δ :

min x max δ D(x, C -δ) subject to        ∈L x = X and ∈L δ = ∆. (5)
A crucial question is whether the solutions of the latter two problems coincide. The following result gives a clue: Proposition 1. The average delay function D is convex both in x and δ.

The proof is available in the technical report [START_REF] Altman | Resilience of routing in parallel link networks[END_REF] A very well studied class of games is that of concave-convex games, for which the maximizing player's optimization function is concave while the minimizer player's is convex. These games are known to have a value, that is their maximin optimization (4) coincide with their minimax [START_REF] Blocq | Worst-case coalitions in routing games[END_REF]. However, in our scenario, the game is convex-convex, and therefore the order at which the players are taking decisions can affect the resulting equilibrium.

In the following, we shall obviously restrict to the case that ∈L C > X+∆ so that there exists a routing strategy with finite cost.

Optimal Attack in response to Link Utilization

In this section, we consider the optimal strategy for player Π δ in response to a given link usage x. That is:

δ * (x) arg max δ≥0 D(x, C -δ), s.t δ = ∆.
The next theorem gives a characterization of the optimal reaction of player Π δ : it is such that only a single link should be attacked, that is, the one inducing the higher throughput degradation:

Theorem 1 (Optimal attack response). For any load vector x, there exists a unique optimal reaction of player Π δ . It is such that:

δ * (x) = ∆ if = * 0 otherwise, with * = arg max ∈L x (C -x )(C -∆ -x ) . (6)
Proof. For a given x vector, note that D is convex in δ and defined on the convex polytope P :

P {(δ ) ∈L , ∀ , 0 ≤ δ ≤ ∆ and δ = ∆}
Define e i the unit vector, i.e. the vector of dimension L with all elements being equal to 0 except for the ith element which is 1. Then P is the convex hull of a set of L + 1 extreme points: {0 ∪ ∆e i , i ∈ L}.

Hence, for any point p of P , there exists non-negative α 0 , ..., α

L such that 1 = ∈L α and p = ∆ α e . Let * = arg max D(x, C -e ). As D is convex, then D(x, C-p) = D(x, C-∆ α e ) ≤ ∆ α D(x, C- e ) ≤ ∆ α D(x, C -e * ) = ∆ D(x, C -e *
) which gives Eq. ( 6). The degradation of the delay induced by attacking link is

x D (∆, x ) -D (0, x ) = x ∆ (C -∆ -x ) (C -x )
which leads to the desired result.

Corollary 1. The degradation induced by player Π δ on the total delay equals to

∆x * (C * -x * )(C * -∆ -x * )
.

Therefore, from Theorem 1, a straightforward algorithm can give the exact optimal attack in 3 × L multiplications and L comparisons.

Optimal Link Utilization as Response to some Attack

We now analyze the optimal link utilization in response to a choice of δ from player Π δ , which is denoted by x * (δ). Then, we seek:

x * (δ) = arg min x≥0 D(x, C -δ), subject to x = X.
The next theorem gives a characterization of x * (δ):

Theorem 2 (Optimal link usage response). There exists a unique real value K such that:

x * (δ) =    C -δ -K C -δ if ∈ X , 0 otherwise, (7)
with

X { , C -δ ≥ K 2 } and K ∈X (C -δ ) -X ∈X C -δ . (8)
(Note the fixed point equation between the set of links used at the optimal response X , and the quantity K.)

Proof. For given C and δ, consider the Lagrangian function

L(λ, x) = D(x, C -δ) -λ ∈L x -X . (9)
Then, the optimal link usage response x * is solution of the optimization problem min

x≥0 L(λ, x). Since x → L(λ, x) is convex, then x * = arg min x L(λ, x) ⇔ ∂L(λ, x * ) ∂x *    ≥ 0 ∀ = 0 if x * > 0. Then          x * = 0 ⇔ ∂ L ∂x (λ, x|x = 0) ≥ 0 ⇔ 1 C -δ ≥ λ, x > 0 ⇔ C -δ -x = C -δ λ (10)
which gives (7) by taking K = 1/ √ λ. Then, summing Eq. (7) over X yields X = ∈X C -δ -K C -δ , which allows us to express K as Eq. (8).

From this theorem, we can then derive the performance achieved at the optimal link usage response: Proposition 2 (Performance at the optimal link usage). At the optimal x * (δ), the total delay on any used link (i.e. such that x * > 0) is given by

x * (δ)D (x , C -δ ) x * (δ) C -δ -x * (δ) = √ C -δ K -1
and the total delay is

D(x * (δ), C -δ) = ∈X C -δ 2 ∈X (C -δ ) -X -|X |. (11)
Proof. From (7), we have:

x * (δ) C -δ -x * (δ) = C -δ -K √ C l -δ K √ C -δ = √ C -δ K -1. Thus D(x * (δ), C -δ) = ∈X √ C -δ K -1 = ( ∈X √ C -δ ) 2 ∈X (C -δ ) -X -|X |.
In order to derive a powerful algorithmic solution, we need the following characterization of the optimal link usage solution:

Proposition 3 (Optimal Usage Characterization). For each link , define the normalized delay as

ND (x, C -δ) = C -δ .D (x, C -δ). ( 12
)
Then, at the optimal x * (δ):

ND (x * (δ), C -δ) = K(X ) if C -δ ≥ 1/K 2 ≥ K(X ) if C -δ ≤ 1/K 2 (13)
Proof. At the optimal response x * (δ), we have, from Eq. ( 7), for any used link:

ND (x * (δ), C -δ) = √ C -δ C -δ -x * (δ) = √ C -δ K √ C -δ = 1/K.
For any unused link, we have:

ND (x * (δ), C -δ) = √ C -δ .D (x * (δ), C -δ) = √ C -δ C -δ = 1 √ C -δ But from Eq. (10), 1 C -δ ≥ λ, i.e. 1 √ C -δ ≥ √ λ = 1/K which concludes the proof.
The proposed water-filling mechanism for the strategy of player Π x is given in Algorithm 1. The links are initially sorted by decreasing capacity. The mechanism gradually increases the amount of x of the various links until reaching X.

More precisely, the algorithm proceeds with initialization of x as zero. At each iteration of the algorithm, the set X is updated by checking for some potential new candidates.

Algorithm 1: The algorithm of player Π x which defines an optimal strategy in response to some attack.

Input: Vector C of channel capacities, Attack vector δ Output: Load vector x 1 Sort links with decreasing capacity C -δ 2 T A ← 0 // The traffic allocated so far 3 Link ← 1 // The link index up to which we inject data 4 ε ← 0.001 // Set it to the desired accuracy 5 x ← 0 // The traffic vector 6 while T A < X do

7 while Link < L and ND1(x, C -δ) ≥ ND Link +1 (x, C -δ) do 8 Link ++ 9 x1 ← x1 + ε 10 K ← √ C1 -δ1 C1 -δ1 -x1 11 for j = 2 to Link do 12 xj ← Cj -δj - √ C j -δ j K 13 Update T A ← ∈L x 14 return x
One can use a direct water-filling algorithm by using ε, a very small quantity, representing the discretization of level increase in the water-filling algorithm. The algorithm would be a direct implementation of Proposition 3, that is, if the current link was filled up to a level (in terms of ND) that is greater or equal than that of next link, then variable Link is incremented so as to start filling the next link. Then, the "for" loop would fill each link j by a small amount η j which is such that ND j (x + η j e j , Cδ) -ND j (x, Cδ) = ε until X is exhausted.

The performance of such algorithm exhibits average performance though, as the numerical precision errors in the level increases of the different links are summed up over the different iterations and can end up in large inaccuracy if the ratio x /η turns out to be large. Performance is significantly improved by using one link (for instance that of greater capacity) as the point of reference of the link level and setting up the other links levels accordingly. We propose another variant of the algorithm where ε represents the discretization of x 1 . Then, at each iteration of the algorithm, x 1 is increased by ε, then K is updated and then all links in X .

Then, the maximal error is Err ≤

∈L |x (x 1 )-x (x 1 +ε)| = ∈L | √ Cj -δj K(x1) - √ Cj -δj K(x1+ε) | = ε √ C1-δ1 ∈L √ C -δ .
Since the links are ordered by decreasing capacity, then the error is bounded by Lε.

Optimal Link Degradation Strategy to Unknown Link Usage

Let us now consider that player Π δ is to choose its attack vector, without knowing the link usage chosen by player Π x . Then, a natural strategy of player Π δ is to choose the attack vector that would guarantee the highest value of delay under any action load x. Such strategy of player Π δ is commonly known as the maxmin strategy and is given by Definition 1.

Definition 1. δ * is a maxmin strategy if it is solution of Mm(C, X, ∆) : max δ≥0 min x≥0 D(x, C -δ), s.t x = X, δ = ∆. ( 14
)
Note that this is equivalent to a two-player sequential game where player Π δ plays first, followed by player Π x , after it observes the action of player Π δ .

Existence and Characterization of the Optimal Strategy

Theorem 3 shows that there exists a unique strategy for player Π δ and provides a characterization of it. Theorem 3. There exists a unique real value α such that the optimal strategy for player Π δ is given by:

δ = C -α if ∈ D, 0 otherwise with α = ∈D C -∆ |D| and D = { |C ≥ α}. ( 15
)
The proof is given in the technical report [START_REF] Altman | Resilience of routing in parallel link networks[END_REF]. Note that the optimal strategy for player Π δ is therefore to attack the links of greater capacity in a way so that their remaining capacities (C -δ ) are all equal to α. Hence, the optimal strategy for player Π δ is independent on the weight X of player Π x .

A decreasing water filling algorithm

Based on Theorem 3, we can derive an algorithm to compute the optimal strategy of player Π δ , which is given in Algorithm 2.

Similarly to Algorithm 1, at each step of the algorithm, the links 1 to Link are being filled. The algorithm ends whenever all links have been attacked or when the attack level ∆ has been exhausted. More precisely, at any stage of the loop, the links 1 to Link are being filled until either the attack has been exhausted (Line 9-10) or the water-level reaches that of the next link (Line 11-12).

Yet, the algorithm differs drastically from Algorithm 1 in its form and complexity. Indeed, from Equation 15, all links ∈ D are such that C -δ 's are Algorithm 2: The algorithm of player Π δ which defines an optimal strategy for unknown link usage.

Input: Vector C of channel capacities, of size L Output: Attack vector δ 1 Sort links with decreasing capacity C 2 Attack ← ∆ // Amount of ∆ left to be allocated 3 Link ← 1 // The link index up to which we attack 4 Diff ← 0 // Extra capacity of the current link to the next 5 η ← 0 // Amount to be allocated in each link 6 δ ← 0 // The attack vector 7 while Link ≤ L and Attack > 0 do equal, which amounts to say that for i, j D, we have δ i -δ j = C i -C j . Hence, the different links are being filled at the same rate η, which allows us to simply derive the level of exhaustion of ∆ or when the set D is to be modified. As opposed to Algorithm 1 which computes the solution with arbitrary precision, Algorithm 2 gives the exact solution. Further, the loop runs for at most L times and the solution is obtained after at most O(L) multiplications and O(L 2 ) additions.

8 Diff = C Link -C Link +1 9 if (Link = L or Attack < Link × Diff) then 10 η ← Attack / Link
Figure 2 shows a typical run of the algorithm, with a set of 5 links. There, the algorithm terminates after 3 loops in the "while" command, as δ = ∆.

Optimal Link Usage Strategy with Unknown Degradation Attack

We finally proceed to the case where player Π x chooses its routing strategy without knowledge of the attack performed by player Π δ . Then, we consider its strategy to be the one that has the best delay guarantee, i.e. the one such that the delay it will suffer from is the lowest possible one in the worst case scenario (i.e. where player Π δ has the strongest attack). The problem is referred to as minmax and given below:

C1 C2 C3 C4 C5 C1 C2 C3 C4 C5 C1 C2 C3 C4 C5 (a) (b) (c)
α Fig. 2: Typical run of the water-filling algorithm. First, in (a), the channel with highest capacity is filled so as to reach the level of the second. Then in (b) the two channels of largest capacities are filled. Finally, in (c) as the total attack is exhausted before reaching the level of C 4 , then channels 1 to 3 are equally filled with the remaining attack.

Definition 2. x * is a minmax strategy if it is solution of: mM(C, X, ∆) : min x≥0 max δ≥0 D(x, C -δ), s.t x = X, δ = ∆. (16)
Note that the minmax strategy is also the optimal one in a scenario of a twoplayer sequential game with perfect information where player Π x plays first followed by player Π δ .

Existence and Characterization of the Optimal Strategy

The following theorem states the uniqueness of the solution and gives a characterization:

Theorem 4. There exists a unique x * solution of Eq. (16). It is such that there exists a unique α and λ such that

x * =        C -∆ 2 + ∆ 2α 1 -4α C ∆ + (α -1) 2 if ∈ C M , C -C /λ if ∈ C I , 0 otherwise. (17)
with

C M ∈ L, C ≥ 1 λ ∆ -α∆ ∆ -α/λ 2 C I ∈ L, 1 λ ∆ -α∆ ∆ -α/λ 2 ≥ C ≥ 1/λ (18)
and the set of optimal responses δ * of player Π δ are: Hence, problem ( 16) is equivalent to the following equivalent constrained optimization problem:

{δ|∃ * ∈ C M , ∀ = * ,
min x,α D(x, C) + α s.t.      ∀x , ∆x (C -∆ -x )(C -x ) ≤ α,
x ≥ 0, and x = X.

(

) 19
The corresponding Lagrangian is

L(x, α, λ, µ) = l x C -x + α + λ(X - x ) - µ α + x C -x - x C -∆ -x with ∀ , x ≥ 0, µ ≥ 0.
Let C M be the set of links for which α + x C -x -

x C -∆-x = 0 and x M the corresponding loads. Then, x M satisfies

x M ∆ (C -∆ -x M )(C -x M ) = α i.e. x M = C - ∆ 2 + ∆ 2α 1 -4α C ∆ + (α -1) 2 .
If / ∈ C M , then the Karush Kuhn Tucker conditions give that µ = 0 and hence the lagrangian reduce to Eq. (9). Then, Eq. 10 leads to Eq. 17.

Finally, ∈ C M iff

x M ≤ x I i.e. x M ∆ (C -∆ -x M )(C -x M ) ≥ α But x M ∆ (C -∆ -x M )(C -x M ) = λ∆ √ C -∆ √ λ √ C -∆ √ λ .
Therefore

∈ C M iff C ≥ 1 λ ∆ -α∆ ∆ -α/λ 2 .
Note that α represents the degradation induced by player Π δ . One can readily check that α = 0 ⇔ C M = ∅ which leads to Eq. 8, that is the optimal strategy for Π δ when there is no capacity degradation. Algorithm 3: A water-distributed algorithm for optimal strategy of player Π x with unknown attack 1 Sort links with decreasing capacity C 2 ε ← 0.01, εα ← 0.1 // Set it to the desired accuracy 3 T A ← 0 // The traffic to be redistributed 4 M ← 1 // The link index up to which we reduce the flow 5 I ← 1 // The link index up to which we increase the flow 6 x ← Solution of Algorithm 1 with no attack (δ = 0)

7 α ← x 1 C 1 -∆-x 1 -x 1 C 1 -x 1 8 value ← x 1 C 1 -∆-x 1 + L l=2 x l C l -x l 9 prec ← value + 1 10 while value < prec do 11 α ← α -εα, prec ← value, I ← M 12 for = 1 to M do // Reduce all Links in M 13 T A ← T A + x 14 x ← C -∆ 2 + ∆ 2α 1 -4α C ∆ + (α -
η ← ε(Cj -xj) 2 Cj + ε(Cj -xj) 24 xj ← xj + η 25 T A ← T A -η 26 value ← x 1 C 1 -∆-x 1 + L l=2 x l C l -x l

An Algorithmic Solution

We use the equivalent optimization problem given in Eq. ( 19). Since it is a convex optimization problem, standard optimization tools (e.g. a projected gradient descent on the Lagrangian) can be used, although they exhibit poor performance, in particular because of the nature of the needed projection and the system's size. Therefore, we propose an algorithm, in a similar vein to water-filling algorithms, which we refer to as water-distributed algorithm. It can be seen as a water-filling algorithm with a top cap on the water level (represented by α).

The mechanism is given by Algorithm 3. We initialize it by using Algorithm 1 to compute the optimal allocation x if there was no attack. We deduce the initial value of α.

We then iteratively decrease the value of α and compute the corresponding allocation x. The algorithm ends when no gain in the delays is obtained.

Instead of computing the whole allocation at each iteration of the algorithm, we compute the amount of flow which is removed to the links of C M as a consequence of the decrease of α (lines 21 to 25) and then redistribute this amount to the other links (Line 26 to 36). As long as ∆ > 0, the set C M is always non-empty (as it should include the link of highest capacity). In contrast, the set C I can be empty or not. Further, the set C M ∪ C I may cover all links or not. The different situations are illustrated in the scenario of Figure 3. The system has a set of two links. In Figure 3a both of them are in C M . In this case, the set C I is empty. Figure 3c shows a scenario where C I is also empty and C M consists of only the link of highest capacity. Finally, Figure 3b shows a case where C M only contains the link of higher capacity, while the other one is in C I .

Conclusion

We have studied in this paper a game between a router that has a fixed demand to ship and a malicious controller that affects the system capacity by some fixed amount and can decide how to split this among different links. It turned out to be a non standard zero-sum game since the cost is convex for both the minimizer and maximizer and thus does not have a saddle point. We thus focused on computing the max-min and the min-max value and proposed efficient solution algorithms. While the max-min problem is solved using a water-filling algorithm, the solution of the minmax problem requires a more complex algorithm which we call water-distributing algorithm. We plan in the future to extend the problem to several players that try selfishly to minimize their cost in the presence of adversarial capacity degradation controller.

  Diff 13 Attack ← Attack -η. Link 14 for j = 1 to Link do 15 δj ← δj + η 16 Link ++ 17 return δ

  δ = 0 and δ * = ∆} . Proof. From Theorem 1, we can write x * = arg min max ∈L D(x, C -∆e ) = arg min D(x, C) + max ∈L D (x , C ) -D (x , C -∆) .

1 ) 2 15T

 12 A ← T A -x 16 while T A > 0 do // Redistribute T A among the links 17 while M < L and α ≤ D(x M +1 , C M +1 -∆) -D(x M +1 , C M +1 ) do 18 M + + 19 I ← M 20 while I < L and ND I (x, C) ≥ ND I +1 (x, C) do 21 I + + 22 for j = M + 1 to I do 23

Fig. 3 :

 3 Fig.3: Different sets C M and C I . In the three examples, C 2 = 30, X = 10 and ∆ = 2. The min-max solution is represented with the black point. In each plot, the two graphs (red and blue) represent the overall delay that would be experienced by the user if the attack was concentrated on a single link (1 and 2 respectively).

6. 3

 3 Different C M and C I Note that the set C M and C I both depend on the parameter X, ∆ and the link capacities C , 1 ≤ ≤ L.


Fabrice Mellor

Ardhuin
email: ardhuin@shom.fr

Alastair D Jenkins

Konstadinos A Belibassakis

Fabrice Ardhuin

Commentary on 'The Three-Dimensional Current and Surface Wave Equations' by George

The lowest order sigma-transformed momentum equation given by Mellor (J. Phys. Oceangr. 2003) takes into account a phase-averaged wave forcing based on Airy wave theory. This equation is shown to be generally inconsistent due to inadequate approximations of the wave motion. Indeed the evaluation of the vertical flux of momentum requires an estimation of the pressure p and coordinate transformation function s to first order in parameters that define the large scale evolution of the wave field, such as the bottom slope. Unfortunately there is no analytical expression for p and s at that order. A numerical correction method is thus proposed and verified. Alternative coordinate transforms that allow a separation of wave and mean flow momenta do not suffer from this inconsistency nor require a numerical estimation of the wave forcing. Indeed, the problematic vertical flux is part of the wave momentum flux, thus distinct from the mean flow momentum flux, and not directly relevant to the mean flow evolution.

Introduction

Wave-induced motions are of prime importance in the upper ocean, and in the coastal ocean (e.g. [START_REF] Ardhuin | Waves and operational oceanography: towards a coherent description of the upper ocean for applications[END_REF] for a recent review). Therefore, usual three-dimensional primitive equations must be modified to account for waves. Among such modified equations, those based on surface-following coordinates provide physically sound definitions of velocities right up to the free surface, allowing a proper representation of surface shears and mixing on a vertical scale smaller than the wave height (i.e. a few meters). Any change of coordinate adds some complexity in the derivation, but the final equations can be relatively simple because part of the advective fluxes are removed, and boundary conditions may be simplified. A new set of such equations was recently derived by [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] using a change of the vertical coordinate only, arguably the simplest possible. [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] set of equations was originally derived for monochromatic waves, but it is easily extended to random waves (e.g. Ardhuin et al. 2004, eq. 8). Unfortunately, we show here that these equations, in the form given by Mellor, are not consistent in the simple case of shoaling waves without energy dissipation. A modification is proposed to solve the problem, but it requires a numerical evaluation of the wave forcing terms. This difficulty is due to the choice of averaging, and the same problem arises with the alternative Generalized Lagrangian Mean equations of Andrews and McIntyre (1978, eq. 8.7a, hereinafter aGLM). Both Mellor's and aGLM equations describe the evolution of a momentum quantity that contains the three-dimensional wave (pseudo)-momentum (hereinafter called 'wave momentum' for simplicity, see McIntyre 1981 for details). Writing an evolution equation for this quantity requires an explicit description of the complex vertical fluxes of wave momentum that are necessary to maintain the vertical structure of the wave field in the surface gravity waveguide.

The problem: wave motions and wave-following vertical coordinates

We discuss here the simple case of monochromatic waves of amplitude a and wavenumber k propagating in the horizontal x direction, with all quantities uniform in the other horizontal y direction. The surface and bottom elevations are η(x, t) and -h(x), respectively, so that the local mean water depth is D(x, t) = h(x) + η(x, t), with the overbar denoting an Eulerian average over the wave phase. We shall assume that the maximum surface slope is a small parameter ε 1 = ka ≪ 1, and that the Eulerian mean current u in the x-direction is uniform over the VOLUME depth. Thus ω will denote the radian wave frequency related to k by the linear wave dispersion relation (e.g. [START_REF] Mei | Applied dynamics of ocean surface waves[END_REF],

ω = ku + σ = ku + [gk tanh(kD)] 1/2 .
(1)

Finally, we assume that the water depth, current and wave amplitude change slowly along the x-axis with a slowness measured by a second small parameter ε 2 taken to be the maximum bottom slope. We thus assume

|(∂D/∂x)| ≤ ε 2 , |(∂a/∂x)| ≤ ε 2 , |(∂u/∂x)/(σ)| ≤ ε 2 , |(∂a/∂t)/(σa)| ≤ ε 2 , |k(∂u/∂t)/σ 2 | ≤ ε 2 , |(∂D/∂t)k/σ| ≤ ε 2 .
The conditions on the bottom slope and current gradients are consistent with the condition on the wave amplitude gradient because in steady conditions the wave amplitude would change due to shoaling over the current and/or bottom.

The vertical coordinate z is implicitly transformed into Mellor's ς coordinate through

z = s (x, ς, t) = η + ςD + s (2)
with s defined by Mellor's eq. ( 23b) as

s = s 0 = aF SS cos (kx -ωt) (3)
and the vertical profile function F SS defined by

F SS = sinh [kD( 1 + ς )] sinh (kD) = sinh [k( z + h )] sinh (kD) + O( a D
).

(4) The coordinate transformation from z to ς has the very nice property of following the vertical wave-induced motion, at least for linear waves on a flat bottom, and to first order in ε 1 . In that case the iso-ς surfaces are material surfaces, and the fluxes of horizontal momentum through one of these surfaces are simply correlations of pressure p times the slope of that surface ∂s/∂x (figure 1.c), which replaces the wave-induced advective flux uw in an Eulerian point of view (figure 1.a). More generally, when averaging is performed following water particles over their trajectory (Lagrangian) or over their vertical displacement (Mellor-sigma), the corresponding advective flux of momentum u i u j is replaced by a modified pressure force (figure 1). 1Using his coordinate transform, [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] obtained a phase-averaged equation for the drift current U = u + u S where u S is the Stokes drift, i.e. the mean velocity of water particles induced by fast wave-induced motions. U is strictly defined as the phase-average particle drift velocity when following the up-and-down wave motion, and u = U -u S is a quasi-Eulerian mean current [START_REF] Jenkins | A theory for steady and variable wind-and waveinduced currents[END_REF](Jenkins , 1987)). Below the wave crests u is equal, to second order in the wave slope, with the Eulerian mean current u (figure 2). Mellor's horizontal mean momentum equation (34a) is reproduced here for completeness, in our conditions with a flow restricted to the vertical x, z plane, a constant water density, no Coriolis force, and no turbulent fluxes and the atmospheric mean pressure set to zero (wind-wave generation due to air pressure fluctuations is absorbed in F x3 ), [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF]. The thick black bars connect the fixed points (x, z) where the average field is evaluated, to the displaced points (x, z) + (ξ 1 , ξ 3 ) where the instantaneous field is evaluated. For averages in moving coordinates the points (x, z) + (ξ 1 , ξ 3 ) at a given vertical level ξ are along the gray lines. The drift velocity is the sum of the (quasi-Eulerian) current and the wave-induced mass transport. In the present illustration an Airy wave of amplitude 3 m and wavelength 100 m in 30 m depth, is superimposed on a hypothetical current of velocity u(z) = -0.5 -0.01z m/s for all z < ζ(x). The quasi-Eulerian current profile is not represented in (c) since it is not directly given in Mellor's theory, although it can obviously be obtained by taking the difference of the other two profiles.

∂DU ∂t + ∂DU 2 ∂x + ∂ΩU ∂ς + gD ∂ η ∂x = F xx + F x3 . (5)
On the right-hand the first other term

F xx = - ∂S xx ∂x = - ∂ ∂x D u 2 + p ∂ s ∂ς (6)
represents the convergence of a horizontal flux of horizontal momentum that accelerates the mean drift velocity U .

The other term

F x3 = - ∂S x3 ∂ς = ∂ ∂ς p∂ s/∂x (7)
represents a similar convergence of a vertical flux of horizontal momentum. Defining g as the acceleration due to the apparent gravity, p and s are of the order of ga and a respectively. In general p and s are almost in phase, thus the flux S x3 is of the order of ga ∂a/∂x, and the force F x3 is of the order of gDka ∂a/∂x. Thus, in the case of shoaling waves, F x3 is of the order of gDε2 1 ε 2 . Mellor estimated the vertical momentum flux S x3 from (3) and the corresponding lowest order wave-induced kinematic pressure on ς levels 2 ,

p = p 0 = ga (F CC -F SS ) cos(kx -ωt), (8)
g is the acceleration due to the apparent gravity, and the vertical profile function F CC is defined by

F CC = cosh [kD( 1 + ς )] cosh (kD) = cosh [k( z + h )] cosh (kD) +O( a D
).

(9) For non-dissipating shoaling waves, the right hand side terms of eq. ( 5) are of order gDε 2 1 ε 2 . The estimation of F x3 thus requires the knowledge of p/(gD) and k s to order ε 1 ε 2 , for which Airy theory is insufficient. In particular, this estimation demands a formal definition of s, not given by [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF]. Further, eq. ( 7) is only valid if the wave-induced velocity Ω through ς levels is zero, or at least, yields a negligible flux Ω u and a negligible mean Jacobian-weighted vertical velocity Ω = Ω/(1 + ∂s/∂ς/D). This is not the case over a sloping bottom with [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] s function.

a. Formal definition of the the coordinate change s

For a general surface ς defined implicitly by z = s(x, ς, t), the ς velocity component Ω is (e.g. Mellor 2003 eq. 20),

Ω = d(z -s) dt - d(z -s) dt = w -u ∂s ∂x -u ∂ s ∂x - ∂ s ∂t , (10)
with s = η + ςD.

In the spirit of [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] derivation, the ς levels should be material surfaces for wave-only motions, so that one may neglect the vertical flux of momentum

(U + u) Ω/(1 + ∂s/∂ς/D).
Using the wave-induced vertical and horizontal displacements, ξ 3 (x, ς, t) and ξ 1 (x, ς, t), defined by

∂ξ i /∂t = u i (x 1 + ξ 1 , z + ξ 3 , t), we redefine the wave part of s, s ′ ≡ ξ 3 -ξ 1 ∂s ∂x . ( 11
)
The first term ξ 3 corresponds to Mellor's definition while the second is a O(ε 2 ) relative correction. This definition yields a wave-induced vertical velocity Ω = -u∂ s ′ /∂x through the iso-ς surfaces redefined by z

= s (x, ς, t) = η + ςD + s ′ . If u ≪ u, as in the examples below, then
Ω is of a higher order compared to that given by [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] s (eq. 2).

b. Wave-induced vertical displacements and pressure over a sloping bottom

A WKBJ approximation using Airy's theory is sufficient for estimating ∂F xx /∂x because the horizontal gradient of any wave-averaged quantity φ is of order ε 2 kφ. On the contrary, the other force F x3 is affected by modifications s ′ 1 and p 1 to the local-flat-bottom solutions s ′ , and p.

For small bottom slopes, s ′ 1 and p 1 are expected to be of the order of ε 2 s ′ and ε 2 p, i.e. of order aε 2 and gaε 1 ε 2 , respectively. Thus ∂ s ′ 1 /∂x is of order kaε 2 = ε 1 ε 2 , and is expected to be in phase with the wave-induced pressure (8), of order ga, giving another term of order gDε 2 1 ε 2 omitted by Mellor in his estimation of ∂( p∂ s/∂x)/∂ς. The modification of the pressure can be obtained from the modification of the velocity potential, and it may be in phase with ∂ s 0 /∂x, thus also contributing at the same order to F x3 .

In order to be convinced of the problem, one may consider the case of steady monochromatic shoaling waves over a slope without bottom friction, viscosity or any kind of surface stress. We also neglect the Coriolis force. In this mathematical experiment, the flow is purely irrotational. We consider that the non-dimensional depth kH is of order 1, and that there is no net mass flux across any vertical section. In that case the mean current and the Stokes drift are of the same order, i.e. of the order Cε 2 1 with C the phase speed. The mean current exactly compensates the divergence of the wave-induced mass transport, and the mean sea level is lower in the area where the wave height is increased (Longuet-Higgins 1967)

η(x) = - kE sinh(2kD) + k 0 E 0 sinh(2k 0 D 0 ) (12)
where the 0 subscript correspond to quantities evaluated at the offshore boundary of the domain. Since wave forcing is steady, the Eulerian mean current response is steady (e.g. Rivero and Sanchez-Arcilla 1994[START_REF] Mcwilliams | An asymptotic theory for the interaction of waves and currents in coastal waters[END_REF], Lane et al. 2006), and thus the Lagrangian mean current is also steady. Thus the first term in ( 5) is zero and the second is of order

DC 2 ε 4 1 ε 2 /D ≃ gDε 4 1 ε 2 .
The vertical mean velocity Ω can be estimated from the steady mass conservation equation, ∂DU ∂x

+ ∂Ω ∂ς = 0 (13)
where the first term is of order DCε 2 1 ε 2 /D and the second is of order Ω. Thus the third term in (5) is of order

C 2 ε 4 1 ε 2 ≃ gDε 4 1 ε 2 .
The remaining terms in (5) are of order ε 2 1 ε 2 , giving the lowest order momentum balance

Feta -D ∂ ∂x (gη) Fxx - ∂S xx ∂x + Fx3 ∂ ∂ς p∂ s/∂x= 0. ( 14
)
For reference the corresponding lowest order Eulerian mean balance is (e.g. [START_REF] Rivero | On the vertical distribution of u w[END_REF]Sanchez-Arcilla 1994, Lane et al. 2006)

- ∂ ∂x gη -w 2 - ∂ u 2 ∂x - ∂ u w ∂z = 0, . (15)
Only the hydrostatic pressure gradient is present in both the Eulerian and Mellor-sigma balances, because the other terms represent a different including wave momentum in the latter (see figure 2). Equation ( 14) is now tested numerically. We take a Roseau-type bottom profile (1976) defined by x and z coordinates given by the real and imaginary part of the complex function

Z(x ′ ) = h 1 (x ′ -iα) + (h 2 -h 1 ) ln(1 + e x ′ -iα ) α .
(16) With α = 15π/180, h 1 = 6 m and h 2 = 4 m (figure 1), and a radian frequency ω = 1.2 rad s -1 (i.e. a frequency f = 0.2 Hz), the non-dimensional water depth varies between 0.85 < kH < 1.1. The reflection coefficient for the wave amplitude is 1.4 × 10 -9 [START_REF] Roseau | Asymptotic wave theory[END_REF], so that reflected waves may be neglected in the momentum balance. We illustrate the force balance obtained for waves with an offshore amplitude a 0 = 0.12 m, which corresponds to a maximum steepness ε 1 = ka = 2.6 × 10 -2 equal to the maximum bottom slope ε 2 = ε 1 . The change in wave amplitude is given by the conservation of the wave energy flux (see [START_REF] Ardhuin | On the momentum balance in shoaling gravity waves: a commentary of shoaling surface gravity waves cause a force and a torque on the bottom by K. E. Kenyon[END_REF] for a thorough discussion), and the wave phase ψ is taken as the integral over x of the local wavenumber, so that ∂ψ/∂x = k. The various terms are then estimated using second order finite differences on a regular grid in ς coordinates, with 201 by 401 points covering the domain shown in figure 3.a. The three terms in eq. ( 14) are shown in figure 3.

We have verified that the depth-integrated forces are in balance, within 0.1% of F eta . However, at most water depths there is a large imbalance, of the order of the individual forces ( i.e. gDε 2 1 ε 2 ) ,up to 180% of F eta . This contradicts the known steady balance obtained from the Eulerian-mean analysis of Rivero and Sanchez-Arcilla (1994).

For the case of shoaling waves without breaking the three-dimensional equations of motion of [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] are not consistent to their dominant order, because of an improper approximation of S x3 . This conclusion holds for any relative magnitude of the wave and bottom slopes ε 1 and ε 2 .

c. Wind-forced waves

Clearly, any deviation of the wave-induced fields s ′ , p, and u from Airy-wave theory may have strong effects on the vertical momentum flux term S x3 . Another example of such a situation, correctly described by Mellor, is the case of wind-wave generation. We briefly address it here because the full solution can has not been given previously. Mellor focused on the wind-wave generation contribution to the vertical momentum flux term p∂ s/∂x term. This equals the wave-supported wind stress at the sea surface, and, below, it explains the growth of the wave momentum profile with the same profile as that of the Stokes drift [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF].

In horizontally uniform conditions, the wave amplitude is a function of time only, and for the sake of simplicity we shall solve the problem in the frame of reference moving at the velocity at which the wave phase is advected by the current. We write the wave-induced non-hydrostatic kinematic Eulerian pressure in the form p E = p E0 + p Ew , the elevation as η = ζ 0 + ζ w and the velocity potential as φ = φ 0 + φ w , in which the 0 subscript refers to the primary waves, and the w subscript refers to the added components in the presence of wind forcing. Taking a primary surface elevation of the form ζ 0 = a cos ψ with the the phase ψ = kx -σt, Mellor considered an atmospheric kinematic pressure fluctuation in quadrature with the primary waves

p a = -gβ ρ w ρ a a sin ψ, (17)
with β a small non-dimensional wave growth factor, and ρ w and ρ a the densities of water and air respectively. He then assumed that the water-side wave-induced pressure was of the form

p Mellor = -gβa cosh [k( z + h )] cosh (kD) sin ψ. ( 18
)
Implicitly s ′ w is zero, and for his purpose φ w was irrelevant. We shall now also determine φ w . The continuity of dynamic pressures at the surface is3

p Ew + g s ′ w = -gβa sin ψ at z = ζ. ( 19
)
A solution is obtained by solving Laplace's equation with proper boundary conditions, to first order in β. The boundary conditions include the Bernoulli equation,

∂φ ∂t = -gζ - ρ a ρ w p a , at z = ζ, (20)
in which non-linear terms have been neglected because they are the sum of products of the form ∇φ 0 • ∇ζ 0 , unchanged from the case without wind, and terms of the form ∇φ w • ∇ζ 0 , which are negligible compared to the left-hand side terms for primary waves of small slope.

Similarly, the surface kinematic boundary condition is linearized as

∂φ ∂z = ∂ζ ∂t at z = ζ. (21)
The combination of both yields

∂ 2 φ ∂t 2 + g ∂φ ∂z = - ρ a ρ w ∂ p a ∂t at z = ζ. ( 22
)
φ w is also a solution of Laplace's equation with the bottom boundary condition ∂φ w /∂z = 0 at z = -h. With the fully resonant atmospheric pressure (17) envisaged by Mellor, one has

ζ 0 = a(t) cos ψ, (23)
φ 0 = ga(t) σ F CC sin ψ, (24)
p E0 = ga(t)F CC cos ψ, (25) da(t) dt = βσa(t) 2 , ( 26
)
φ w = βg a 2σ F CC [A cos ψ + B cos ψ ′ ] ( 27
)
with ψ ′ = (kx + σt). The elevation and under-water non-hydrostatic pressure corresponding to φ w are given by ( 21) and the linearized Bernoulli equation

∂ ∂t (φ 0 + φ w ) = -p E0 -p Ew (28) VOLUME yielding ζ w = β a 2 [(1 -A) sin ψ + B sin ψ ′ ] (29) p Ew = gβ a 2 F CC [-(1 + A) sin ψ + B sin ψ ′ ] . (30)
Mellor's expression for p w (eq. 18) is obtained by replacing ζ w and p Ew in ( 19), giving A = 1. One may take B = 1 to have φ w = 0 at t = 0, or more simply B = 0, which gives ζ w = 0, and p Ew = F CC p a ρ a /ρ w . The choice of B has no dynamical effect. In the present case φ w should give a contribution to S xx because it is in phase with φ 0 , but this is a relative correction of order β, thus negligible. To the contrary, the contribution of p w to ( p∂ s/∂x) is quite important, because for uniform horizontal conditions this flux is otherwise zero.

A solution to the problem ?

Contrary to that particular wind-forcing term, there is no simple asymptotically analytical correction for p and s ′ that can account for the bottom slope and wave field gradient. A major problem in this situation is that the wave velocity potential becomes a non-local function of the water depth. The velocity potential and pressure fields may only be investigated analytically for plane beds (e.g. [START_REF] Ehrenmark | An alternative dispersion equation for water waves over an inclined bed[END_REF] or specific bottom profiles. Numerical solutions for the three-dimensional wave motion are generally found as infinite series of modes (e.g. Massel 1993). The velocity potential for any of these modes satisfies Laplace's equation with a local vertical profile F n proportional to cos(k n z + k n h) and a dispersion relation σ 2 = gk n tan(k n D). The local amplitudes of these modes are non-local functions of the water depth, and may be obtained numerically with a coupled-mode model [START_REF] Massel | Extended refraction-diffraction equation for surface waves[END_REF]. This non-local dependance of the wave amplitude on the water depth arises from the elliptic nature of Laplace's equation, satisfied by the velocity potential in irrotational conditions. The series of modes can be made to converge faster by adding a 'sloping bottom mode' that often accounts for a large part of the correction and is a local function of the depth and bottom slope. It is thus of interest to see if that correction only, without the infinite series, may provide a first order analytical correction to Mellor's momentum flux S x3 .

Following [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small amplitude water waves over variable bathymetry regions[END_REF], one may define the velocity potential for that mode as

φ 1 = - dh dx aσDF (z) cos ψ, (31)
In order to satisfy the bottom boundary condition w = ∂φ 1 /∂z = -dh dx ∂φ 0 /∂x, the function F should verify DdF (-h)/dz = 1/ sinh(kH) and the satisfaction of the surface boundary condition may be obtained with F (0) = dF (0)/dz = 0. Athanassoulis and Belibassakis (1999) have used

F = F AB ≡ 1 sinh(kD) z -ζ D 3 + z -ζ D 2 , (32)
and [START_REF] Chandrasekera | Linear refractiondiffraction model for steep bathymetry[END_REF] have used

F = F Ch ≡ 1 kD sinh 2 (kD)
1 -cosh(kz -kζ) .

(33) With these choices φ 1 does not satisfy exactly Laplace's equation, and thus requires further corrections in the form of evanescent modes. An infinite number of other choices is available, either satisfying Laplace's equation or the surface boundary conditions, but never both, so that each of these solutions is only approximate, and the exact solution is given by the infinite series of modes, which can be computed numerically for any bottom topography (e.g. [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small amplitude water waves over variable bathymetry regions[END_REF][START_REF] Belibassakis | A coupled-mode model for the refraction-diffraction of linear waves over steep three-dimensional bathymetry[END_REF], Magne et al. 2006).

The vertical displacement and Eulerian pressure corrections are given by time integration of the vertical velocity and the linearized Bernoulli equation,

ξ 31 = dh dx aD dF dz sin ψ (34) p E1 = dh dx aDF (z) sin ψ (35)
Thus, in absence of wind forcing but taking into account the 'sloping bottom mode' to first order in the bottom slope, the wave-induced flux of momentum through iso-ς surfaces is

p ∂ s ′ ∂x = (F CC -F SS ) ga 2 ∂(aF SS ) ∂x + gka 2 2 dh dx D dF dz + ςF CS + gka 2 2 dh dx -F SC F + F SS D dF dz + ςF CS , (36)
with F CS = cosh [kD( 1 + ς )]/sinh (kD). The first line is the term given by [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF]. The second line arises from the correction due to the difference between s ′ and s, and the third line arises due to corrections p 1 = p E1 -gs ′ 1 to the pressure on ς levels. These additional term are of the same order as the first term, and have no flux at the bottom and surface. Thus the depthintegrated equations including that term also comply with known depth-integrated equations (e.g. [START_REF] Smith | Observed variability of ocean wave Stokes drift, and the Eulerian response to passing groups[END_REF]).

In the case chosen here F Ch gives a net momentum balance closer to zero than [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] original expression (figure 4). However, the remaining error is significant. Thus one cannot use only that mode, and the contribution of the evanescent modes have to be computed, which can only be done numerically. 14) for steady shoaling waves over a smooth bottom profile. The net force has been integrated over x and normalized by a similar integration of the the hydrostatic pressure force Feta. Several solutions are obtained. One corresponds to Mellor's orginal expression, one possible analytical correction using F Ch , and numerical estimations using the NTUA model, with various numbers of modes.

A numerical evaluation of the forces was performed using the NTUA model [START_REF] Athanassoulis | A consistent coupled-mode theory for the propagation of small amplitude water waves over variable bathymetry regions[END_REF]. The NTUA solution was obtained in a domain with 401 points in the horizontal dimension. For the small bottom slope used here, the model contains a numerical reflection R = 0.002 much larger than the analytical value given by [START_REF] Roseau | Asymptotic wave theory[END_REF]. However, this only introduces a modulation, in the x direction, of the estimated forces. This modulation is significant but still relatively smaller than the average. The net force estimated from NTUA results is found to converge to the expected force balance described by eq. ( 14) as the number of evanescent modes is increased (figure 4). In this calculation the values of F xx do not differ significantly from those estimated using Mellor's analytical expressions, as expected. The only significant difference between the NTUA numerical result with 10 modes, and Mellor's analytical expression is found in F x3 , with a much stronger value near the surface in the numerical result, allowing a balance with the strongly sheared F xx (figure 4). [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] changed the vertical coordinate from ς to z, using an implicit function s in two parts, z = s(x, y, ς, t) + s(x, y, ς, t) with s changing only slowly in space and time and s representing the faster waveinduced change of vertical coordinate. If the ς levels are material surfaces, then the momentum flux S x3 = p ξ ∂ s/∂x is the surface-following coordinate counterpart of the Eulerian vertical momentum flux term u w discussed by [START_REF] Rivero | On the vertical distribution of u w[END_REF], with p ξ the wave-induced pressure at the displaced position (in the surfacefollowing coordinates). However, p ξ ∂ s ′ /∂x and u w do not represent the same physical quantity since the former contains wave momentum, which is not included in the latter.

Conclusions

Just like the Eulerian momentum flux u w is modified by the bottom slope, wave amplitude gradients, windwave generation, boundary layers, or vertical current shears, these effects also modify S x3 . But in these situations, the ς levels as defined by [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF] are not material surfaces, and a missing Eulerian-like flux term Ω w would have to be added to correct the momentum equations, with Ω the wave-induced velocity across ς levels. Alternatively, we propose to replace s with s ′ , defined by eq. ( 11) such that ς levels are closer to material surfaces, i.e. so that Ω is of a higher order.

Whether the original s or our corrected s ′ is used, the wave-induced momentum flux S x3 must be estimated to first order in the bottom slope ε 2 for consistency. This requires an O(ε 2 ) estimation of both p ξ and s or s ′ . Unfortunately there is no analytical O(ε 2 ) expression for the wave motion. Thus Mellor's equations, even when corrected, require a computer-intensive solution that is generally not feasible. For example, Magne et al. (2006) only included a total of five modes in their calculation of wave propagation over a submarine canyon. In an example shown here, this small number of modes is insufficient for an accurate estimation of wave-forcing terms.

The trouble with these equations can be avoided by using, instead, equations of motion for the quasi-Eulerian velocity u = U -u S [START_REF] Jenkins | A theory for steady and variable wind-and waveinduced currents[END_REF](Jenkins , 1987(Jenkins , 1989)). Such equations have been obtained in the limit of vanishing wave amplitude using an analytical continuation (e.g. using a Taylor expansion) of the current profile across the surface [START_REF] Mcwilliams | An asymptotic theory for the interaction of waves and currents in coastal waters[END_REF]. A general and explicit solution can also be obtained from the exact Generalized Lagrangian Mean (GLM) equations of Andrews and McIntyre (1978a) expanded to second order in the surface slope ε 1 (Ardhuin et al., manuscript submitted to Ocean Modelling, see [http://arxiv.org/abs/physics/0702067]). In these, the equation for the horizontal quasi-Eulerian momentum involves no flux term like p∂ s/∂x because this corresponds to the flux p ξ ∂ξ 3 /∂x[1 + O(ε)] of wave momentum u S (Andrews and McIntyre 1978b, eq. 2.7b), not directly relevant to the problem of mean flow evolution (see also [START_REF] Jenkins | Interaction of ocean waves and currents: How different approaches may be reconciled[END_REF]. This flux of wave momentum only appears in evolution equations for the total momentum U , such as given by [START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF], or the 'alternative' form of the GLM equations (Andrews and McIntyre 1978, eq. 8.7a).

For that reason, the equations for the quasi-Eulerian velocity u are simple and consistent in their adiabatic form (without wave dissipation), at least to lowest order in wave slope and current vertical shear, for which analytical expressions exist for the wave forcing terms. Further details on the relationships between all these equa-VOLUME tions, and further validation against numerical solutions of Laplace's equation can be found in Ardhuin et al. (submitted manuscript).

Acknowledgments. This work was sparked by a remark of the late and much missed Tony Elfouhaily. The open collaboration of George Mellor and his criticisms of an early draft greatly improved the present note, together with the feedback of anonymous reviewers and many colleagues. Among them Dano Roelvink suggested the evaluation of the torque resulting from Mellor's equations in the vertical plane, and Jacco Groeneweg made a detailed critique of early drafts. Initiated during the visit of A.D.J. at SHOM, this research was supported by the Aurora Mobility Programme for Research Collaboration between France and Norway, funded by the Research Council of Norway (NFR) and The French Ministry of Foreign Affairs. A. D. J. was supported by NFR Project 155923/700.

  FIG. 1. Wave-induced fluxes of horizontal momentum in Eulerian, Generalized Lagrangian, and Mellor-sigma averages of the flow. Viscous or turbulent fluxes are neglected for simplicity. Distorted squares represent an elementary fluid volume and its position at four phases of the wave cycle, and the large arrow indicate the local wave orbital velocity. The horizontal and vertical fluxes of the horizontal momentum are represented by smaller arrows. Their expression are given to lowest order, without Jacobian corrections due to a change of volume (e.g. this results in u 2 + p becoming Mellor's u 2 + p∂s/∂ς).

  FIG.2. Averaging procedures (left) and examples of resulting velocity profiles (right) in the case of (a) Eulerian averages, (b) the Generalized Lagrangian Mean, (c) sigma transforms[START_REF] Mellor | The three-dimensional current and surface wave equations[END_REF]. The thick black bars connect the fixed points (x, z) where the average field is evaluated, to the displaced points (x, z) + (ξ 1 , ξ 3 ) where the instantaneous field is evaluated. For averages in moving coordinates the points (x, z) + (ξ 1 , ξ 3 ) at a given vertical level ξ are along the gray lines. The drift velocity is the sum of the (quasi-Eulerian) current and the wave-induced mass transport. In the present illustration an Airy wave of amplitude 3 m and wavelength 100 m in 30 m depth, is superimposed on a hypothetical current of velocity u(z) = -0.5 -0.01z m/s for all z < ζ(x). The quasi-Eulerian current profile is not represented in (c) since it is not directly given in Mellor's theory, although it can obviously be obtained by taking the difference of the other two profiles.

  FIG. 3. (a)Snapshot of the pressure field for a slowly varying Airy wave over a the bottom topography given by eq. (16). The forces in the balance (14) are shown in panels b, c and d, with their sum in panel e, all estimated from Mellor's analytical expressions. All forces have been normalized by gDε 3 . N.B. in the case shown here ε 1 = ε 2 = ε.

  FIG. 4. Net forces in the momentum balance (14) for steady shoaling waves over a smooth bottom profile. The net force has been integrated over x and normalized by a similar integration of the the hydrostatic pressure force Feta. Several solutions are obtained. One corresponds to Mellor's orginal expression, one possible analytical correction using F Ch , and numerical estimations using the NTUA model, with various numbers of modes.

For the Generalized Lagrangian Mean (GLM) only the contributions to lowest order in ε 1 are indicated. Indeed, in GLM the waveinduced advective flux is not strictly zero, but of higher order, since the average only follows a zero-mean displacement with a residual advection, contrary to a truly Lagrangian mean with zero advection (e.g.[START_REF] Jenkins | A theory for steady and variable wind-and waveinduced currents[END_REF]).

This pressure includes a hydrostatic correction due to the vertical displacement.

Here the pressure is Eulerian. For correspondance to Mellor's pressures on ς levels, one should take p = p E -g s ′ .