Didier Clamond
email: diderc@unice.fr

Denys Dutykh
email: denys.dutykh@univ-savoie.fr

Dimitrios Mitsotakis

Conservative modified Serre-Green-Naghdi equations with improved dispersion characteristics

Keywords: Shallow water waves, improved dispersion, energy conservation MSC: [2010]74J15 (primary), 74S10, 74J30 (secondary) PACS: [2010]47.35.Bb (primary), 47.35.Fg, 47.85.Dh (secondary)

For surface gravity waves propagating in shallow water, we propose a variant of the fully nonlinear Serre-Green-Naghdi equations involving a free parameter that can be chosen to improve the dispersion properties. The novelty here consists in the fact that the new model conserves the energy, contrary to other modified Serre's equations found in the literature. Numerical comparisons with the Euler equations show that the new model is substantially more accurate than the classical Serre equations, specially for long time simulations and for large amplitudes.

Introduction

Water waves in channels and oceans are usually described by the Euler equations. Due to their complexity, several approximate models have been derived in various wave regimes. Once a new mathematical model is proposed, the limits of its applicability have to be determined. In shallow water, the main restriction comes from the ratio between the characteristic wavelength λ and the mean water depth d, the so-called shallowness parameter σ = d/λ ≪ 1. Restrictions on the free surface elevation are characterised by the dimensionless parameter ε = a/d, where a is a typical amplitude (εσ = a/λ is a wave steepness). Many approximate equations have been derived for waves in shallow water, such as the Korteweg-de Vries (KdV) equation [START_REF] Korteweg | On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves[END_REF] for unidirectional waves, the Saint-Venant equations [39] for bidirectional non-dispersive waves and many variants of the Boussinesq equations [START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media. I: Derivation and linear theory[END_REF]6] for dispersive waves propagating in both directions. In addition to shallowness (σ ≪ 1), KdV and Boussinesq equations assume small amplitudes (e.g. ε = O(σ 2 )) [START_REF] Johnson | A Modern Introduction to the Mathematical Theory of Water Waves[END_REF].

Considering long waves propagating in shallow water but without assuming small amplitudes (i.e. σ ≪ 1 and ε = O(1)), Serre [35] derived a so-called fully-nonlinear weakly dispersive system of equations [40,41] which, after further approximations, include the Korteweg-de Vries, Saint-Venant and Boussinesq equations as special cases. For steady flows, these equations were already known to Rayleigh [26]. The Serre equations were independently rediscovered by Su and Gardner [START_REF] Su | KdV equation and generalizations. Part III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF], and again later by Green, Laws and Naghdi [START_REF] Green | On the theory of water waves[END_REF]. These nowadays popular equations constitute an asymptotic fully nonlinear model including all the terms up to order σ 3 into the momentum equation. Serre's equations represent a substantial improvement with respect to the Boussinesq theory [START_REF] Castro-Orgaz | Boussinesq-and Serre-type models with improved linear dispersion characteristics: applications[END_REF], but many shallow water phenomena involve significant dispersive effects that are not well described by Serre's equations.

One possibility to improve the Serre model is by including the terms of order σ 5 (and higher-order terms) into the momentum equation. This program was accomplished for Boussinesq equations in, e.g., [28]. However, this modification makes the fifth-order (and higher-order) derivatives to appear in the model equations, making its numerical resolution (and thus its applicability) rather challenging. Actually, the numerical resolution of these high-order Boussinesq-like equations is slower (and less accurate) than that of the Euler equations, at least for simple (periodic) domains.

Another way of improving the classical model was coined by Bona and Smith [5] (see also [START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media. I: Derivation and linear theory[END_REF][START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF]). The idea consists in introducing a free parameter into the model that can be appropriately chosen to improve some of the desired properties. This can be achieved by replacing the depth-averaged velocity variable by the velocity of the fluid evaluated at at a certain depth in the bulk of fluid [29]. Most often, this parameter is chosen to optimise somehow the dispersion relation properties [START_REF] Kim | Extended Boussinesq equations for rapidly varying topography[END_REF]27]. In general, the community is rather focused on various linear properties of the model, even if it is used later to simulate nonlinear waves. Similar ideas have also been applied to the Serre equations with flat [START_REF] Dias | On the fully-nonlinear shallow-water generalized Serre equations[END_REF][START_REF] Dutykh | Adaptive modeling of shallow fully nonlinear gravity waves[END_REF][START_REF] Liu | Two sets of higher-order Boussinesq-type equations for water waves[END_REF] and varying bottoms [START_REF] Antunes | Boussinesq and Serre type models with improved linear dispersion characteristics: Applications[END_REF][START_REF] Antunes | Extended Serre Equations for Applications in Intermediate Water Depths[END_REF][START_REF] Castro-Orgaz | Boussinesq-and Serre-type models with improved linear dispersion characteristics: applications[END_REF]. Free parameters are obtained from arbitrarily-weighted averages of different (but of same order) approximations of some quantities. However, these modified Serre equations invalidate one fundamental physical property: the conservation of the energy. Some variants also invalidate the Galilean invariance [START_REF] Duran | On the Galilean Invariance of Some Nonlinear Dispersive Wave Equations[END_REF], that is a severe drawback for many applications. The same remarks apply as well to previous works on the improved Boussinesq equations [3, 30]. Thus, one may improve the dispersive properties of the model but, on the other hand, loses the energy conservation property. For many applications, specially in the case of long time simulations, the disadvantages can be crucial, overriding all the possible advantages. In the present paper, we address this issue, proposing a method for deriving an improved version of the Serre equations that preserves the aforementioned nice properties of the original Serre model. For the sake of simplicity, the method is illustrated for 2D waves over a horizontal bottom, but generalisations to 3D and varying bottom can be obtained in an analogous manner.

The present paper is organised as follows. In Section 2 a simple derivation of the classical Serre equations is presented. It is followed by the derivation of an already known oneparameter generalisation of these equations and its shortcomings are explained. In Section 3 we derive a new one-parameter generalisation of the Serre equations that conserve the energy. In Section 4 we discuss several criteria for choosing the free parameter. Some numerical results are provided in Section 5 demonstrating the advantages of the new Serrelike equations. Finally, the main conclusions, possible generalisations and perspectives are discussed in Section 6.

Classical and modified Serre's equations

Here, we derive the classical Serre equations and a modified version of these equations with an additional free parameter using variational methods.

Classical Serre's equations

In order to model irrotational two-dimensional long waves propagating in shallow water over horizontal bottom, one can approximate the velocity field by the ansatz

u(x, y, t) ≈ ū(x, t), v(x, y, t) ≈ -(y + d) ūx , (2.1)
where d is the water depth and ū is the horizontal velocity averaged over the water column i.e., ū = h -1 ´η -d u dy, h = η + d the total water depthy = η and y = 0 being the equations of the free surface and of the still water level, respectively; a sketch of the fluid domain is depicted in the Figure 1. The horizontal velocity u is thus (approximately) uniform along the water column and the vertical velocity v is chosen so that the fluid incompressibility is valid. Note that the vorticity ω = v x -u y ≈ -(y + d) ūxx is not exactly zero, meaning that the potential flow is approximated by a rotational velocity field. The fact that the irrotationality is violated should not be more surprising than the violation of other relations, such as the isobarity of the free surface.

With the ansatz (2.1), the vertical acceleration is

Dv Dt = ∂ v ∂t + u ∂ v ∂x + v ∂ v ∂y ≈ -v ūx -(y + d) Dū x Dt = γ y + d h ,
where γ is the vertical acceleration at the free surface, i.e.,

γ = Dv Dt y=η ≈ h ū 2 x -ūxt -ū ūxx = 2 h ū 2 x -h ∂ x [ ūt + ū ūx ] . (2.
2)

The kinetic and potential energies per water column, respectively K and V , are similarly easily derived

K ρ = ˆη -d u 2 + v 2 2 dy ≈ h ū2 2 + h 3 ū 2 x 6 , V ρ = ˆη -d g (y + d) dy = g h 2 2 ,
where ρ > 0 is the (constant) fluid density and g > 0 is the acceleration due to gravity (directed downward). A Lagrangian density L can then be introduced as the kinetic minus the potential energies plus a constraint for the mass conservation

L ρ = h ū2 2 + h 3 ū 2 x 6 - g h 2 2 + { h t + [ h ū ] x } φ, (2.3)
where φ is a Lagrange multiplier. Physically, φ is a (approximated) velocity potential expressed at the free surface. The Euler-Lagrange equations for the action functional S = L dx dt are

δ φ : 0 = h t + [ h ū ] x , (2.4
)

δū : 0 = φ h x -[ h φ ] x -1 3 [ h 3 ūx ] x + h ū, (2.5
)

δh : 0 = 1 2 ū2 -g h + 1 2 h 2 ū 2 x -φt + φ ūx -[ ū φ ] x , (2.6)
thence

φx = ū -1 3 h -1 [ h 3 ūx ] x , (2.7) φt = 1 2 h 2 ū 2 x -1 2 ū2 -g h + 1 3 ū h -1 [ h 3 ūx ] x . (2.8)
Differentiating (2.8) with respect to x and using (2.7) in order to eliminate φ, after some elementary algebra, one obtains the surface tangential momentum equation (cf. Appendix B for the physical interpretation of this equation)

∂ t ū -1 3 h -1 (h 3 ūx ) x + ∂ x 1 2 ū2 + g h -1 2 h 2 ū 2 x -1 3 ū h -1 (h 3 ūx ) x = 0, (2.9)
that can be written into the more familiar non-conservative form

ūt + ū ūx + g h x = -1 3 h -1 ∂ x h 2 γ . (2.10)
After multiplication by h and exploiting (2.4), we derive the conservative equations for the momentum

∂ t [ h ū ] + ∂ x h ū2 + 1 2 g h 2 + 1 3 h 2 γ = 0, (2.11)
that can be rewritten

∂ t h ū -1 3 (h 3 ūx ) x + ∂ x h ū2 + 1 2 g h 2 -2 3 h 3 ū 2 x -1 3 h 3 ūū xx -h 2 h x ūū x = 0.
(2.12)

From the two conservative equations (2.4) and (2.11), an equation for the energy conservation can also be derived in the form

∂ t 1 2 h ū2 + 1 6 h 3 ū 2 x + 1 2 g h 2 + ∂ x ( 1 2 ū2 + 1 6 h 2 ū 2 x + g h + 1 3 h γ ) h ū = 0. (2.13)
In summary, we have just derived the following system of equations

h t + ∂ x [ h ū ] = 0, (2.14) ∂ t [ h ū ] + ∂ x h ū2 + 1 2 g h 2 + 1 3 h 2 γ = 0, (2.15)
together with (2.2), that are the Serre equations [35] and are a special case of the Green-Naghdi equations [START_REF] Green | On the theory of water waves[END_REF]. Here, we refer to these equations as the classical Serre-Green-Naghdi (cSGN) equations.

The derivation presented in this section is straightforward but, for readers more familiar with small parameter expansions, some additional material is given in Appendix A. Other variational derivations can be found, for instance, in [START_REF] Fedotova | Variational principle for approximate models of wave hydrodynamics[END_REF][START_REF] Kim | A derivation of the Green-Naghdi equations for irrotational flows[END_REF].

Physically, equations (2.14) and (2.15) describe, respectively, the mass and momentum flux conservations. The physical interpretation of (2.9) have been debated in the literature. Actually, (2.9) is a conservative equation for the surface tangential momentum, as shown in Appendix B (see also [17] for an alternative derivation).

Modified Serre's equations

The cSGN equations describe long waves in shallow water, thus the horizontal and temporal derivatives are small quantities, i.e., ∂ x ∝ O(σ) and ∂ t ∝ O(σ), where σ ≪ 1 is of the order of the water depth divided by the characteristic wavelength (Appendix A). As a consequence, the left-hand side of (2.10) is of order one, while its right-hand side is of order three (γ is of order two), while higher-order terms are neglected in this equation.

The definition (2.2) of γ involving parts of the equation (2.10), the latter can be used to derive another approximation of the vertical acceleration at the free surface. Indeed, substituting the relation

ūt + ū ūx = -g h x -1 3 h -1 ∂ x h 2 γ
, into the definition of γ, one obtains to the same order of approximation (Appendix A)

γ = 2 h ū 2 x + g h h xx + O σ 4 , (2.16)
which is another expression for the vertical acceleration consistent with the order of approximation. It is however possible to obtain a more general system averaging the two expressions (2.2) and (2.16) as

γ = 2 h ū 2 x + β g h h xx + (β -1) h ∂ x [ ūt + ū ūx ] + O σ 4 , (2.17)
where β is a parameter at our disposal. Thus, the modified Serre equations are

h t + ∂ x [ h ū ] = 0, (2.18)
∂ t [ h ū ] + ∂ x h ū2 + 1 2 g h 2 + 1 3 h 2 γ = 0, (2.19
)

2 h ū 2 x + β g h h xx + (β -1) h ∂ x [ ūt + ū ūx ] = γ, (2.20)
and the cSGN equations are recovered if β = 0. These modified SGN equations (mSGN), or similar ones, have been derived and studied before in the literature (see, e.g., [START_REF] Antunes | Boussinesq and Serre type models with improved linear dispersion characteristics: Applications[END_REF][START_REF] Dias | On the fully-nonlinear shallow-water generalized Serre equations[END_REF][START_REF] Liu | Two sets of higher-order Boussinesq-type equations for water waves[END_REF]).The free parameter β involved in the mSGN equation can be chosen such that it optimises the linear dispersion relation (see Section 4), but other criteria can be used to define β [START_REF] Dutykh | Adaptive modeling of shallow fully nonlinear gravity waves[END_REF]. From the modified Serre equations, one can derive a secondary relation

∂ t h ū + 1 3 (β -1) (h 3 ūx ) x + ∂ x h ū2 + 1 2 g h 2 + 1 3 β g h 3 h xx + 2 3 (2β -1) h 3 ū 2 x + 1 3 (β -1) h 3 ū ūxx + (β -1) h 2 h x ū ūx = 0, (2.21)
that is the modified counterpart of (2.12). However, equations analogous to (2.9) and (2.13) cannot be obtained if β = 0. This means, in particular, that the energy is not conserved for the mSGN equations and that a variational principle cannot be obtained (if β = 0). This is a serious drawback, specially for long time simulations and for theoretical investigations.

It should be noted that instead of using the generalised acceleration (2.17) into the momentum equation (2.15), it could be used into the energy equation (2.13). This yields other modified Serre's equations conserving the energy but not the momentum. Thus, the shortcomings of the mSGN are not addressed that way (if β = 0).

The mSGN equations (2.19) -(2.21) are asymptotically consistent with their cSGN counterparts for all values of the parameters β. However, the mSGN equations for the tangential momentum and for the energy are worse than asymptotically inconsistent: they are not conservative if β = 0. Even worse, a variational principle for the mSGN equations cannot be obtained if β = 0 . Thus, seeking for an improved shallow water model tuning the momentum equation, some physically crucial properties are lost, although the modifications are consistent from an asymptotic viewpoint.

In the next section, we derive another set of modified Serre's equations that have a free parameter and that conserve both the energy and the momenta.

Consistent modified Serre's equations

In order to derive a better variant of the SGN equations, we modify the Lagrangian (2.3) appropriately instead of modifying directly the cSGN equations, this modified Lagrangian being asymptotically consistent with (2.3). Thus, the existence of a variational principle for the resulting approximate equations is automatically ensured, as well as conservative equations for the momentum, energy and tangential momentum at the surface.

Modified Lagrangian

The Lagrangian density (2.3) for the cSGN equations can be rewritten replacing hū 2

x by its expression in terms of γ according to (2.2) such that:

L ρ = h ū2 2 + h 2 γ 12 + h 3 12 [ ūt + ū ūx ] x - g h 2 2 + { h t + [ h ū ] x } φ.
Substituting the generalised expression (2.20) for γ, we obtain the new Lagrangian L ′ 0 as

L ′ 0 ρ = h ū2 2 + h 3 ū 2 x 6 + β h 3 12 [ ūt + ū ūx + g h x ] x - g h 2 2 + { h t + [ h ū ] x } φ = L ρ + β h 3 12 [ ūt + ū ūx + g h x ] x . (3.1)
Obviously, L and L ′ 0 differ by fourth-order terms -the squared bracket in (3.1) -that are neglected in these approximations. These two Lagrangians have therefore the same order of approximation according to the asymptotic analysis (Appendix A). It should be noted that the squared bracket in (3.1) has a simple generalisation in three dimensions, thus allowing a straightforward three-dimensional extension of this modified Lagrangian (Appendix C).

Integrating by parts the terms involving ūt , ūū x and h xx , and using the mass conversation in order to remove the resulting term h t , one gets

h 3 [ ūt + ū ūx ] x = ∂ t [ h 3 ūx ] -3 h 2 h t ūx + h 3 ū 2 x + h 3 ū ūxx = ∂ t [ h 3 ūx ] + ∂ x [ h 3 ū ūx ] + 3 h 3 ū 2 x , h 3 h xx = ∂ x [ h 3 h x ] -3 h 2 h 2 x .
Substituting these relations into (3.1) and omitting the boundary terms

∂ x [• • • ] and ∂ t [• • • ]
(since they do not contribute to the variational principle), we obtain the simplified modified Lagrangian

L ′ ρ = h ū2 2 + (2 + 3β) h 3 ū 2 x 12 - g h 2 2 - β g h 2 h 2 x 4 + { h t + [ h ū ] x } φ. (3.2)
With this Lagrangian, it follows that the modified kinetic and potential energy densities are, respectively,

K ′ ρ = h ū2 2 + (2 + 3β) h 3 ū 2 x 12 , V ′ ρ = g h 2 2 1 + β h 2 x 2 ,
to be compared with the corresponding quantities in the cSGN (obtained when β = 0).

In the cSGN, the kinetic energy is approximate (compared to the exact equations), while the potential energy is exact. On the other hand, both energies are approximate in (3.2) if β = 0 . Having one (or more) exact relation in an approximated model is not necessarily a good thing; what really matters is to have a better approximation of the overall system of equations, as advocated in [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF]. It is more important to have a better approximation of the complete system of equations than, e.g., an exact potential energy.

Improved Serre's equations

The Euler-Lagrange equations for the modified Lagrangian (3.2) are

δ φ : 0 = h t + [ h ū ] x , (3.3
)

δū : 0 = h ū + φ h x -h φ x -1 3 + 1 2 β h 3 ūx x , (3.4
)

δh : 0 = 1 2 ū2 -g h -φt + φ ūx -[ ū φ ] x + 1 2 + 3 4 β h 2 ū 2 x -1 2 β g h h 2 x + 1 2 β g [ h 2 h x ] x , (3.5) thence φx = ū -1 3 + 1 2 β h -1 h 3 ūx x , (3.6) φt = 1 2 ū2 -ū φx -g h + 1 2 + 3 4 β h 2 ū 2 x + 1 2 β g h [ h h x ] x . (3.7)
After differentiation of (3.7) with respect to x and using (3.6), we obtain the improved Serre-Green-Naghdi (iSGN) equations, written as

h t + ∂ x [ h ū ] = 0, (3.8)
q t + ∂ x ū q -1 2 ū2 + g h -1 2 + 3 4 β h 2 ū 2 x -1 2 β g (h 2 h xx + hh 2 x ) = 0, (3.9)
where q = φx is given by (3.6). Equation (3.9) is the iSGN counterpart of the equation (2.9) for the cSGN, i.e. (3.9) is a conservative equation for the surface tangential momentum.

Analogs to the cSGN equations can also be derived.

A non-conservative momentum equation is obtained directly from (3.9)

ūt + ū ūx + g h x + 1 3 h -1 ∂ x h 2 Γ = 0, (3.10)
where

Γ = 1 + 3 2 β h ū 2 x -ūxt -ū ūxx -3 2 β g h h xx + 1 2 h 2 x .
Notice that Γ = γ if β = 0 . After multiplication by h, the equation (3.10) yields at once the conservation of the momentum flux

∂ t [ h ū ] + ∂ x h ū2 + 1 2 g h 2 + 1 3 h 2 Γ = 0 . (3.11)
Finally, the requested equation for the energy conservation can also be derived in the form

∂ t 1 2 h ū2 + ( 1 6 + 1 4 β) h 3 ū 2 x + 1 2 g h 2 1 + 1 2 βh 2 x + ∂ x 1 2 ū2 + ( 1 6 + 1 4 β) h 2 ū 2 x + g h 1 + 1 4 βh 2 x + 1 3 h Γ h ū + 1 2 β g h 3 h x ūx = 0. (3.12)
Thanks to the variational principle, we have derived a one-parameter generalisation of the Serre equations that retains important physical properties such as the momentum and energy conservations, as well as the Galilean invariance. These modified equations have a variational principle that is consistent with the original one according to the asymptotic analysis.

Steady waves

Since the iSGN equations are Galilean invariant, we can consider now steady solutions (i.e., independent of time in the frame of reference moving with the wave). For 2L-periodic solutions, the mean water depth d and the mean depth-averaged velocity -c are

d = h = 1 2L ˆL -L h dx, -c d = h ū = 1 2L ˆL -L h ū dx, (3.13)
thus c is the wave phase velocity observed in the frame of reference without mean flow. With h = h(x) and ū = ū(x) , the mass conservation (3.8) straightforwardly yields ū = -c d / h , and substitution into (3.9) and (3.11) gives, respectively,

(2 + 3β) F (2hh xx -h 2 x ) 12 (h/d) 2 - β (hh xx + h 2 x ) 2 (d/h) + F d 2 2 h 2 + h d = 1 + F 2 + C 1 , (3.14) (2 + 3β) F (hh xx -h 2 x ) 3 (h/d) - β (2hh xx + h 2 x ) 2 (d/h) 2 + 2 F d h + h 2 d 2 = 1 + 2 F + C 2 , (3.15)
F = c 2 /gd being a Froude number squared and where C n (n = 1, 2) are dimensionless integration constants to be determined from the relations (3.13) (C 1 = C 2 = 0 for solitary waves). Eliminating h xx between (3.14) and (3.15), after some straightforward algebra, one obtains the first-order ordinary differential equation

d h dx 2 = F -(1 + C 2 + 2F) (h/d) + (2 + 2C 1 + F) (h/d) 2 -(h/d) 3 1 3 + 1 2 β F -1 2 β (h/d) 3 . (3.16)
The left-hand side being positive, so is the right-hand side. Its numerator and the denominator are therefore necessarily of the same sign.

The equation (3.16) can be solved in parametric form with the help of Jacobian elliptic functions. However, we consider here only solitary waves, i.e., h(∞) = d and ū(∞) = -c thence C 1 = C 2 = 0 . Thus, with h(x) = d + η(x) , we obtain the differential equation

d η dx 2 = (F -1) (η/d) 2 -(η/d) 3 1 3 + 1 2 β F -1 2 β (1 + η/d) 3 . Introducing the change of independent variable x → ξ such that x(ξ) = ˆξ 0 (β + 2/3)F -β h 3 (ξ ′ )/d 3 (β + 2/3) F -β 1/2
dξ ′ , (3.17)

we obtain the classical equation for solitary waves

d η dξ 2 = (F -1) (η/d) 2 -(η/d) 3 1 3 + 1 2 β F -1 2 β
, with a solution given in the parametric form

η(ξ) d = (F -1) sech 2 κ ξ 2 , (κd) 2 = 6 (F -1) (2 + 3β) F -3 β , (3.18)
the wave amplitude being a = η(0) = (F -1)d . Substitution of (3.18) into (3.17) yields an expression for x(ξ). This expression being very complicated, though an explicit closed-form, it is not given here since it is of little practical interest (numerical quadrature is more efficient).

It should be noted that the cSGN and iSGN equations have the same relation between the phase velocity and the wave amplitude for steady solitary waves (i.e., the relation c(a) is independent of β), but other relations are generally not independent of β .

Criteria for choosing the free parameter

We present here some criteria for choosing the free parameter β The possibilities listed below are by no mean the only ones that could be considered.

Linear approximation

For infinitesimal waves, η and ū being both small, the iSGN equations yield the linear system of equations

η t + d ūx = 0, (4.1) ūt -1 3 + 1 2 β d 2 ūxxt + g η x -1 2 β g d 2 η xxx = 0. (4.2)
Seeking for traveling waves of the form η = a cos k(x -ct) , one obtains the (linear) dispersion relation

c 2 g d = 2 + β (kd) 2 2 + ( 2 3 + β) (kd) 2 ≈ 1 - (kd) 2 3 + 1 3 + β 2 (kd) 4 3 ,
that should be compared with the exact relation

c 2 g d = tanh(kd) kd = 1 - (kd) 2 3 + 2 (kd) 4 15 - 17 (kd) 6 315 + • • • . (4.3)
One can choose the parameter β so that the linear dispersion relation is correct up to the highest possible order of its Taylor expansion around k = 0 . Thus, with β = 2/15 ≈ 0.1333 the exact and iSGN dispersion relations coincide up to k 4 (but only up to k 2 if β = 2/15). It should be noted that this choice improves the dispersion relation only in the vicinity of kd = 0 , i.e., for very long infinitesimal waves. Another possibility is to minimise some norm of the difference between the exact and approximate linear dispersion relations over an interval k ∈ [0, k max ] with k max chosen a priori.

Decay of a solitary wave

In the far field, a solitary wave decays like e -κ|x| , the trend parameter κ > 0 being related to the Froude number by the "dispersion" relation *

c 2 g d = 2 -β (κd) 2 2 -( 2 3 + β) (κd) 2 = 1 + (κd) 2 3 + 1 3 + β 2 (κd) 4 3 + • • • , ( 4
.4) * Recall that we are considering the frame of reference moving with the wave. In a fixed frame, the wave travels with a constant phase speed c and solitary waves decays, both in space and time, like e -κ|x-ct| . Thus, relations (4.4) and (4.5) are not dispersion relations for the evanescent modes, the latter decaying in space and being periodic in time.

to be compared with the exact relation [START_REF] Mccowan | On the solitary wave[END_REF][START_REF] Stokes | The outskirts of the solitary wave[END_REF]

c 2 g d = tan(κd) κd = 1 + (κd) 2 3 + 2 (κd) 4 15 + 17 (κd) 6 315 + • • • . (4.5)
The main advantage of this approach compared to the previous one is that (4.5) is exact and valid for all amplitudes, while (4.3) is exact only in the limit a → 0 . Note that (4.5) can be simply deduced from (4.3) by setting k = iκ . The tangent function is meromorphic with single poles at κd = ±π/2, ±3π/2, • • • . The relation (4.4) has a single pole at κd = ± 6/(2 + 3β) . This pole is at κd = ±π/2 if

β = 2 3 12 π -2 -1 ≈ 0.
1439 . This choice of β is close to the optimal value 2/15 proposed above. The fact that both methods (Taylor expansion and meromorphic interpolation) produce similar optimal values for β is an indication that β = 2/15 is a good choice.

Limiting wave

For exact steady solutions, we focus now on the vicinity of the crest of the highest waves. To this end, we consider the solitary wave solution given above without loss of generality since the analysis is local. Taylor expansions of (3.17) and (3.18) around the crest ξ = 0 give

x(ξ) = (2 + 3β) F -3 β F 3 (2 + 3β) F -3 β 1 2 ξ + O(ξ 3 ), η(ξ) = a 1 -1 4 (κξ) 2 + O(ξ 4 ).
Peaked crests are found if dx/dξ = 0 at ξ = 0 . That happens if

F 2 = 1 + 2 3 β -1 .
With this peculiar limiting Froude number, the slope at the crest is

d η dx x=0 ± = d η dξ d x dξ ξ=0 ± = ∓ κ a √ F 2 + F + 1 √ 3 F = ∓ a 2 / 3β a + d .
In other words, if β = 0 , we have a family of limiting waves parametrised by the parameter β . Various inner angles are obtained depending on the choice of β . For β = 2/15 , we have

F = √ 6 ≈ 2.45 , a/d = √ 6 -1 ≈ 1.
45 and the crest has (approximately) a 28.4 • inner angle. This is therefore not a good approximation of the limiting solitary wave.

For β = 2 3 (12π -2 -1) , the crest inner angle is about 37.3 • . This is thus a somehow better approximation but still far from the exact 120 • angle [START_REF] Stokes | Supplement to a paper on the theory of oscillatory waves[END_REF].

From these peculiar choices, it can be seen that the iSGN system optimised for linear waves cannot describe equally well the highest wave. A 120 • inner angle is obtained if and in terms of β , this condition becomes

κ a √ F 2 + F + 1 = F or a 2 / β = a + d,
3 β 4 = β + √ 2 β -1 β - √ 8 β + 2 .
This relation gives β ≈ 0.34560 thence a/d ≈ 0.711 , that is much closer to the exact value a/d ≈ 0.8332 [31]. However, with this value of β , the linear dispersion relation of infinitesimal waves is not at all improved.

Remark 1. The iSGN model provides limiting values for wave heights, unlike the cSGN one, that is an interesting feature in some situations. However, this model is not designed for modelling extreme waves, choosing β to improve the linear dispersion relation is probably the best choice for most simulations.

Numerical examples

In this section, we present few numerical tests illustrating the performance of the iSGN model. The numerical methods used are briefly described in corresponding sections below. The aim here is just to show that the new model outperforms the classical one.

Steady solitary waves

For the first test, we compare steady solitary wave solutions of the cSGN and the iSGN (with β = 2/15) models with the exact solution of the full Euler equations [9,14]. The comparison of surface elevations for three different amplitudes is given in Figure 2. It can be seen that both models are good approximations for small amplitudes. For large amplitudes, however, the iSGN model is significantly better than the cSGN one. This result may be surprising because the value β = 2/15 is chosen to improve the dispersion of sinusoidal infinitesimal waves (see Section 4.1). This is not so surprising considering that this value of β is close to the one optimising the the relation (4.4) relating the phase speed with the trend parameter (see Section 4.2).

Random wave evolution

We consider now a random initial condition (see upper Figure 3(a)) which is generated from the Gaussian Fourier power spectrum with random phases distributed uniformly in ]0, 2π] (see [START_REF] Dutykh | Adaptive modeling of shallow fully nonlinear gravity waves[END_REF] for details where the same type of initial conditions were used). This initial wave is moderately steep (a/d = 0.1) but rather dispersive (σ = 0.25). Thus, this test focuses on the dispersive properties of the iSGN model. In order to solve the iSGN equations on a periodic domain, we use the continuous Galerkin/Finite-Element method detailed in [START_REF] Mitsotakis | On the Galerkin/Finite-Element Method for the Serre Equations[END_REF]. The unsteady Euler equations are solved numerically using conformal mapping and pseudo-spectral discretisation [7,[START_REF] Li | A Numerical Study of the Exact Evolution Equations for Surface Waves in Water of Finite Depth[END_REF].

The evolution of the generated random initial condition is shown in Figure 3. In Figure 3(a) one can see that all three models are initialised with the same condition. During the evolution of this initial condition, the cSGN model starts lagging behind the reference solution given by the full Euler model. In Figure 4, we show a magnification of a portion of the computational domain at the final simulation time, where one can see that the iSGN equations commit a much lower phase error. It results a better prediction in the positions of the wave crests, that is important in applications focused on extreme waves.

Discussion

In this paper, we presented a novel model for fully nonlinear long waves in shallow waters. These equations generalise the Serre equations in the sense that the new model contains a free parameter. The value of this parameter can be used to optimise, for instance, the linear dispersion relation properties. However, it is possible to use this extra degree of freedom to optimise the linear shoaling [3] or any other desired characteristics.

In all the modified shallow water models we are aware of, the introduction of an extra parameter leads to the violation of the energy conservation. For some models, the Galilean invariance is violated too. Here, we succeeded in deriving a new model that preserves these properties, thanks to the use of a variational formalism based on Hamilton principle [START_REF] Clamond | Practical use of variational principles for modeling water waves[END_REF]. Instead of tweaking the system of PDEs, we modified the corresponding Lagrangian functional. This approach allows to preserve the underlying variational structure as well. Moreover, we showed that this model repeats another feature of the full Euler equations: the existence of limiting waves that are not included in the Serre equations.

We performed some unsteady simulations that clearly show the importance of the iSGN equations compared to the classical ones. The fact that the energy is conserved in the iSGN model is very important for long-time simulations. Also, we have shown that the improved dispersion relation leads to a better description of large solitary waves.

The main goal of this study was to make a proof of principle. However, the developments and ideas presented above can be generalised to non-flat bottom and in three dimensions (see Appendix C), for example. The derivations and examples presented here are sufficient to demonstrate the advantages of our approach.

Tweaking the vertical acceleration at the surface, we easily obtained a one-parameter generalisation of the classical Serre equation. Modifying other quantities, in the spirit of [START_REF] Dias | On the fully-nonlinear shallow-water generalized Serre equations[END_REF], one may obtain a multi-parameter generalisation of these equations. These extra free parameters can be tuned to improve even further the model equations. These generalisations are also left to future investigations.

We conclude this discussion by noting that, in addition to the Galilean invariance and the energy conservation, the dispersion-improved iSGN equations have the same order of derivatives and are thus not more difficult to solve numerically (similar algorithms and running times) compared to the classical Serre equations, unlike the high-order Boussinesq-like equations. This is an interesting feature for practical applications. energies are

K ρ = ˆη -d u 2 + v 2 2 dy = h ū2 2 + h 3 ū 2 x 6 + O σ 4 , V ρ = ˆη -d g (y + d) dy = g h 2 2 .
The incompressibility and the bottom impermeability being identically satisfied, the Hamilton principle can be reduced to the Lagrangian density

L ρ = h ū2 2 + h 3 ū 2 x 6 - g h 2 2 + { h t + [ hū ] x } φ + O σ 4 ,
which, after neglecting higher-order terms, corresponds exactly to the Lagrangian density in the action integral (2.3).

It should be noted that instead of enforcing the mass conservation, one could enforce the impermeability of the free surface. These two approaches are equivalent here, however. Indeed, the incompressibility and the bottom impermeability are identically fulfilled with (2.1), fulfilling the surface impermeability yields the mass conservation, and vice versa.

B. Tangential momentum at the free surface

Let φ be the velocity potential of an irrotational motion, i.e., u = φ x , v = φ y . Denoting with 'tildes' the quantities written at the free surface * and exploiting the surface impermeability, we have the relation φ t = φt -η t ṽ = φt -ṽ2 + η x ũ ṽ = φt -ṽ2 + φx -ũ ũ.

(B.1)

The Bernoulli equation at the free surface can then be rewritten φt + φx -ũ ũ + g η +

C. Generalisation in three dimensions

The generalisation of our approach in three dimensions is straightforward. Let x = (x 1 , x 2 ) the horizontal Cartesian coordinates and u = (u 1 , u 2 ) the horizontal velocity where ∇ is the horizontal gradient. From this ansatz, one can derive the 'irrotational' Green-Naghdi equations [9,[START_REF] Kim | A derivation of the Green-Naghdi equations for irrotational flows[END_REF]

h t + ∇ • [ h ū ] = 0, ūt + ∇ | ū| 2 2 + g ∇h + ∇ • [ h 2 γ ] 3 h = ū • ∇h 3 ∇[ h ∇ • ū ] -ū • ∇ h∇ • ū 3
∇h,

where γ = h{(∇ • ū) 2 -∇ • ūt -ū • ∇[∇ • ū]
} is the vertical acceleration at the free surface. These equations can be straightforwardly derived from the Lagrangian density

L ρ = h | ū| 2 2 + h 3 (∇ • ū) 2 6 - g h 2 2 + { h t + ∇ • [ h ū ] } φ.
Along the lines of the two-dimensional case (see equation (3.1)), an obvious alternative Lagrangian is

L ′ = L + 1 12 ρ β h 3 ∇ • ūt + 1 2
∇| ū| 2 + g ∇h , where β is a dimensionless constant at our disposal. According to the asymptotic analysis of Appendix A, L is obtained neglecting terms of order σ 4 (and higher) and we have ∇ • ūt + 1 2 ∇| ū| 2 + g ∇h = 0 + O(σ 4 ). Thus, clearly, L and L ′ differ by fourth-order terms, so they are both consistent to the same order of approximation considered here. From L ′ , one can easily derive modified equations of motion with a free parameter that can be chosen to improve the linear dispersion relation, for example. These equations yield automatically the energy conservation. The derivations are left to the reader.

There are several 3D generalisations of the Serre equations that have been debated in the literature [START_REF] Dell'isola | Variational Models and Methods in Solid and Fluid Mechanics[END_REF]. We picked one of them in order to illustrate that the improvement method can be easily used in three dimensions; our goal here is not to explore in details all the possible 3D cases. It is also clear that the approach presented here can be easily used to improve models in presence of a varying bottom and surface tension, for example.

Contents 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Classical and modified Serre's equations . . . . . . . . . . . . . . . . . . . . . 2.1 Classical Serre's equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Modified Serre's equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Consistent modified Serre's equations . . . . . . . . . . . . . . . . . . . . . . . 3.1 Modified Lagrangian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Improved Serre's equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Steady waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Criteria for choosing the free parameter . . . . . . . . . . . . . . . . . . . . . 4.1 Linear approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Decay of a solitary wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Limiting wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Steady solitary waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Random wave evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A Asymptotic derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B Tangential momentum at the free surface . . . . . . . . . . . . . . . . . . . C Generalisation in three dimensions . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 .

 1 Figure 1. Sketch of the domain.

Figure 2 .

 2 Figure 2. Comparison of steady solitary waves. '-': Euler; '--': iSGN (β = 2/15); '-• -': cSGN (β = 0).

Figure 3 .

 3 Figure 3. Evolution of a random initial condition in three different models. '-': Euler; '--': iSGN (β = 2/15); '-• -': cSGN (β = 0).

Figure 4 .

 4 Figure 4. Magnification of Figure 3(c). '-': Euler; '--': iSGN (β = 2/15); '-• -': cSGN (β = 0).

  field. A shallow water ansatz fulfilling the fluid incompressibility and the (horizontal) bottom impermeability isu(x, y, t) ≈ ū(x, t), v(x, y, t) ≈ -(y + d) ∇ • ū, (C.1)

  Equation (2.9) is recovered using the ansatz (2.1) and the relation (2.7). Therefore, Equation (2.9) represents the horizontal derivative of the Bernoulli equation at the free surface, i.e., the surface tangential component of the momentum (Euler) equation for irrotational flow. Alternative derivations of this result can be found in[17].

			1 2	ũ2 -1 2	ṽ2 = 0.	(B.2)
	After differentiation with respect of x, this equation becomes
	∂ t φx + ∂ x g η + 1 2	ũ2 -1 2	ṽ2 + φx -ũ ũ = 0.	(B.3)

Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgments D. Dutykh & D. Clamond acknowledge the support from CNRS under the PEPS InPhyNiTi project FARA. D. Mitsotakis was supported by the Marsden Fund administered by the Royal Society of New Zealand.

A. Asymptotic derivation

For the horizontal velocity u, the most general solution of the Laplace equation satisfying the seabed impermeability ([23]) is

where ȗ(x, t) = u(x, y = -d, t) is the velocity at the bottom. Assuming long waves in shallow water means that ∂ x = O(σ) and, these waves having finite phase velocities, ∂ t = O(σ). Thence, the relation defining the depth-averaged horizontal velocity ū, i.e.

) , can be solved for ȗ as

) . Thus, the horizontal velocity (A.1) can be represented equivalently to O(σ 4 ) as

Similarly, fulfilling the fluid incompressibility and the bottom impermeability, the vertical

, thence at the free surface ṽ = -hū x -1 3 h 2 h x ūxx + O(σ 5 ) . Assuming finite velocity and surface elevation, we take ū = O(σ 0 ) -thus φ = O(σ -1 ) because u ≈ φx ) and η = O(σ 0 ) -and the depth-integrated kinetic and potential


INTRODUCTION

Power grid reliability requirements enacted by international and national level standards require avoiding electrical service outages. The outage of power delivery must be in the range of no more than 5 minutes per year to achieve higher availability (99.999%).

In [START_REF] Altaher | Evaluation of Time-Critical communications for IEC 61850-Substation Network Architecture[END_REF], experimental setup was conducted to evaluate critical communication in a testbed of process-level network. The evaluation is done, with framework of the IEC 61850 standard, by implementing communication protocols. In this paper we aim to investigate dependability of bay-level hardware architectures inside an IEC 61850-based distribution substation. Our approach uses reliability block diagrams and the safety integrity approach framework from the international standards IEC 61508, to investigate reliability, availability and safety availability in this context.

This paper is organized in five sections. Second section identifies briefly the electrical power grid. Third section introduces the IEC 61850 substation communication standard, its features and concepts; in sub-section 3 we emphasis a distribution substation. In section 4, we detail the dependability concept and its attributes. Proposed architecture with reliability and availability calculations are mentioned in sub-section 3 whereas sub-section 4 proposes an evaluation of the functional safety. Section 5 concludes the paper.

THE ELECTRICAL POWER GRID

The electrical power delivery follows many stages from power generation through transmission and distribution to the ultimate main load centers. Practically the electrical power grid composes of power stations, cabling system, transmission and distribution substations and related control centers. The interconnection of these subsystems becomes progressively intelligent forming the smart grid. The grid exchanges automatically real time, electrical power, parameters enabling utilities to manage and control remotely all components and to manipulate the grid reliability and safety [START_REF] Mcdonald | Distribution Systems, Substations, and Integration of Distributed Generation[END_REF]). Modern grids have telecommunication networks supporting distributed protection and control functions inside the substations. Further, the power system includes extra-high voltage (EHV) networks used to transmit electricity from power generation stations to distribution networks known as high-voltage (HV) networks. End consumers are connected to the distribution facilities, i.e. substations. Electrical substations build coordinated transmission and distribution nodes, which are configured by control centers Dependability Optimization of Process-level Protection in an IEC-61850-Based Substation A. Altaher, S. Mocanu & J.-M Thiriet Univ. Grenoble Alpes, GIPSA-lab, F-38000 Grenoble, France ABSTRACT: Power substations are intensively renovated toward using information and communication technologies such as object oriented modeling and Ethernet networks. In the last two decades, Substation automation systems used capabilities of network communication services adopted from sophisticated international standardizations such as IEC 61850. Distributed safety related functions take advantage of these technologies to protect the process-level equipment. Substation devices such as intelligent electronic devices, measurement units and circuit breaker controllers, with new capabilities, i.e. enabling IEC 61850, are integrated to build the protection and control functions that form the safety-related system. The objective of this research is to evaluate quantitatively the dependability for transformer protection architectures in the bay level. Safety integrity levels model, described in both IEC 62061 and IEC 61508, gives measurements for safety integrity levels according to the probability of failure. The determination of these levels is an approach to estimate system dependability.

to maintain the health of the power grid. Smart technologies and intelligent devices protect and control substation functionalities either locally or remotely by means of information and communication technologies. The IEC 61850 standard, communication networks and systems in substations, drives interoperability between several substation devices, e.g. intelligent electronic devices (IEDs), from different vendors. Additionally, the standard provides flexible assignment of different functions that enhance safety-related functions [START_REF] Brand | Design of Based Substation Automation Systems According to Customer Requirements[END_REF].

THE IEC 61850 STANDARD

The technical committee 57 (TC57), which belongs to The International Electro-technical Commission (IEC), has released the IEC 61850 standard to enforce interoperability and to enable abstraction of communication services (IEC TC57. 2010). Its parts manage standardization of different services dedicated for vertical and horizontal communications inside substation levels. It enables delivery of highspeed peer-to-peer messages for status and events exchange. The standards includes at least 10 parts, the first five parts identify general and specific functional requirements. Part 6 covers an engineering tools such as substation configuration language to manage the design configurations by allowing description of relations between substation functions, substation primary or/and secondary equipment. Part 7 includes 4 sub parts identifying abstraction and mapping of data services to communication protocols. Part 8 defines manufacturing message services (MMS) mapping to data services. Part 9 also has sub parts, part 9-1 defines mapping of sampled values (SV) and 9-2 deals with the Process-bus. Part 10 emphasizes procedures of conformance testing.

IEC 61850 Standard features

The substation automation systems benefit from the standard parts at different levels (station, bay and process levels). Inherited engineering models are enforced by the standard providing new capabilities such as object-oriented data modeling, communication mapping services, configuration tools such as substation configuration language, and new processlevel instrumentation technologies. These capabilities paved the way for horizontal communications between distributed devices inside substation levels. In result, the IEDs cooperate in real-time, to enable distribution of functionalities, by using Ethernet network technology. The substation primary equipment is controlled by dedicated formal device function numbers (IEEE std. C37.2 2008). The IEC 61850 part 5 introduces standardized logical nodes (LN) to define these functions. Logical nodes are embedded into the IEDs to form logical devices that are defined by specific requirements. Obviously, the IEDs are programmable devices that provide protection and control functions. These devices contain logic solvers, input and output ports, network interfaces to gather information about the primary equipment via communication protocols. GOOSE (general object oriented substation events) messages are multi-casted from publishers (IEDs) to subscribers (IEDs) in real-time pattern. These messages are important for delivering status (data set) about primary equipment parameters, e.g. circuit breaker switch position and status. Other specific messages called sampled values (SV) are used to inform baylevel devices about the process-level measured physical parameters, i.e. voltage and current. Innovative merging units (MU) are used to convert analog parameters and to transmit synchronously the sampled values with precise time-synchronization (IEC 61850 part 9-2).

The IEC 61850 object modeling

The standard defines object models to enforce interoperability between IEDs from different vendors. These models help reducing costs and time required for configuration of SAS implementations, and improving engineering documentation. As mentioned earlier logical nodes are models that form basic functionalities. A logical node is the smallest part of a function that exchange data (IEC TC 57. 2010). It contains data object and methods used to create functional components of different classes. These classes provide different protection, measurement, control, monitoring functions in SAS operations. Among these IEC 61850 defined classes are XCBR circuit breaker, XSWI isolator or earth switch, TCTR Current transformer, YLTC power transformer, PTOC time overcurrent protection and ATCC automatic tap changer controller. Interoperability between logical nodes achieved through standardized data objects included in every logical node. Three levels of data modeling and services are considered. The first level is the abstract communication service interface (ACSI) that allocates models and associated services for accessing data and object elements. The second level is common data classes (CDC) that specify data attributes while the third level defines logical node classes and data classes. These levels represent hierarchical levels (Fig. 1).

Alongside this concept, an intelligent electronic device, i.e. physical device PHD, could be constructed by one or several logical devices.

The distribution substation under study

We choose the transformer bay (Fig. 2), in distribution substation, as a case study to analyze the functional components according to the IEC 61850 standard. We start by identifying specific components such as primary and secondary equipment. The primary equipment is the main process-level circuits composed of a bus bar, power lines, feeders and transformers, while secondary devices are bay-level auxiliary devices such as monitoring, protection and control devices. The bay-level in the substation is used logically between station-level and processlevel equipment. The station-level consists of engineering computers, human machine interfaces, supervisory control and data acquisition, tele-control access, whereas the process-level is composed of electrical power equipment, e.g. circuit breakers.

Planning a SAS design requires drawing of a single line diagram (Fig. 2) that shows process-level components. Electromechanical equipment such as disconnectors and CBs are shown. The power transformer, i.e. converting high voltage into a lower voltage levels, converts 34.5 kV into 13.8 kV. This transformer bay forms a small distribution (D1-2) substation architecture (IEC 61850-1 2010). In other hand, additional components: bus bar, line, breakers and disconnectors are interconnected to construct the primary switchgear. These components are controlled by local/remote commands via a local Ethernet network. The bay-level would include protection and control IEDs that handle functionalities of the process-level and gather information about the equipment. The protection and control IEDs are interconnected via communication network (Fig. 2).

DEPENDABILITY OF SAS SYSTEMS

Dependability is defined as the trustworthiness of a computer system such that reliance can justifiably be placed on the service it delivers. The service delivered by a system is its behavior as it is perceived by its users; a user is another system (human or physical) which interacts with the former [START_REF] Laprie | Dependability: basic concepts and terminology in English, French, German, Italian and Japanese[END_REF]). Additionally, it is defined as the ability to perform as and when required . Dependability is a collective set of time-related performance characteristics that coexist with other requirements such as output, efficiency, quality, safety, security and integrity (Hardeveld & Kiang 2012).

Reliability and availability of process level

Reliability is defined as continuity of service, in other definition it is defined as a measure of continuous delivery of correct service or, equivalently, of time to failure [START_REF] Laprie | Dependability: basic concepts and terminology in English, French, German, Italian and Japanese[END_REF]) while IEC TC56 committee defines reliability as the probability that an item fulfils the required functions for the required duration.

In power utility communication systems, IEC 61850 part 3 section 4 considers reliability as quality requirement by focusing on communications for substation automation networking services. The standard provides a reference for other standards such as IEC 60870-4 that specifies performance requirements for tele-control. Further, IEC 61850 identifies the reliability of communications, inside the substation different levels, as data exchange without failure, loss or delay of critical messages. Precisely, there should be no single point of failure (SPOF) in substation networks. If failures exist, outcomes may cause damage to substation equipment. The standard insists that communication reliability is needed as a requirement for substation automation; therefore, it recommends fail-safe design that should be handled to avoid undesired control action.

The availability is defined as readiness of usage, and as a measure of the delivery of correct service with respect to the alteration of correct and incorrect service [START_REF] Laprie | Dependability: basic concepts and terminology in English, French, German, Italian and Japanese[END_REF]. IEC TC56 illustrates: "the availability is extent to which an item is operational and able to perform any required function or set of functions if a demand is placed on it". Evidently, one can recognize the relation between the dependability attributes. The availability is derived from reliability, safety and integrity (Hardeveld & Kiang 2012). Hence, the dependability has different attributes [START_REF] Laprie | Dependability: basic concepts and terminology in English, French, German, Italian and Japanese[END_REF]) such as reliability, availability, safety and security.

In addition, IEC TC56, which is responsible for dependability management standards, defines core dependability attributes that include reliability, availability, maintainability and maintenance support. The current practice in the field of power transmission and distribution requires dependable networks to deliver continuous electrical energy. Power substations contribute directly to the power system reliability and availability. Transmission and distribution substations process-level equipment is the main concern for evaluating substation role on general power system reliability. The substation availability depends mainly on primary and secondary equipment failure rates, communication topology, redundancy, maintenance procedures and support…

In previous work [START_REF] Altaher | Evaluation of Time-Critical communications for IEC 61850-Substation Network Architecture[END_REF], the authors investigate reliability of process-level network within an IEC 61850 enabled communication services. In this paper, we aim to evaluate hardware architectures considering redundant components.

Safety and functional safety

Safety is avoidance of catastrophic consequences on people, property and environment; hence, reducing risk frequency and consequence is primary objective for safety systems. In this manner, safety is a measure of continuous safeness, or equivalently, of time to catastrophic failure [START_REF] Laprie | Dependability: basic concepts and terminology in English, French, German, Italian and Japanese[END_REF]). Functional safety is a part of overall safety that depends on the correct functioning of the process or equipment in response to its inputs (IEC TC 65. 2010). Protection functions in substations were found to be safety related with varying levels of risk [START_REF] Purewal | Functional safety in application of programmable devices in power system protection and automation[END_REF]). These functions construct principal protection layer to prevent hazardous. Among these hazardous events are; short-circuits, arc flash and inter-phase short-circuits. Switchgear equipment faults could lead to critical failures such that failing to force sequential clearance of faults. In result, these events cause hazard consequences against substation technicians.

For practical design reasons, disconnectors in process-level are not able to switch on/off higher current loads [START_REF] Zurakowski | Safety and Security Issues in Electric Power Industry[END_REF]. Technically, these disconnectors are used for isolation clearance between disconnected elements, which cannot be achieved by circuit breakers. Switching-off a line requires sequential procedure performed in timely adjusted mechanism to avoid any damage. Perfect example is when tripping circuit breaker, in case of failing or delayed command that might cause late reaction. Consequence could be an electrical arc accompanied by high rate optic and acoustic phenomena, causing spray of melted iron and flashovers.

Another consequences may include inter-phase short-circuit leading to damage of equipment and injury of nearby people. The loss of disconnector would cost procurement of new equipment, replacement time, partial power outage and probability of critical accidents. E.g. the wrong switching-off of a loaded line by the related disconnector. The potential consequences will cost damage of process-level equipment and switchgear devices. In addition, substation disturbance could cause partial interruption of power distribution.

In [START_REF] Prata | Failure Risk Associated with Different Substation and HV network Configurations[END_REF]) failure mode analysis was performed on MV substation to identify the components that were more likely to cause an interruption. That analysis has shown the most critical component in the studied substation is power transformer. Another failure mode analysis was done by [START_REF] Zurakowski | Safety and Security Issues in Electric Power Industry[END_REF] to evaluate risks and their consequences on HV power substations. In our study, we consider both; a) the above mentioned failure modes that cause risks for substation staff and switchgear equipment, and b) the risk matrix which is elaborated by [START_REF] Prata | Failure Risk Associated with Different Substation and HV network Configurations[END_REF].

Integration and separation

Many standards recommend separation of control systems from safety functions, e.g. Institute of Electrical and Electronics Engineering standards (IEEE) recommend that the safety system design shall be such other systems failures shall not prevent the safety system form meeting its requirements. Furthermore, ANSI/ISA (American National Standards Institute/International Society of Automation) mentioned that separation between basic process control system (BPCS) and safety instrumentation system (SIS) functions reduce the probability that both control and safety functions become unavailable at the same time, or that inadvertent changes affect the safety functionality of the SIS. Contrary, in many substation design requirements, protection and control systems are integrated. Functional safety imposes separation of protection functions from control system to avoid common cause failures.

In the protection and control, the process of specifying the degree of targeted safety should concern the associated risk assessment. Concerning the safety probability per function, designers of SAS safety-related functions estimate a tolerable failure with probability of 10 -5 /h to 10 -6 /h [START_REF] Brand | Design of Based Substation Automation Systems According to Customer Requirements[END_REF].

One measure for the functional safety is the probability of failure on demand when a system fails to respond to a demand for action arising from a potentially hazardous condition (IEC TC 65. 2010). This parameter increases during mission time or test interval.

Reliability, availability and safety metrics

For dependability evaluation, the combinatory model of reliability block diagram (RBD) is used to illustrate the functional components, and for analyzing different system architectures such as parallel, series, non-parallel and series. A system works if there is path of functioning components. Comparing other methods such as fault trees (FT), reliability graphs, Bayesian networks and Markov models, this model is an effective tool that provide flexibility to determine reliability of system by considering its components. Evaluation of dependability concerns; hardware components, communication network as component and redundancy of critical components to avoid single point of failure. Protection and control subsystems are allocated in series arrangements, while redundant components are represented by parallel arrangements. These arrangements are required when one subsystem replace another subsystem. Assuming useful life period and constant failure rates (i.e. exponentially distributed lifetime), the following formulas give expressions for calculating reliability by given metrics such as mean time to failure (MTTF), mean time to repair (MTTR), i.e. average time from detecting the failure of a component until its replacement. Eq.1 gives an expression of failure rate as reciprocal of MTTF

i MTTF i 1 = λ (1)
Where, λi is failure rate of i th component. The reliability function of component i is given by Eq.2:

t i λ ) ( - = e t i R
(2) Where t is the mission time for a component in useful life period. Series system reliability, Rs(t) is given by Eq.3, assuming independent failure of individual components.

)t λ ( 1 i 1 ) ( ) ( ∑ = = = - = ∏ n i e t R t R n i i s
(3) Parallel system reliability (non-repairable components), Rp(t), correspondingly, is given by Eq.4:

∏ = - = n i i p t Q t R 1 ) ( 1 ) ( (4)
Where Qi(t) is the unreliability of i th component, and it is given by Eq.5:

t λ i 1 ) ( - - = e t Q i (5)
Furthermore, to evaluate availability, following formulas use repair (µ) and failure (λ) rates to calculate the availability. Repair rate (replacing a component) is given by Eq.6:

i i MTTR 1 = µ (6)
Then, availability Ai of i th component will be calculated by using MTTF and MTTR metrics as given in Eq.7:

i i i i MTTR MTTF MTTF A + = (7)
Similarly, the availability formula, using repair and failure rates, is given in Eq.8:

i i i i A λ µ µ + = (8)
To evaluate the safety availability, tolerable risk/failure, non-protected failure is required to de-termine the risk reduction factor RRF. Eq.9 gives an expression to find the RRF:

t np F F RRF = (9)
Where Fnp is the unprotected risk frequency that equals dangerous undetected failures, i.e. PFD for low demand and PFH for high demand, and Ft is tolerable risk frequency which is mentioned early (see section 4.3). PFDavg is given by Eq.10:

RRF avg PFD 1 = (10)
where RRF represents the risk reduction factor. In Eq.11, safety availability is given by knowing the PFDavg.

avg PFD SA - = 1 (11)
The following table (Table 1) contains the MTTF values that are obtained from [START_REF] Brand | Substation automation handbook[END_REF][START_REF] Lindquist | Circuit breaker failure data and reliability modelling[END_REF], these figures will be used for reliability, availability and safety availability calculations. From the given metrics, reliability, i.e. mission time is 1000 hours, and availability are calculated for the three architectures (table 3). Average Probability of Failure on Demand-Low rate demand >=10 -2 to 10 -1 >=10 -3 to 10 -2 >=10 -4 to 10 -3 >=10 -5 to 10 -4

Failure rate (λ) per hour -high rate demand >=10 -6 to 10 -5 >=10 -7 to 10 -6 >=10 -8 to 10 -7 >=10 -9 to 10 -8

The results are given in table 3 as following: The IEC 61508 standard sets SIL level according to dangerous failure rates as explained in table 3. We assume as in the standard (IEC TC 65. 2010) that statistically only every other failure is a potentially dangerous failure. For complex devices failure modes are assumed and divided as 50% safe and 50% dangerous. This assumption is explained in Eq.12.

λ λ λ × = = 2 1 s D (12)
The Eq.1 thus holds for λD, which equals in this case 1/MTTFD, for determining both dangerous undetected λDU failure rate and dangerous detected λDD, Eq.13 and Eq.14 are used: For redundant components with MooN structure, designers should consider common cause failures, represented by β factor (IEC TC 65. 2010). For 1oo2 structure, we assume that β=0.04, and βD=0.02,i.e. dangerous common cause factor, for both redundant Ethernet switches and bay IEDs. PFH is given by:

) 2 1 ( 2 ) ) 1 ( ) 1 (( 2 MTTR T DU MTTR DD D GE t CE t DU DD D PFH + × × + × × + × × × - + × - × = λ β λ β λ β λ β (17)
Where tGE (group down time) is given by Eq.18: and DC is the (automatic) diagnostic coverage (assumed 90% coverage) and T1 is the proof-test interval (assumed 1 month). Eq.19 gives the total failure rate of the SIF system (series):

∑ = = n i i ser 1 λ λ (19)
For two identical components with constant failures, the total failure rate, i.e. 1oo2 redundancy is given by Eq.20 The above-mentioned formulas are used for estimating the PFH of the three safety related architectures (Fig. 3) and the results are tabulated in (table 4). From table 3 andtable 4, we can notice clearly that two architectures with redundant components are satisfying IEC 61850 requirement, i.e. avoiding single point of failure, the architecture 3 has a higher reliability protection function regarding its lower dangerous failure and overall failure rate. These architectures are all suitable for SIL 1 level (high demand). Figure 4 illustrates decreasing of the reliability with mission time increasing. We evaluate three process-level architectures by classical reliability methods. Safety availability of these architectures is evaluated by using functional safety approach.

First architecture is basic well-known processlevel/bay-level topology, which already exists in many distribution substations. The second proposed architecture benefits from Ethernet switch redundancy to eliminate SPOF in case of communication faults, furthermore, the third proposed architecture uses redundant bay device to increase availability of main-bay protection and control. In future research, we will evaluate availability of shared logical nodes, therefore a detailed analysis of Ethernet communication network traffic, and the switching mechanism in relation with worst-case load will be performed

Figure 1 .

 1 Figure 1. Hierarchical object modeling

Figure. 2

 2 Figure.2 Substation communications among different levels, Logical Nodes within IEDs

Figure. 3 :

 3 Figure.3: three different process/bay level architectures

  channel equivalent mean down time (tCE) is calculated according to the standard in Eq.15

  Calculating functional safety metrics

Figure. 4 :

 4 Figure.4: reliability of the three architectures

Table 1 :

 1 Components MTTF and MTTR values

	Component	MTTF	MTTR
		(Years)	(Hours)
	Bay IED	150	8
	Ethernet Switch	50	4
	Merging Unit	150	8
	CB IED	100	8
	Transformer IED	150	8
	Time source	150	4

Table 2 :

 2 Safety Integrity Levels and the Probability of Failures on Demand

	Safety Integrity Levels	1	2	3	4
	Safety Availability	90%-99%	99%-99.9%	99.9%-99.99% Non relevant
	Risk Reduction Factor	10 to 100	100 to 1000	1000 to 10,000 10,000 to 100,000

Table 3

 3

	: Protection and control architectures reliability and
	availability
	Architecture	Reliability %	Availability %
	Basic	99.24185353	99.99512951
	Redundant Ethernet	99.31740880	99.99604270
	Redundant Bay IED	99.34260667	99.99665150
	Transformer protection and control function is a
	SRS function that simultaneously controls line, bus
	bar and transformer bays equipment, including con-
	trolling transformer tap position, CB and discon-
	nectors switches tripping (opening). This function
	does operate in continuous (high mode). Considering
	(table 2), high demand rate requires probability of
	failure rates per hour (PFH) which is the average
	dangerous frequency of failure per hour of compo-
	nents or subsystems.

Table 4 :

 4 probability of failure per hour (PFH) and safety availability of SRS functions for proposed architectures

	Architecture	PFH	Safety Availability %
	Basic	3.8E-06	99.9996
	Redundant Ethernet	2.7E-06	99.9997
	Redundant Bay IED 2.3E-06	99.9998