Yann Guitton

Marielle Schleiche

Fabrice Monteau

Anne-Lise Royer

Nathalie De

La Cotte

Karine Chalvet-Monfray

Olivier Plantard

Lipids profiles characterization of the tick Ixodes ricinus in HPLC-HRMS

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Poster presented at Analytics, Nantes, France, 5-8 september 2022.

Introduction Results

Sample preparation  Several protocols were evaluated and the one below was selected.  Due to the small weight of samples (80 µg to 300 µg) the drying step have been optimized  Grinding can be done with a Precellys

Method

The methodology is now in place and can be applied to new studies

Biomarker identification

The formal identification of the statistically significant ions is not done yet.

Total lipids profiles of individual ticks (especially nymphs) were analyzed by LC-HRMS (QExactive Orbitrap). By Principal Component Analysis based on potential triglycerides signals, we were able to reveal a gradual change in lipid composition between nymphs aged from 0 to 3 months since the last blood meal, validating the proof of concept of this method for estimating tick age. Potential applications of this method will be investigated, as well as its utility for tick-borne disease epidemiology and vector control.

  Conclusion www.laberca.org b INRAE, Oniris, BIOEPAR, 44300 Nantes, France C INRAE, VetAgro Sup, Univ Lyon, EPIA, 69280 Marcy-L'Etoile, France a Oniris, INRAE, LABERCA, 44300 Nantes, France Fax : +33 (0)2 40 68 78 78 -Tél : +33 (0)2 40 68 78 80 laberca@oniris-nantes.fr Yann Guitton a , Marielle Schleicher ab , Fabrice Monteau a , Anne-Lise Royer a , Nathalie Delacotte b , Karine Chalvet-Monfray c , Olivier Plantard b Lipids profiles characterization of the tick Ixodes ricinus in HPLC-HRMS supernatant to an HPLC vial ISTD addition & evaportation Resuspention in ACN:IPA:water (65:30:5) LC-HRMS The tick Ixodes ricinus (Arthropoda; Chelicerata; Arachnida; Acari) is an obligate hematophagous ectoparasite of vertebrates, and the main vector of tick-borne diseases for humans and animals in Europe. The time a tick spends searching for a suitable host can vary from days to years (from 2 to 6 years). Knowing the age of ticks is thus an important data for understanding the epidemiology of tick-borne diseases. Unfortunately, there is currently no appropriate method for determining tick age that has been experimentally validated. Ixodes ricinus female adult (A0:year 2014, A1:year 2015, host: roe deer) and nymph of known age from 0 to 3 month after moulting (N0, N1, N2, N3 year 2015, host: Gerbil) were obtained after their metamorphosis in the laboratory and individualized.

  Developments stage Ixodes ricinus (InVS, 2015, from Parola et al ) Because between each blood meal, the tick uses its lipid reserves while waiting for a new host, in this study we investigated potential triglycerides, a major constituent of tick lipids and a non-renewable energy source, as a marker of aging (i.e., time since last blood meal). Principal Component Analysis of TAG portion of LC-HRMS lipids profils for Adults ((PCA) were applied to describe the patterns observed in the data. For all multivariate statistics, the data wee log10 transformed and Pareto scaled. A total of 52 nymphs (13 for each the 4 groups N0, N1, N2, N3) and 10 adults (5 for each two groups A2 and A14) were analyzed. The PCA was performed on the data acquired in ESI+ mode with R considering the totality of the extracted MS signals in the TAG portion of the chromatogram. With our data processing workflow 1510 variables (features) were extracted for the nymphs and 1595 for the adults. The number of features was finally reduced to 457 and 526 ions using quality filters (reproducibility of the signal in quality controls, chromatographic peak shape...).


G L Richard

J P Vila

Depth-averaged model for open-channel flows in the smooth turbulent case and reconstruction of the 3D-fields

A two-dimensional depth-averaged model is derived for open-channel flows in the smooth turbulent case. The derivation is consistently obtained with a method of matched asymptotic expansions in the outer and inner layers using a mixing length model of turbulence including the free-surface reduction of the eddy viscosity. The shear effects are taken into account by an extra tensor variable of the model called enstrophy tensor. The friction coefficient is an explicit expression of the water depth. The three-dimensional (3D) velocity field and the friction velocity can be reconstructed from the values of the depthaveraged quantities. Numerical simulations show that the enstrophy can be used to evaluate the development of the turbulent boundary layer. In the case of subcritical unsteady flows, the reconstructed velocity can be described with a logarithmic law modified by Coles' wake function with apparent von Kármán constant, integration constant and wake-strength parameter, which differ from their values in steady flows. In the viscous sublayer, the steady-state relation between the velocity and the vertical coordinate, in the inner scaling, is not valid for unsteady flows. Large errors on the calculation of the von Kármán constant can be made if the validity of the steady-state relation is assumed for unsteady flows.

Introduction

The usual equations of open-channel hydraulics in the unsteady case are the one-dimensional (1D) Saint-Venant equations, also called the nonlinear 1 shallow water equations. They can be easily extended to the two-dimensional (2D) case. These equations are derived with the shallow-water assumption i.e. the water depth is small compared to the characteristic length in the direction parallel to the bottom. As a consequence, at the leading order, the pressure is hydrostatic. Another assumption is that all shear effects are neglected, which means that the velocity is supposed to be uniform over the depth. With these assumptions, the Saint-Venant equations can be derived from the Euler equations of incompressible and inviscid fluids with a depthaveraging procedure.

Except in the case of discontinuities, which are created in finite time due to the hyperbolic structure of the equations, there is no inherent dissipative effects in this approach, which implies that they must be added empirically, most often as an empirical friction force. The Kármán-Prandtl relation for smooth pipes can be extended with slightly different numerical values to the case of smooth open channels but the friction coefficient is found only implicitly. Approximate relations were proposed to obtain an explicit expression of the friction coefficient. More details can be found in [START_REF] Chow | Open-channel hydraulics[END_REF] or [START_REF] Yen | Open channel flow resistance[END_REF] for example.

To find the expression of the friction force, or more generally of the dissipative terms, as part of the derivation process of the depth-averaged equations implies to take into account the mean flow and turbulence structure of the flow. Experimental investigation for open-channel flows is more recent than for turbulent boundary layers in close channels because turbulence measurements is more difficult in water than in air flows and it actually started with the advent of Laser Doppler Anemometers [START_REF] Steffler | LDA measurements in open channel[END_REF][START_REF] Nezu | Open-channel flow measurements with a laser Doppler anemometer[END_REF]. The structure of fully developed open-channel flows is similar to boundary layers and pipe flows, with an inner region controlled by the kinematic viscosity ν and by the friction velocity u b = τ b /ρ where τ b is the shear stress at the bottom and ρ the fluid density, and an outer region controlled by the water depth h and the maximum velocity. These regions overlap in a layer where the logarithmic law holds. Denoting by u the mean velocity and by z the vertical coordinate, this log law can be written

u + = 1 κ ln z + + B (1)
where u + = u/u b , z + = zu b /ν, κ is the von Kármán constant and B the integration constant. In the outer layer, a deviation from the log-law can be taken into account by Coles' wake function. However the wake strength parameter Π is smaller than for zero-pressure-gradient boundary layers and is nearly equal to zero at a relatively low Reynolds number [START_REF] Nezu | Open-channel flow measurements with a laser Doppler anemometer[END_REF]. [START_REF] Cardoso | Uniform flow in a smooth open channel[END_REF] found only a weak wake and noted that an apparent log-law can approximate the entire velocity profile. For subcritical flows, the measured values of κ and B are respectively 0.412 and 5.29 [START_REF] Nezu | Open-channel flow measurements with a laser Doppler anemometer[END_REF]. Very close values were obtained by [START_REF] Cardoso | Uniform flow in a smooth open channel[END_REF].

In the case of supercritical flows, the same value of κ 0.41 was measured but it was found that B decreases if the Froude number increases above 1 [START_REF] Tominaga | Velocity profiles in steep openchannel flows[END_REF][START_REF] Prinos | Uniform flow in open channels with steep slopes[END_REF]. [START_REF] Miguntanna | Re-examining log law velocity profile in smooth open channel flows[END_REF] found that the integration constant B is a function of the channel aspect ratio.

In the framework of the eddy viscosity assumption, the mixing length approach (Prandtl 1925) has been extensively applied to open-channel flows and is widely recognized as able to provide an accurate description of the flow over a smooth plane. In the inner layer, the mixing length satisfies the classical linear layer modified by Van Driest damping function [START_REF] Van Driest | On turbulent flow near a wall[END_REF]). An expression of the mixing length for open-channel flows, including the wake strength parameter, was obtained by [START_REF] Nezu | Open-channel flow measurements with a laser Doppler anemometer[END_REF] giving a reduction of the eddy viscosity near the free surface where the mixing length is equal to zero.

In the unsteady case the difficulty is the determination of the friction velocity. Various methods were used and, in particular, the friction velocity can be extracted from velocity measurements in the viscous sublayer assuming the validity for unsteady flows of the law u + = z + which is found in the viscous sublayer in the steady case. The value of the von Kármán constant remains close to κ 0.41 for weakly unsteady flows [START_REF] Nezu | Turbulent structure in unsteady depth-varying open-channel flows[END_REF]) but can deviate from the steady-case value for a strong unsteadiness (Onitsuka & Nezu 2000[START_REF] Nezu | Flow Characteristics in Open-Channel Flows with Strong Unsteadiness[END_REF]. Considerable variations of the integration constant B and also of the wake-strength parameter Π were found.

The present study is a continuation of a previous work (Richard et al. 2017) where a new model for open-channel flows was derived using a mixing length model of turbulence and a method of matched asymptotic expansions. In this article, this work is improved and extended on the following points:

-The mixing length expression of [START_REF] Nezu | Open-channel flow measurements with a laser Doppler anemometer[END_REF] with the freesurface damping effect is used.

-This expression of the mixing length enables an accurate reconstruction of the velocity field from the bottom to the free surface using the calculated depth-averaged quantities.

-The model is extended to the case of 3D-flows, leading to a 2D depthaveraged model.

-The effects of the corrective first-order terms obtained consistently by an asymptotic method are evaluated in unsteady flows with comparisons to experimental results from the literature on the development of the turbulent boundary layer and on unsteady velocity profiles.

The governing equations, the assumptions and the scaling are given in §2. The asymptotic expansions in the outer and inner layers and the matching procedure are presented in §3. The depth-averaged model is consistently derived in §4 using the asymptotic expansions. The method to reconstruct the bottom friction and the 3D velocity fields is given in §5. Numerical simulations are presented in §6 to study the development of the turbulent boundary layer and the velocity profiles in unsteady situations. Technical details are given in appendices A-D.

Governing equations 2.1 Turbulence model

We study a turbulent flow on a sloping channel with a smooth bottom. The angle between the channel and a horizontal plane is θ. The basis for the coordinates x, y and z is (e x , e y , e z ). The angle between the axis Ox and the fall line is α and the axis Oz is normal to the bottom (see figure 1). In these axes the components of the gravity acceleration are g = g (sin θ cos α, sin θ sin α, -cos θ) T .

The turbulence is modelled with the mixing-length model. The viscous stress tensor is written τ = 2ρ (ν + ν T ) D where ρ is the fluid density, ν its kinematic viscosity, ν T the turbulent viscosity. The tensor D is the strainrate tensor defined by D = [grad v + (grad v) T ]/2 where v is the mean velocity field. The turbulent viscosity is given by ν 1956), if the wake-strength parameter is zero or can be neglected, is given by the expression [START_REF] Nezu | Open-channel flow measurements with a laser Doppler anemometer[END_REF],

T = √ 2 L 2
L m = κz 1 -e -z + /A + 1 - z h ( 2
)
where κ is the von Kármán constant (κ 0.41), A + is a dimensionless constant with the usual value A + = 26 and h is the fluid depth. The dimensionless variable z + is the viscous or wall coordinate defined by z + = zu b /ν where the shear or friction velocity u b is related to the bottom shear stress τ b by u b = τ b /ρ. The factor 1 -z/h was absent in the expression of the mixing length used by Richard et al. (2017) and consequently the velocity profile was accurate only in the inner layer. We define the effective viscosity as ν eff = ν + ν T . The constitutive law can thus be written τ = 2ρν eff D.

The wake-strength parameter Π of Coles' law of the wake was found to be considerably smaller in the case of open-channel flows than in the case of zero-pressure-gradient boundary layers where the value Π = 0.55 is observed. [START_REF] Nezu | Open-channel flow measurements with a laser Doppler anemometer[END_REF] found that Π is near zero for Re 10 4 and increases to a maximum of 0.2 for Re 2.5 × 10 4 (our definition of the Reynolds number is Re = hU/ν, different from the definition of Nezu & Rodi 1986). [START_REF] Cardoso | Uniform flow in a smooth open channel[END_REF] found a wake of limited strength (Π 0.08) in the core of the outer region but they found that the wake effect is partly compensated in the near-surface zone by a retarding flow, such that an apparent logarithmic law can approximate the entire velocity profile, explaining why the logarithmic law is often used with success in open channel flow up to the water surface. They also highlighted that the outer region of open-channel flow may not have a universal structure, possibly depending on secondary currents, flow history and inactive turbulence components. Given the small importance of the wake function in open-channel flows and the large increase of complexity needed to take it into account, the wake function is neglected. However we will show in the following that, although no wake function is included in the description of uniform and steady flows, an apparent wake function appears in the unsteady case.

The mass conservation equation in the incompressible case is div v = 0. The components of the velocity field are denoted by v = (u, v, w) T . The components of the viscous stress tensor are denoted by τ xx , τ yy , τ zz , τ xy , τ xz and τ yz and p denotes the pressure. The momentum balance equation is

ρ ∂v ∂t + div (v ⊗ v) = ρg -grad p + div τ (3)
The no-penetration and no-slip conditions at the bottom imply that v(0) = 0. At the free surface, the kinematic boundary condition is

w(h) = ∂h ∂t + u(h) ∂h ∂x + v(h) ∂h ∂y (4)
and the dynamic boundary condition gives the following equations

[p(h) -τ xx (h)] ∂h ∂x -τ xy (h) ∂h ∂y + τ xz (h) = 0 (5) [p(h) -τ yy (h)] ∂h ∂y -τ xy (h) ∂h ∂x + τ yz (h) = 0 (6) -τ xz (h) ∂h ∂x -τ yz (h) ∂h ∂y -p(h) + τ zz (h) = 0 (7)

Shallow-water scaling

The equations are written in dimensionless form using a characteristic depth h 0 , a characteristic length L in the Ox-direction and a characteristic velocity u 0 with the shallow-water hypothesis

ε = h 0 L 1 (8)
The dimensionless quantities are denoted with a prime and are defined as

x = x L ; y = y L ; z = z h 0 ; u = u u 0 ; v = v u 0 ; w = w εu 0 ; p = p ρgh 0 ; t = tu 0 L ; h = h h 0 ; L m = L m κh 0 ; ν T = ν T κ 2 h 0 u 0 ; τ xz = τ xz ρκ 2 u 2 0 ; τ yz = τ yz ρκ 2 u 2 0 ; τ xx = τ xx ερκ 2 u 2 0 ; τ yy = τ yy ερκ 2 u 2 0 ; τ zz = τ zz ερκ 2 u 2 0 ; τ xy = τ xy ερκ 2 u 2 0 ; (9)
A characteristic turbulent viscosity is ν e = κ 2 h 0 u 0 . We define the Froude number F , the Reynolds number Re and the mixing-length Reynolds number Re ML as

F = u 0 √ gh 0 ; Re = h 0 u 0 ν ; Re ML = h 0 u 0 ν e = 1 κ 2 (10)
There is no assumption on the Froude number i.e. F = O(1). We then define the ratio

η = Re ML Re = 1 κ 2 Re = ν ν e ( 11
)
This number is usually very small in open-channel hydraulics. We will assume that

η = ε 2+m , m > 0 (12)
The dimensionless mixing length is L m z √ 1 -s where s = z/h and the effective viscosity is ν T = ν/ν e + ν T = η + ν T which leads to

ν eff = z 2 (1 -s) ∂u ∂z 2 + ∂v ∂z 2 + O(ε 2 ) ( 13
)
We define

λ = sin θ κ 2 F 2 (14)
The molecular viscosity is negligible in this scaling. In this scaling the mass balance equation writes

∂u ∂x + ∂v ∂y + ∂w ∂z = 0 (15)
Defining the 2D-vectors u = (u , v ) T , λ = λ(cos α, sin α) T and τ sh = (τ xz , τ yz ) T , the momentum balance equation in the Oxy-plane becomes

ε κ 2 ∂u ∂t + div (u ⊗ u ) + ∂w u ∂z = λ + ∂τ sh ∂z - ε κ 2 F 2 grad p + O(ε 2 ) (16)
In the Oz-direction the momentum balance can be written

∂p ∂z = -cos θ + O(ε) (17)
The dynamic boundary condition at the free surface (5)- [START_REF] Onitsuka | Unsteadiness effects on inner-layer characteristics in smooth open-channel flows[END_REF] reduce to

p (h) = O(ε) ; τ xz (h) = O(ε 2 ) ; τ yz (h) = O(ε 2 ) (18)
As in Richard et al. (2017), in this scaling the boundary condition at the bottom cannot be satisfied. It is necessary to use another scaling in an inner layer near the bottom wall where the molecular viscosity is included.

Viscous scaling

This scaling is a zoom of the shallow-water scaling using the small parameter η. Dimensionless quantities in this scaling are denoted by a tilde. Some dimensionless quantities are not changed and some other are magnified. We define

The expression of the dimensionless mixing length in the viscous scaling is

Lm = z√ 1 -s 1 -exp - z√ τb κA + (20)
In this scaling the exponential term is not negligible. The dimensionless strain-rate tensor is D = ηD . The dimensionless effective viscosity is νeff = 1 + νT . This implies that the molecular and turbulent viscosities are of the same order of magnitude in this scaling. The mass conservation is not changed and writes

∂ ũ ∂ x + ∂ṽ ∂ ỹ + ∂ w ∂ z = 0 (21)
The momentum balance equation gives

∂ τxz ∂ z = O(η) ; ∂ τyz ∂ z = O(η) ; ∂ p ∂ z = O(η) (22)
3 Asymptotic expansions

The methodology is formally the same as in Noble & Vila (2013) for powerlaw laminar flows and in [START_REF] Richard | A threeequation model for thin films down an inclined plane[END_REF] for laminar Newtonian flows and was detailed in Richard et al. (2017) in the case of two-dimensional flows. This method is extended to the case of three-dimensional flows. Each variable is expanded with respect to the small parameter ε as

X = X 0 + εX 1 + O(ε 2 ) (23)
for any variable X. The expansion of the components of the viscous stress tensor will be denoted as

τ xz = τ (0) xz + ετ (1)
xz + O(ε 2 ). The expressions of the variables are obtained at order zero and then at order one.

Order zero

In the shallow-water scaling, the momentum balance equation [START_REF] Richard | A threeequation model for thin films down an inclined plane[END_REF] gives

∂τ (0) xz ∂z = -λ cos α ; ∂τ (0) yz ∂z = -λ sin α (24)
and the boundary conditions [START_REF] Teshukov | Gas-dynamic analogy for vortex free-boundary flows[END_REF] lead to τ (0) sh (h) = 0. The integration gives τ (0)

sh = λh (1 -s) (25)
The constitutive law τ = 2ν eff D gives τ (0) sh = ν eff ∂u 0 /∂z which leads to

z 2 ∂u 0 ∂z 2 + ∂v 0 ∂z 2 ∂u 0 ∂z = λh (26)
This equation gives the norm

∂u 0 ∂z = √ λh z (27)
The components of u 0 can be integrated between the free surface and an arbitrary depth to obtain

u 0 = u 0 (h) + √ λh λ λ ln s (28)
The expression (25) does not diverge when s → 0 but the expression (28) diverges for s → 0. It is thus necessary to use the viscous scaling to find the expression of the velocity in an inner layer near the bottom. Then a matching procedure will be used in an overlap region to fit the expression of the velocity in the outer layer (with the shallow-water scaling) and in the inner layer (with the viscous scaling).

In the viscous scaling, (22) implies that τ (0) xz and τ (0) yz are constant in the inner layer and thus equal to their values at z = 0. Since τ = τ and because the expressions of τ sh (0) = λh . We have also τb = λh . The constitutive law is integrated in the viscous scaling in order to find the velocity in the inner layer. We define ξ = 2 √ λh z and A = 2κA + . The constitutive law gives

1 + z2 (1 -s) 1 -e -ξ/A 2 ∂ ũ0 ∂ z ∂ ũ0 ∂ z = λh (29)
From this relation, we can show that

∂ ũ0 ∂ z = 2λh 1 + 1 + ξ 2 (1 -s) (1 -e -ξ/A ) 2 (30)
This leads to

∂ ũ0 ∂ z = 2λh 1 + 1 + ξ 2 (1 -s) (1 -e -ξ/A ) 2 (31)
The integration of these equations between the bottom and an arbitrary depth gives

ũ0 = √ λh λ λ - ξ 1 + 1 + ξ 2 + ln ξ + 1 + ξ 2 + R(ξ) (32)
where the function R is defined by

R(ξ) = ξ 0 dξ 1 + 1 + ξ 2 (1 -e -ξ/A ) 2 - ξ 0 dξ 1 + 1 + ξ 2 + O( √ η) (33)
The limit of this function for ξ → ∞ is denoted by R i.e.

R = ∞ 0 dξ 1 + 1 + ξ 2 (1 -e -ξ/A ) 2 - ∞ 0 dξ 1 + 1 + ξ 2 (34)
The vector u 0 in the outer layer and the vector ũ0 in the viscous layer are fitted by the matching procedure. We write that both velocities coincide in an overlap region which is at a very small depth of order √ η written z = √ ηbh where b = O(1). The matching relation can be written

u 0 (s = √ ηb) = ũ0 ξ = 2b √ λh 3 √ η + O( √ η) (35)
The term of O( √ η) is smaller than ε because of [START_REF] Prinos | Uniform flow in open channels with steep slopes[END_REF]. This enables to obtain consistency at order 1. This procedure gives the values of the velocity at the free surface u 0 (h) . The expression of u 0 (h) can be written

u 0 (h) = √ λh λ λ (R -1 + ln 2 + ln M -ln η) + O( √ η) (36)
where

M = 2 √ λh 3 (37)
As in Richard et al. (2017) we introduce the small parameter

µ = - 2 ln η (38)
with ε < µ < 1. The main small parameter ε is smaller than µ p for any positive integer p if ε is small enough (Richard et al. 2017). Each term of the asymptotic expansions of a given order with respect to ε is further expanded in a power series of this second small parameter µ.

The expression (36

) shows that u 0 (h) is of O( √ λ/µ). Assuming that u 0 (h) is of O(1) implies that λ = O(µ 2 ). We write λ = µ 2 λ 0 with λ 0 = O(1).
The expression of u 0 (h) can be written

u 0 (h) = λ 0 h λ λ [2 + µ (R -1 + ln 2 + ln M )] + O( √ η) (39)
This gives the complete expressions of u 0 as

u 0 = λ 0 h λ λ [2 + µ (R -1 + ln 2 + ln M + ln s)] (40)
At order zero, the velocity has the well-known logarithmic profile. In the 1Dcase, reverting to dimensional quantities and introducing the friction velocity, which is u b = √ gh sin θ, the fluid velocity can be written at order zero

u 0 u b = 1 κ ln zu b ν + 1 κ (R -1 + 2 ln 2 + ln κ) (41)
which is the usual log-law (1) with the inner variables u + = u/u b and z + = zu b /ν. The expression of the integration constant B is

B = 1 κ (R -1 + 2 ln 2 + ln κ) (42)
The values κ = 0. Close to the wall, ξ → 0 and ũ0 ∼ ξ √ λh /2. This yields the relation u + = z + , which is valid in the viscous sublayer.

Even if the expressions of u 0 and v 0 diverge for z → 0, they are integrable functions on [0, h] and their depth-averaged values can be calculated. For any quantity X, its depth-averaged value is defined as

X = 1 h h 0 X dz (43)
The depth-averaged velocity at order zero can be calculated from (40). Using the notation U = u = (U, V ) T , we obtain

U 0 = λ 0 h λ λ [2 + µ (R -2 + ln 2 + ln M )] (44)
We define the quantity C(µ) as

C(µ) = U 0 √ λ 0 h cos α = V 0 √ λ 0 h sin α ( 45
)
Its expression is

C(µ) = 2 + µ (R -2 + ln 2 + ln M ) (46)
The velocity in the Oz-direction can be found from the mass conservation equation [START_REF] Richard | Consistent equations for open-channel flows in the smooth turbulent regime with shearing effects[END_REF]. Taking into account the kinematic boundary condition, the integration of (15) leads to

w = ∂h ∂t - ∂ ∂x z h u dz - ∂ ∂y z h v dz (47)
The depth-averaged mass conservation equation is

∂h ∂t + ∂hU ∂x + ∂hV ∂y = 0 (48)
With this equation, the derivative of h with respect to time can be estimated as

∂h ∂t = - ∂h U 0 ∂x - ∂h V 0 ∂y + O(ε) (49)
At order zero, we have

w 0 = - ∂h U 0 ∂x - ∂h V 0 ∂y + ∂ ∂x h 1 s u 0 ds + ∂ ∂y h 1 s v 0 ds ( 50
)
which leads to

w 0 = -s λ 0 h cos α ∂h ∂x + sin α ∂h ∂y 1 + µ 2 (R -1 + ln 2 + ln M + ln s) (51)
The last quantity to calculate at order zero is the pressure. It is found from [START_REF] Steffler | LDA measurements in open channel[END_REF]. The integration is straightforward and gives

p 0 = (h -z ) cos θ (52)
In the inner layer, (22) implies that p0 is constant. The connection with the expression (52) in the outer layer gives simply p0 = h cos θ.

Order 1

In the shallow-water scaling, at order 1, the momentum balance equation becomes

∂u 0 ∂t + div u 0 ⊗ u 0 + ∂w 0 u 0 ∂z = κ 2 ∂τ (1) sh ∂z - 1 F 2 grad p 0 (53)
with the boundary condition τ

sh (h) = 0. The integration of (53) in the Ox-direction, taking into account the boundary conditions leads to

κ 2 τ (1) xz = ∂ ∂t z h u 0 dz + ∂ ∂x z h u 2 0 dz + ∂ ∂y z h u 0 v 0 dz + u 0 w 0 + 1 F 2 ∂ ∂x z h p 0 dz (54)
The calculation in the Ox and Oy directions gives the result

τ (1) sh = λ 0 h κ 2 λ λ λ λ • gradh 1 -s + µT 1 (s) + µ 2 T 2 (s) + h κ 2 gradh cos θ F 2 (s -1) (55)
with

T 1 (s) = (1 -s) (R + 1 + ln 2 + ln M ) + s ln s 2 (56) T 2 (s) = R 2 2 + R - 7 2 + ln 2 R + 1 + ln 2 2 + ln M R + 1 + ln 2 + ln M 2 1 -s 2 + s ln s 2 - 5 2 + R 2 + ln 2 2 + ln M 2 + ln s (57)
At the bottom we obtain

τ (1) sh (0) = λ 0 h κ 2 λ λ λ λ • gradh 1 + µT 1 (0) + µ 2 T 2 (0) - λ 0 h κ 2 gradh cos θ λ 0 F 2 (58) with T 1 (0) = R + 1 + ln 2 + ln M and T 2 (0) = R 2 4 + R 2 - 7 4 + ln 2 2 R + 1 + ln 2 2 + ln M 2 R + 1 + ln 2 + ln M 2 (59)
As for the order zero, the velocity is obtained with the constitutive law τ = 2ν eff D . Technical details are given in Appendix A.

The integration of the constitutive law gives

u 1 = u 1 (h) + √ λ 0 h 2κ 2 λ λ λ λ • gradh × ln s µ +T 1 (0) ln s+ Li 2 (1 -s) 2 +µT 2 (0) ln s+ µ 4 (R -5 + ln 2 + ln M ) Li 2 (1-s) + µ 2 -ln(1 -s) ln 2 s -2 ln sLi 2 (s) + 2Li 3 (s) -2ζ(3) - √ λ 0 h κ 2 cos θ λ 0 F 2 ln s µ gradh - 1 2 λ λ λ λ • gradh (60)
In this expression, Li n denotes the polylogarithm function of order n. For n = 2 the dilogarithm is defined as

Li 2 (s) = - s 0 ln (1 -t) t dt (61)
For n = 3 the trilogarithm can be defined as

Li 3 (s) = s 0 Li 2 (t) t dt (62)
The function ζ is Riemann zeta function and

ζ(3) = Li 3 (1) 1.20 is Apéry's constant.
Because the expression of u 1 above diverges when s → 0, a matching procedure is necessary to connect this expression with the expression in the inner layer found with the viscous scaling. This procedure yields the expression of u 1 (h).

In the viscous scaling, the momentum balance equation ( 22) implies that ∂ τ (1) sh /∂ z = 0. The matching procedure for τ

sh is thus straightforward:

τ (1) sh = τ (1)
sh (0). The velocity field is obtained from the constitutive law. Technical details are given in Appendix B. The integration leads to

ũ1 = √ λ 0 h 2µκ 2 λ λ λ λ • gradh T 0 (0) + µT 1 (0) + µ 2 T 2 (0) ln ξ + 1 + ξ 2 + R 1 (ξ) - √ λ 0 h µκ 2 gradh cos θ λ 0 F 2 R(ξ) - ξ 1 + 1 + ξ 2 + ln ξ + 1 + ξ 2 - √ λ 0 h 2µκ 2 λ λ λ λ • gradh cos θ λ 0 F 2 2ξ 1 + 1 + ξ 2 -ln ξ + 1 + ξ 2 +R 1 (ξ)-2R(ξ) +O(µ 2 ) ( 63
)
where the function R 1 is defined by

R 1 (ξ) = ξ 0 dξ 1 + ξ 2 (1 -e -ξ/A ) 2 - ξ 0 dξ 1 + ξ 2 (64)
The matching procedure follows the same principle as for order zero i.e.

u 1 (s = b √ η) = ũ1 ξ = 2b √ λh 3 √ η + O ( √ η) (65)
We obtain

u 1 (h) = √ λ 0 h 2κ 2 λ λ λ λ • gradh × 2 µ 2 + 1 µ (R 1 + ln 2 + ln M + 2T 1 (0))- π 2 12 +2T 2 (0) +T 1 (0) (R 1 + ln 2 + ln M ) - √ λ 0 h κ 2 cos θ λ 0 F 2 gradh 2 µ 2 + 1 µ (R -1 + ln 2 + ln M ) - √ λ 0 h 2κ 2 cos θ λ 0 F 2 λ λ λ λ • gradh - 2 µ 2 + 1 µ (2 + R 1 -2R -ln 2 -ln M ) +O(µ) (66)
where

R 1 = +∞ 0 dξ 1 + ξ 2 (1 -e -ξ/A ) 2 - +∞ 0 dξ 1 + ξ 2 (67)
is the limit of R 1 (ξ) when ξ → +∞. This procedure yields the complete expression of u 1 . Denoting U 1 = (U 1 , V 1 ) T , the depth-averaged velocity at order 1 is then

U 1 = √ λ 0 h 2κ 2 λ λ λ λ • grad h 2 µ 2 + 1 µ (R 1 -1 + ln 2 + ln M + 2T 1 (0))- 1 2 +2T 2 (0) +T 1 (0) (R 1 -1 + ln 2 + ln M ) - √ λ 0 h κ 2 cos θ λ 0 F 2 grad h 2 µ 2 + 1 µ (R -2 + ln 2 + ln M ) - √ λ 0 h 2κ 2 cos θ λ 0 F 2 λ λ λ λ • grad h - 2 µ 2 + 1 µ (3 + R 1 -2R -ln 2 -ln M ) +O(µ) (68)
4 Depth-averaged equations

Mass and momentum balance equations

The depth-averaged mass conservation equation is given above (48). It can be written in vector form using the two-dimensional divergence operator

∂h ∂t + div (hU ) = 0 (69)
Averaging over the depth the momentum balance equation in dimensionless form in the shallow-water scaling leads to

∂h U ∂t + div (h u ⊗ u ) + grad h 2 2F 2 cos θ = κ 2 ε h λ -τ sh (0) + O(ε)
(70) The expressions (25) at order zero gives τ (0) sh (0) = h λ. The depth-averaged momentum balance equation becomes

∂h U ∂t + div (h u ⊗ u ) + grad h 2 2F 2 cos θ = -κ 2 τ (1) sh (0) + O(ε) (71)
To calculate the term u ⊗ u , we define the tensor

ϕ = (u -U ) ⊗ (u -U ) h 2 (72)
which is conveniently called enstrophy tensor as in Richard et al. (2019) because it has the same dimension as the square of a vorticity. By definition we have the equality

u ⊗ u = U ⊗ U + h 2 ϕ (73)
The depth-averaged momentum balance equation can be written

∂h U ∂t + div h U ⊗ U + h 3 ϕ + grad h 2 2F 2 cos θ = -κ 2 τ
(1)

sh (0) + O(ε) (74) with ϕ = ϕh 2 0 /u 2 0 .
The enstrophy can be expanded as ϕ = ϕ (0) + εϕ (1) + O(ε 2 ). The calculation of

ϕ (0) = 1 h 2 1 0 u 0 -U 0 ⊗ u 0 -U 0 ds (75) yields ϕ (0) = 1 h λ ⊗ λ λ (76)
Writing λ = µ 2 λ 0 where λ 0 is of O(1), the expression of the enstrophy tensor at order zero can be written

ϕ (0) = µ 2 λ 0 h λ ⊗ λ λ 2 = O(µ 2 ) ( 77
)
The expressions at order 1 are found from the integral

ϕ (1) = 1 h 2 1 0 u 0 -U 0 ⊗ u 1 -U 1 + u 1 -U 1 ⊗ u 0 -U 0 ds (78) This gives ϕ (1) = λ 0 κ 2 h λ ⊗ λ λ 2 λ λ • grad h × 1 + cos θ λ 0 F 2 + µ (R + 2 + ln 2 + ln M -ζ(3)) - λ 0 κ 2 h cos θ λ 0 F 2 λ λ ⊗ grad h + grad h ⊗ λ λ + O(µ 2 ) (79)
With all expressions of the asymptotic expansions at order 0 and order 1, τ

sh (0) can be consistently written

τ (1) sh (0) = 1 -α 1 µ C(µ) µ 2 C 2 (µ) U 0 U 1 + U 0 U 0 • U 1 U 0 - α κ 2 κ -α 1 κµ C(µ) κµ C(µ) h 2 λ λ trϕ (1) + α 1 µ C(µ) h 2 ϕ (1) • λ λ + O(µ 3 ) (80)
Quantities of order 1 appear at the right-hand side of this equation. We have

U 0 U 1 + U 0 U 0 • U 1 U 0 = 1 ε U U -U 0 U 0 + O(ε) (81)
and

ϕ (1) = 1 ε ϕ -ϕ (0) + O(ε) (82)
Consequently τ

(1) sh (0) can be written as a sum of relaxation terms as

τ (1)
sh (0) = 1 -α 1 µ C(µ) µ 2 εC 2 (µ) U U - C 2 (µ) µ 2 h λ - α εκ 2 κ -α 1 κµ C(µ) κµ C(µ) h λ λ (h trϕ -λ)+ α 1 ε µ C(µ) h h ϕ • λ λ -λ +O(µ 3 ) (83) with α = R 1 -R + 1 and α 1 = R 1 -R + 1 - 1 2 (ζ(3) -1) (84)
With κ = 0.41 and

A + = 26, R = 2.67, R 1 = 4.82, α = 3.
C f = µ 2 κ 2 C 2 (µ) (85)
The usual Darcy coefficient is f = 8C f . At equilibrium, for a uniform and stationary flow, this coefficient satisfies the implicit relation

1 √ f = ln 10 2κ √ 2 lg Re H f + 1 2κ √ 2 R -2 - 3 2 ln 2 + ln κ (86)
where the Reynolds number Re H defined with the hydraulic diameter is Re H = 4Re. This relation is similar to the Kármán-Prandtl law for pipes flows with smooth surfaces. This inconvenient of this relation is that the friction coefficient is found only implicitly. However, in the general case (i.e. equilibrium or non-equilibrium flows), the relation (85) leads to the explicit relation

κ C f = 2κ √ 2 √ f = R -2 + 2 ln 2 + ln κ + ln gh 3 sin θ ν (87)
With the expression (83) the depth-averaged momentum balance equation ( 74) is obtained in a closed conservative form with source terms, which is

∂hU ∂t +div hU ⊗ U + h 3 ϕ +grad gh 2 2 cos θ = 1 - α 1 κ C f (ĝh -C f U U ) + α κ -α 1 C f h C f ĝ ĝ htr ϕ - ĝ κ 2 -κα 1 h C f hϕ • ĝ ĝ - ĝ κ 2 (88)
where ĝ denotes the projection of the vector g on the plane of the bottom i.e. ĝ = (g sin θ cos α, g sin θ sin α) T and ĝ = g sin θ. It remains to find an evolution equation for the enstrophy tensor.

Enstrophy equation

The momentum balance equation in dimensionless form in the shallow-water scaling can be written

∂u ∂t + div (u ⊗ u ) + ∂w u ∂z + 1 F 2 grad p = κ 2 ε λ + ∂τ sh ∂z + O(ε) (89)
Forming u ⊗(89) + (89)⊗u and averaging the obtained equation over the depth, taking into account the boundary conditions and neglecting all terms of O(µ 3 ) because of the approximation of weakly-sheared flows, leads to the equation of the enstrophy tensor. Details on this derivation are given in Appendix C. The result can be written

∂h ϕ ∂t + div (h ϕ ⊗ U ) -2h ϕ div U + grad U • h ϕ + h ϕ • (grad U ) T = κ 2 ε 1 h 2 U ⊗ τ sh (0) + τ sh (0) ⊗ U -2W + O(µ 3 ) + O(ε) (90)
where W is the dissipation tensor defined by

W = h 0 ν eff ∂u ∂z ⊗ ∂u ∂z dz ( 91
)
The dissipation tensor is expanded as W = W 0 + εW 1 + O(ε 2 ). These asymptotic expansions are given in Appendix C and enable to write the right-hand side of (90) as a sum of relaxation terms. Reverting to dimensional quantities, this equation writes

∂hϕ ∂t + div (hϕ ⊗ U ) -2hϕ div U + grad U • hϕ + hϕ • (grad U ) T = α 2 κ C f h 2 U ⊗ (C f U U -ĝh) + (C f U U -ĝh) ⊗ U -αα 2 C f h U ⊗ ĝ ĝ + ĝ ĝ ⊗ U h tr ϕ - ĝ κ 2 -κα 2 C f h U ⊗ hϕ • ĝ ĝ - ĝ κ 2 + hϕ • ĝ ĝ - ĝ κ 2 ⊗ U (92) where α 2 = 1 2 (ζ(3) -1) (93)
With κ = 0.41 and A + = 26, α 2 = 2.47. Note that α 1 = α -α 2 . The final system of equations is composed of the mass conservation equation ( 69), the momentum balance equation ( 88) and the enstrophy equation (92). Except for the source terms at the right-hand side of these equations, the system has the same mathematical structure as the system derived by [START_REF] Teshukov | Gas-dynamic analogy for vortex free-boundary flows[END_REF] who gave the proof of its hyperbolicity.

Two-dimensional Saint-Venant equations

As implied by (77), the enstrophy is of O(µ 2 ) + O(ε). Furthermore, since

U 1 = O(1/µ 2 ) we can write τ (1) sh (0) = µ 2 C 2 (µ) U 0 U 1 + U 0 U 0 • U 1 U 0 + O(µ) (94)
Consequently the expression (83) of τ

sh (0) can be written

τ (1) sh (0) = µ 2 εC 2 (µ) U U - C 2 (µ) µ 2 h λ + O(µ) (95)
Neglecting terms of O(µ), the dimensional depth-averaged momentum balance equation becomes in dimensional form

∂hU ∂t + div (hU ⊗ U ) + grad gh 2 2 cos θ = ĝh -C f U U (96)
At this level of approximation, there is no enstrophy balance equation and the system reduces to the two-dimensional Saint-Venant equations. The friction term is consistently rather than empirically introduced. Keeping terms up to O(µ 2 ) and neglecting terms of O(µ 3 ) gives the complete system {(69), (88), (92)}.

Energy equation

The system admits an energy balance equation. Taking half the trace of Equation (137) in dimensional form gives the energy balance equation

∂he ∂t + div (heU + U • Π) = 1 - α κ C f ĝh -C f U U • U +α κ -α C f h C f h trϕ - ĝ κ 2 ĝ ĝ •U -καh C f hϕ • ĝ ĝ - ĝ κ 2 •U ( 97
)
where the specific energy e is

e = U • U 2 + gh 2 cos θ + h 2 2 trϕ ( 98
)
and where the tensor Π is

Π = gh 2 2 cos θ I + h 3 ϕ (99)
In this expression, I is the identity tensor. The terms of the right-hand side of the energy equation are relaxation terms due to the dissipative effects in the flow. The expression of the turbulent energy of the system is h 2 trϕ/2. In the particular case of the Saint-Venant equations where the terms of O(µ) are neglected, the specific energy reduces to

e = U • U 2 + gh 2 cos θ (100)
the tensor Π reduces to Π = (gh 2 /2) cos θ I and the energy balance equation reduces to

∂he ∂t + div (heU + U • Π) = ĝh -C f U U • U (101)
5 Reconstruction of the 3D fields

The three-dimensional fields can be reconstructed from the values of the depth h, of the depth-averaged fluid velocity U and of the enstrophy tensor as a function of the applicate z or of s = z/h. The expression of the shear stress at the bottom can be found from the expressions (25) at order zero and (58) at order 1. At order zero, the expression τ (0) sh (0) = λh can be written τ (0)

sh (0) = µ 2 U 0 U 0 /C 2 (µ).
The shear stress at order 1 has been already consistently written as a sum of relaxation terms in (83). Forming τ sh (0) = τ (0) sh (0) + ετ [START_REF] Cardoso | Uniform flow in a smooth open channel[END_REF] sh (0) and reverting to dimensional quantities leads to the expression of the shear stress at the bottom

τ sh (0) ρ = C f U U - α 1 κ C f (C f U U -ĝh) -α κ -α 1 C f C f ĝ ĝ h h trϕ - ĝ κ 2 + κα 1 C f h hϕ • ĝ ĝ - ĝ κ 2 (102)
which is a function of the depth h, the average velocity U and the enstrophy tensor ϕ, with relaxation terms on these quantities but without any derivative. This expression enables to calculate very easily the bottom shear stress with the correction of order 1.

From the expressions (40) and (44), we can write u 0 = U 0 [1+(µ/C(µ))(1+ ln s)]. At order 1, the expressions (66) and (68) lead to

u 1 = U 1 1 + µ C (1 + ln s) -(1 + ln s) α µ 2 C 2 U 1 • λ λ λ λ +α 2 Li 2 (1 -s) + 1 - π 2 6 × µ C U 1 • λ λ λ λ -h trϕ 1 2 trϕ 0 λ λ -α µ 2 C 2 U 1 • λ λ λ λ + O(µ) (103)
Forming u = u (0) + εu (1) and reverting to dimensional quantities gives the 3D-reconstruction of the horizontal velocity field in the outer layer, accurate at order 1,

u = U 1 + C f κ 1 + ln z h -α C f κ 2 C f U • ĝ ĝ -ĝh ĝ ĝ 1 + ln z h +α 2 1 - π 2 6 + Li 2 1 - z h C f κ U • ĝ ĝ -h √ trϕ -α C f κ 2 C f U • ĝ ĝ -ĝh ĝ ĝ (104)
This expression enables to reconstruct the 3D-profile of the horizontal velocity in the outer layer from the quantities h, U and ϕ calculated with the resolution of the 2D-model A similar method is conducted in the inner layer. The expressions of the velocity at order zero and one are given in Appendix D. Note that the expression of ũ1 in the inner layer has to be accurate to within O(µ 2 ) in order to obtain a matching with the expression of u 1 in the outer layer accurate to within O(µ) when ξ → ∞. In dimensional form, the 3D-reconstruction of the horizontal velocity in the inner layer is

u = 1 κ C f U + 1 -2 α κ C f + 2 αα 1 κ 2 C f C f U • ĝ ĝ -ĝh +2α 2 C f C f κ U • ĝ ĝ -h √ trϕ ĝ ĝ R(ξ) - ξ 1 + 1 + ξ 2 + ln ξ + 1 + ξ 2 + 1 κ 1 - α κ C f + αα 1 κ 2 C f C f U • ĝ ĝ -ĝh + α 2 C f C f κ U • ĝ ĝ -h √ trϕ ĝ ĝ × × R 1 (ξ) -2R(ξ) + 2ξ 1 + 1 + ξ 2 -ln ξ + 1 + ξ 2 (105) with ξ = 2κ ν z gh sin θ. ( 106
)
This expression is less convenient than the expression in the outer layer because the functions R and R 1 are not explicit and need to be numerically calculated but the full 3D-velocity profile, from the bottom to the free surface, can be calculated with the depth-averaged quantities h, U and ϕ using the expression (105) in the inner layer and the expression (104) in the outer layer. These expressions connect asymptotically in the overlap layer with an accuracy of O(µ). At the equilibrium, the relaxation terms are equal to zero, and these expressions reduce to

u = U 1 + C f κ 1 + ln z h (107)
in the outer layer, and to

u = C f κ U R(ξ) - ξ 1 + 1 + ξ 2 + ln ξ + 1 + ξ 2 (108)
in the inner layer. In the viscous sublayer, ξ → 0 and the velocity in the inner layer is equivalent to a linear function of z. Defining the friction velocity u b = τ b /ρ, u + = u/u b and z + = zu b /ν and taking τ b = τ sh (0), where τ sh (0) is given by (102), we obtain

u + = z + √ gh sin θ τ sh (0)/ρ C f U + 1 -2 α κ C f + 2 αα 1 κ 2 C f C f U • ĝ ĝ -ĝh + 2α 2 C f C f κ U • ĝ ĝ -h √ trϕ ĝ ĝ (109)
In an equilibrium situation, the relaxation terms are equal to zero, in particular C f U = √ ĝh, and the previous expression gives u + = z + which corresponds to the usual law in the viscous sublayer u + = z + . In a nonequilibrium situation, in particular for an unsteady flow, the relation between u + and z + is still a linear relation but it is more complex and u + /z + = 1.

6 Numerical simulations

Numerical scheme

The system of equations ( 69), ( 88) and ( 92) is a hyperbolic system with relaxation source terms. In the 1D-case it can be written

∂U ∂t + ∂F ∂x = S (110)
where U = [h, hU, he] T , F = [hU, hU 2 + Π, hU e + ΠU ] T and

S =        0 1 - α 1 κ C f (ĝh -C f U |U |) + κ (α -α 1 ) -α 1 C f h C f hϕ - ĝ κ 2 1 - α κ C f (ĝh -C f U |U |) U -α 2 hC f hϕ - ĝ κ 2 U (111)
The energy of the system is e = (U 2 +gh cos θ+h 2 ϕ)/2 and Π = (gh 2 cos θ)/2+ h 3 ϕ. The characteristics of the system are λ 1,2 = U ± gh + 3h 2 ϕ and λ 3 = U .

This system is solved with a classical explicit Godunov-type finite-volume method and a Rusanov Riemann solver. The time step is calculated with a Courant-Friedrichs-Lewy (CFL) condition. At each time step, the enstrophy is calculated from the energy.

The system is solved for the simulation of a subcritical flow in an open channel. At the entrance the discharge is prescribed and the flow is supposed to be non-developed. This means that the enstrophy can be taken equal to zero since the velocity is uniform in the depth. If ϕ = 0 at the entrance, the system reduces to the system of Saint-Venant. The depth in a ghost cell at the entrance is then calculated from the conservation of the Riemann invariant of the Saint-Venant system U -2 √ gh. The end of the channel is treated as a sharp-crested weir as in Richard & Gavrilyuk (2013): if the depth h N in the last cell is smaller than some height d, which corresponds to the height of the weir, then the discharge q N +1 in a ghost cell after the last cell is zero, otherwise q [START_REF] Henderson | Open channel flow[END_REF]. Neumann boundary conditions are taken for the depth and the enstrophy.

N +1 = (2/3)C d [2g (h N -d) 3 ] 1/2 with C d = π/(π +2)+0.08(h N -d)/d

Development of the boundary layer

Simulations are performed for a uniform flow in a steady case, the so-called normal conditions. The value of the Reynolds number Re = hU/ν is chosen. The kinematic viscosity is fixed at ν = 1.0 × 10 -6 m 2 • s -1 . This gives the value q n of the discharge q = hU . The friction coefficient C fn for a uniform and steady flow is then calculated with (86). The value of the Froude number F = U/(gh) 1/2 is chosen and the angle θ is then calculated by sin θ = F 2 C fn in order to have a normal flow. The average velocity of the normal flow is then found by U n = (F 2 gνRe) 1/3 and the normal depth h n is determined as well. The height of the weir d is calculated from the resolution of the equation (gh

3 n sin θ/C fn ) 1/2 = (2/3)C d [2g(h n -d) 3 ] 1/2
, C d being calculated with the normal depth.

The discharge q n is prescribed at the entrance. After a transient regime the system reaches a steady state where the depth and the velocity are h n and U n everywhere except near the beginning of the channel. Since ϕ = 0 is prescribed at the entrance of the channel because the flow is supposed to be non-developed, the enstrophy relaxes towards its equilibrium value. The enstrophy approaches asymptotically its normal value ϕ n = g sin θ/(κ 2 h n ) and, after some distance, the flow is indistinguishable from a normal flow with a fully developed boundary layer. The value of the enstrophy can be used as an evaluation of the development of the boundary layer at the beginning of the channel, with the value ϕ = 0 for a non-developed flow and the value ϕ n for a fully developed flow. Assuming that, for a partially developed boundary layer, the velocity profile satisfies the usual logarithmic law below the boundary layer thickness δ and that it is uniform above this limit up to the free surface, we have and

u 0 = λ 0 h 2 + µ R -1 + ln 2 + ln M + ln z δ if z δ (112)
u 0 = λ 0 h [2 + µ (R -1 + ln 2 + ln M )] if δ z h (113)
Denoting ϕ max = ϕ(δ = h), this yields

ϕ ϕ max = 2 - δ h δ h ( 114
)
As the boundary layer thickness approaches its fully developed value asymptotically, it is difficult to define precisely where the flow becomes fully developed and several definitions were proposed. In our case, the goal is only to check whether the model gives the right order of magnitude of the length of the flow developing zone L i.e. the distance from the entrance of the channel beyond which the flow is fully developed. A reasonable criterion is to take δ/h > 0.99 for a fully developed boundary layer. The value δ/h = 0.99 corresponds to ϕ/ϕ max = 0.9999 according to (114). In the following we use the criterion ϕ/ϕ max > 0.9999 to define a fully developed flow. Numerical simulations were conducted for values of the Reynolds number between 10 4 and 10 6 and for values of the Froude number between 0.1 and 0.8. The case Re = 10 5 and F = 0.5 is presented on Figure 3 where the black curve is ϕ/ϕ max and the red curve is δ/h n calculated from (114), both given as a function of the normalized abscissa along the channel x/h n counted from the entrance. The slope corresponding to these values of the Reynolds and Froude numbers is sin θ = 1/2162 4.6 × 10 -4 (note that tgθ is practically equal (green ); Re = 5 × 10 5 (red ) and Re = 10 6 (black ). The value of L/h n depends mainly on the Reynolds number but weakly on the Froude number. For a given value of Re, it is larger when F becomes close to 1 and slightly larger for very small values of F . For a given value of the Froude number, L/h n increases with the value of Re. The variation of L/h n with the Reynolds number for a Froude number equal to 0.5 is shown on Figure 4(b) in a logarithmic plot. In the case F = 0.5, it is very close to the law L/h n 25.8Re 1/10 (dashed line).

It is very difficult to make comparisons with experimental results due to the fact that the length of the flow developing zone is defined differently, that only relatively small values of the Reynolds number can be studied in laboratory channels and that the channels used in the experiments have a finite width. The goal here is only to check that the order of magnitude of the calculated length L is reasonable. Kırkgöz & Ardıçhoglu (1997) conducted experiments in a smooth channel 0.3 m wide. Due to the relatively small value of the channel width, many experiments are in fact 2D-situations and cannot be considered for a comparison with a 1D-model (the ratio width/depth is as low as 1.50 in an experiment). Therefore only the cases where the ratio width over depth is larger than 4 are considered therafter (and a ratio of 4 is already quite small). The remaining measurements have values of the Reynolds numbers between 7 × 10 3 and 2.1 × 10 4 and values of the Froude number between 0.30 and 0.72. The authors presented the values of L/h as a function of the ratio 4Re/F (or Re H /F where Re H = 4Re) and proposed the empirical law L/h = 76 -0.0001(4Re/F ). We can remark that this law gives obviously wrong results if the Reynolds number is high enough since the predicted value of L/h becomes negative. For values of the Froude number equal to 0.1, 0.5 and 0.8, L/h becomes negative if Re is higher than 19 000, 95 000 and 152 000 respectively. Furthermore the lowest measured values of L/h were found in the cases of narrow channels (when the ratio width/depth is smaller than 4).

The calculated values of L/h n as a function of 4Re/F are presented on Figure 5 (red •) together with the measured values of Kırkgöz & Ardıçhoglu (1997) (black ) and the values calculated from their empirical law (blue ). The range of the Reynolds and Froude numbers values was restricted to the range of the experiments with the largest ratios width/depth (i.e. Re = 10 4 and 0.25 F 0.8) In spite of all above reservations about this comparison, the order of magnitude of the length of the flow developing zone seems reasonable for these values of the Reynolds and Froude numbers.

Case

Re

F sin θ T d (s) (m) q p /q b x (m) C3S30 10 5 0. [START_REF] Kırkgöz | Velocity profiles of developing and developed open channel flow[END_REF]

Unsteady flows

The numerical simulations of unsteady flows are inspired by the experiments of [START_REF] Nezu | Turbulent structure in unsteady depth-varying open-channel flows[END_REF]. The discharge q 0 prescribed at the entrance of the channel is sinusoidal after a delay time t R large enough for the base flow to be in steady-state conditions. At time t R the discharge is increased from a base value q b to a peak value q p after a time T d then decreased to the base value q b after the same duration. Therefore

q 0 = q b if t t R or if t t R +2T d .
Otherwise q 0 is given by

q 0 = q b + q p -q b 2 1 -cos π(t -t R ) T d if t R t t R + 2T d ( 115
)
The flow is supposed to be non-developed at the beginning of the channel. This means that ϕ = 0 is prescribed at the entrance. The flow is studied far enough from the entrance, at an abscissa x, for the flow to be fully developed (x > L). Simulations were performed for values of the Reynolds number equal to 10 4 , 10 5 and 10 6 , values of the Froude number equal to 0.18, 0.5 and 0.8 and values of T d equal to 30 s (a strongly unsteady case), 120 s and 240 s.

In addition the ratio q p /q b was set to 4 for Re = 10 4 and to 5 otherwise and various channel lengths were considered. The base flow is in the normal conditions and this prescribes the value of sin θ. The various parameters of the simulations are gathered in table 1. The variation of the depth h against the average velocity U shows the characteristic loop diagram observed for rivers in flood. The peak velocity appears before the peak depth. The cases C7M24 (low Froude number F = 0.18, weakly unsteady T d = 240 s, black curve), C7M30 (low Froude number F = 0.18, strongly unsteady T d = 30 s, blue curve) and C4M24 (larger Froude number F = 0.8, weakly unsteady T d = 240 s, red curve) are presented on figure 6(a) with Re = 10 5 (see table 1) where the loops are run counter- clockwise. The loop is wider if the flow is more strongly unsteady (C7M30) or if the Froude number increases (C4M24). The evolutions of the depth (black curve), the average velocity (red curve) and the enstrophy (blue curve) are presented on figure 6(b) (C7M24), (c) (C4M24) and (d) (C7M30). Weakly unsteady cases (figures 6(b) and (c) are closer to a kinematic wave with only slight shape changes during the propagation, whereas, in a strongly unsteady case (figure 6(d), the front of the wave steepens with a tendancy to take a sawtooth shape. The evolution of the enstrophy depends on the case: for a small Froude number, the enstrophy increases in the wave (figures 6(b) and (d)) while it decreases for a larger Froude number (figure 6(c)). For intermediate values of the Froude number, the enstrophy increases in the early stages of the wave and then decreases (figure 7(a) for F = 0.5).

The influence of the downstream boundary condition can be important in some cases, particularly during the falling stage of the wave, as there is some reflection on the weir. This phenomenon can lead to a complex behaviour at the end of the falling stage or shortly after that as in the C7M30 case (figure 6(d)). The abscissa x considered to study the flow was usually chosen far from the weir (which is at an abscissa ) but the effects of the interactions with the weir are not trivial as it can be seen from the comparison on figure 7(a) between the cases C3S30 (x = 17 m and = 20 m, solid curves) and C3M30 (x = 17 m and = 30 m, dashed curves) where in both cases Re = 10 5 and F = 0.5. The graphs of h/h n (black curves), U/U n (red curves) and ϕ/ϕ n (blue curves) show that the end of the falling stage is more complex when the distance to the weir is larger. Even the amplitude of the wave is modified by the distance to the weir. Because of this sensitivity to the downstream boundary condition, it is not possible to make precise comparisons to experimental results without a precise knowledge of the weir used at the end of the channel and, more generally, of the precise hydraulic conditions of the experiments, such as the slope.

The reconstruction of the bottom shear stress with the expression (102) is presented on figure 7(b) for the cases C7M24 (black curves), C7M30 (blue curves) and C4M24 (red curves). The solid curves show the ratios of τ sh (0), denoted by τ b , over its value for the base flow, denoted by τ bn , and the dashed curves show the ratios of h/h n . As for the average velocity, the peak value of the bottom shear stress is attained before the peak depth. This is in accordance with the results of [START_REF] Nezu | Turbulent structure in unsteady depth-varying open-channel flows[END_REF]. The difference between these two peaks increases for a strongly unsteady case or for a larger Froude number. In a strongly unsteady case (C7M30), the graph of the bottom shear stress approaches a sawtooth shape and the end of the falling stage is complex due to reflections on the weir.

The reconstruction of the velocity profile is presented on figure 8(a) in the case C7M24 and on figure 8(b) in the case C7M30, with u + = u/u b , z + = zu b /ν, u b being the friction velocity. Even in the case of a relatively weak unsteadiness (C7M24), the velocity profile during the wave is modified with respect to the velocity profile of the steady case (black curve) both in the inner layer and in the outer layer. The red curve shows the velocity profile at the peak depth. In a strongly unsteady case (C7M30), the evolution of the velocity profile is more complex. The black curve is the profile of the steady flow and the green, red and blue curves are the profiles during the early part of the rising stage, at the peak depth and at the end of the falling stage respectively.

These curves can be approximately interpreted with the same laws as in the steady case but the constants have apparent values which are different from the steady-case values. In the outer layer, a log-law is approximately satisfied with an apparent von Kármán constant κ app and an apparent integration constant B app . Furthermore, in many cases, a deviation from this apparent log-law can be interpreted as a wake function as in [START_REF] Coles | The law of the wake in the turbulent boundary layer[END_REF] with an apparent wake-strength parameter Π app . The velocity profiles can be described in the outer layer with the law

u + = 1 κ app ln z + + B app + 2Π app κ app sin 2 πz 2h (116)
Note that this relation is only a convenient description of the actual velocity profile in the outer layer, which is in fact given by (104). The variations of the apparent von Kármán constant is presented on figure 9(a) for the weakly unsteady cases (T d = 240 s) C7M24 (F = 0.18, black ), C3M24 (F = 0.5, blue ) and C4M24 (F = 0.8, red •) and on figure 9(b) for the strongly unsteady cases (T d = 30 s) C7M30 (F = 0.18, black ), C3M30 (F = 0.5, blue ) and C4M30 (F = 0.8, red •). For F = 0.8, the value of the apparent von Kármán constant remains close to the value κ = 0.41 but for smaller Froude numbers, the difference between κ app and κ can be important. For the weakly unsteady cases, the overall evolution is that κ app increases in the rising stage, reaching a value κ app 0.435 for F = 0.18, and decreases in the falling stage. For the strongly unsteady cases, the evolution of κ app is more complex, especially because the interaction with the weir can have a strong effect at the end of the falling stage. The value of κ app increases at the beginning of the rising stage, reaching κ app 0.46 for F = 0.18, The apparent wake-strength parameter Π app is shown on figure 9(e) for the cases C7M24 (black ), C3M24 (blue ) and C4M24 (red •) and on figure 9(f) for the cases C7M30 (black ), C3M30 (blue ) and C4M30 (red •). Its value remains close to 0 for F = 0.8 but is larger for smaller Froude numbers, particularly for F = 0.18 where it reaches a maximum of Π app 0.10 for T d = 240 s and Π app 0.4 for T d = 30 s. In most cases, Π app is positive but it can take negative values at the beginning of the rising stage for F = 0.18 and T d = 30 s or at the end of the falling stage in some cases.

The graphs of figure 9 show that the constants of the apparent law (116) depend strongly on the Froude number. On the contrary, simulations performed for Re = 10 4 , Re = 10 5 and Re = 10 6 at a given Froude number show that these apparent constants depend weakly on the Reynolds number.

In the viscous sublayer, the expression (109) can be used to reconstruct the velocity profile. The ratio u + /z + is calculated with this expression and the results are presented on figure [START_REF] Nezu | Open-channel flow measurements with a laser Doppler anemometer[END_REF] for the cases C7M24 (black ), C7M30 (blue ) and C4M24 (red ). In the steady state, u + /z + = 1 but for unsteady flows, this ratio is smaller. For F = 0.8, u + /z + remains close to 1 but for F = 0.18, u + /z + decreases below 0.8. For a weakly unsteady case (T d = 240 s), this ratio decreases in the rising stage and increases in the falling stage but in a strongly unsteady case (T d = 30 s), the minimum value of u + /z + is reached before the peak depth and there is some further perturbations at the end of the falling stage due to interactions with the weir.

The calculation of the von Kármán constant by Onitsuka & Nezu (2000) and [START_REF] Nezu | Flow Characteristics in Open-Channel Flows with Strong Unsteadiness[END_REF] in unsteady flows used an evaluation of the friction velocity assuming the validity in unsteady situations of the law u + = z + in the viscous sublayer. One of the results of the present work is that this law is not valid in an unsteady case where u + is still a linear function of z + but with u + /z + < 1. As it can be seen on figure 10, the difference can be important, particularly for low Froude numbers, with values of u + /z + as small as 0.7. Consequently, using the relation u + = z + in the viscous sublayer to evaluate the friction velocity in unsteady situations can entail a large error on the calculation of the apparent von Kármán constant and also on the apparent integration constant.

We have calculated the apparent von Kármán constant by using u + = z + to calculate the friction velocity instead of using the relation There is also, particularly for a strongly unsteady case, a sudden increase of κ app near the beginning of the rising stage before a rapid decrease, and this feature was noted by Onitsuka & Nezu (2000). Because of the similarity between our calculations of κ app using u + = z + in the viscous sublayer, even if this relation is not valid in our approach, and the results of Onitsuka & Nezu (2000) and [START_REF] Nezu | Flow Characteristics in Open-Channel Flows with Strong Unsteadiness[END_REF], we think that the qualitative discrepancies between our really predicted values of κ app (figures 9(a) and (b)) and the calculations made from experimental results in the literature are due to the wrong assumption that u + = z + is valid in unsteady situations, which leads to large errors in the evaluation of κ app . The evaluation of the apparent integration constant is also flawed if the validity of u + = z + is assumed in unsteady flows.

Conclusion

A consistent 2D-depth-averaged model for open-channel flows in the smooth turbulent case is derived with a matched asymptotic method and a mixing length model of turbulence including the free surface damping effect but without the wake function. The model can predict accurate velocity profiles in the inner layer and in the outer layer. It can be used in unsteady situations to reconstruct the bottom shear stress and the 3D velocity profile, where the effects of the first-order corrections can be clearly seen. The friction coefficient has an explicit form and can be consistently calculated from the water depth. Shearing effects are taken into account with the variable enstrophy. The development of the turbulent boundary layer can be evaluated from the model's enstrophy. Numerical simulations with the 1D-model show that the predicted length of the flow developing zone has a correct order of magnitude. The ratio of this length over the normal depth increases with the Reynolds number but depends weakly on the Froude number in the case of subcritical flows.

Numerical simulations were conducted for unsteady flows in the subcritical case with a rising stage followed by a falling stage and a sinusoidal hydrograph. The simulations show the characteristic loop diagram observed for river floods. The peak value of the velocity and of the bottom shear stress is attained before the peak depth and the delay of the peak depth is larger for a larger Froude number or for a stronger unsteadiness.

The velocity profile in unsteady flows can be described by an apparent logarithmic law with an apparent von Kármán constant and an apparent integration constant. In many cases, a deviation from this log law can be described by Coles' wake function with an apparent wake-strength parameter. The variations of these apparent constants depend weakly on the Reynolds number but strongly on the Froude number. The variations are large for small Froude numbers and very small for Froude numbers close to 1. The apparent von Kármán constant increases at the beginning of the rising stage and decreases during the falling stage for weakly unsteady flows or before the peak depth for strongly unsteady flows where it can become smaller than the steady value at the end of the falling stage. The variations of the apparent integration constant and of the apparent wake-strength parameter are qualitatively similar.

In the viscous sublayer, the law u + as a function of z + can be studied since the friction velocity can be calculated with the model. It is found that u + in unsteady flows is a linear function of z + but that u + /z + < 1 if the flow is unsteady. The value of u + /z + remains close to 1 for the larger subcritical Froude numbers but can be as small as 0.7 for small Froude numbers. This ratio decreases at the beginning of the rising stage and increases at the peak depth and in the falling stage for a weak unsteadiness or before the peak depth for a strong unsteadiness. Consequently our model predicts that the law u + = z + is not valid in unsteady flows. Assuming the validity of this law to evaluate the friction velocity can lead to large errors in the calculation of the von Kármán constant and of the integration constant in unsteady situations. Indeed, if we assume the validity of this law to calculate the friction velocity instead of using the value predicted by the model, we find values of the apparent von Kármán constant and of the integration constant which are completely different from the consistent predicted values but which are rather similar to the values obtained from experimental measurements by authors who used this method of calculation of the friction velocity.

Acknowledgments

B Integration of the constitutive law in the inner layer

The constitutive law gives

τsh = νeff ∂ ũ ∂ z ( 124
) with νeff = 1 + z2 (1 -s) 1 -e -ξ/A 2 ∂ ũ ∂ z (125)
At order 1, we obtain Because of [START_REF] Prinos | Uniform flow in open channels with steep slopes[END_REF], this expression can be reduced to ∆ 1 + ξ 2 1 -e -ξ/A 2 (132)

τ (1) sh = ∂ ũ1 ∂ z +z 2 (1 -s) 1 -e -ξ A 2 ∂ ũ0 ∂ z ∂ ũ1 ∂ z + ∂ ũ0 ∂ z -1 ∂ ũ0 ∂ z • ∂ ũ1 ∂ z ∂ ũ0 ∂ z (126) We can write ∂ ũ0 ∂ z -1 ∂u 0 ∂z • ∂u 1 ∂z = λ λ • ∂ ũ1 ∂ z (127)
With these relations, (126) becomes

∂ ũ1 ∂ z = 1 √ ∆ λ 0 h κ 2 λ λ λ λ • gradh T 0 (0) + µT 1 (0) + µ 2 T 2 (0) - 2 1 + √ ∆ λ 0 h κ 2 gradh cos θ λ 0 F 2 - 1 - √ ∆ √ ∆ 1 + √ ∆ λ 0 h κ 2 λ λ λ λ • gradh cos θ λ 0 F 2 (133)
The integration from the bottom to an arbitrary depth gives (63). Averaging this equation, taking into account the boundary conditions and the expression (52) of the pressure at order zero, leads to

C Derivation of the enstrophy equation

∂ ∂t (h u ⊗ u )+div (h u ⊗ u ⊗ u )+U ⊗grad h 2 2F 2 cos θ +grad h 2 2F 2 cos θ ⊗U = κ 2 ε (h U ⊗ λ + h λ ⊗ U -2W ) + O(ε) (135)
where the expression of the dissipation tensor W is (91). The averaged quantity u ⊗ u is expressed with the enstrophy tensor in (73) and u ⊗ u ⊗ u can be written

u ⊗ u ⊗ u = U ⊗ U ⊗ U + h 2 U ⊗ ϕ + h 2 ϕ ⊗ U + u * ⊗ U ⊗ u * + u * ⊗ u * ⊗ u * (136)
where u * = u -U . Since u * 0 = O(µ), we have u * ⊗ u * ⊗ u * = O(µ 3 ) + O(ε). All terms of O(µ 3 ) are neglected in the approximation of weakly sheared flows. We obtain

∂ ∂t h U ⊗ U + h 3 ϕ +div h U ⊗ U ⊗ U + h 3 U ⊗ ϕ + h 3 ϕ ⊗ U +div h 3 ϕ ⊗U +h 3 ϕ •(grad U ) T +U ⊗grad h 2 2F 2 cos θ +grad h 2 2F 2 cos θ ⊗U = κ 2 ε (h U ⊗ λ + h λ ⊗ U -2W ) + O(µ 3 ) + O(ε) (137)
Equation ( 70) is written The difference (137)-(139) leads to the evolution equation of the enstrophy tensor

∂h U ∂t + div h U ⊗ U + h 3 ϕ + grad h 2 2F 2 cos θ = κ 2 ε h λ -τ sh (0) + O(ε) (138
∂h ϕ ∂t + div (h ϕ ⊗ U ) -2h ϕ div U + grad U • h ϕ + h ϕ • (grad U ) T = κ 2 ε 1 h 2 U ⊗ τ sh (0) + τ sh (0) ⊗ U -2W + O(µ 3 ) + O(ε) (140)
The direct integration of (91) is not possible but the asymptotic expansion of the dissipation tensor can be calculated with (137). We obtain

W 0 = h 2 U 0 ⊗ λ + λ ⊗ U 0 (141)

Figure 1 :

 1 Figure 1: Definition sketch.

  where the colon denotes the double dot product. For open-channel flows, the mixing length L m , in the model of Prandtl modified by the damping term of Van Driest(

  x = x ; ỹ = y ; t = t ; ũ = u ; ṽ = v ; p = p ; z = xy ; τxz = τ xz ; τyz = τ yz ; τxx = τ xx ; τyy = τ yy ; τzz = τ zz ; ν = ν η ;

  yz do not diverge for z → 0, we have simply τ (0) sh = τ (0)

  41 and A + = 26 give B = 5.28. These values agree with the value B = 5.29 ± 0.47 (and κ = 0.412 ± 0.011) found by Nezu & Rodi (1986) and with the value B = 5.10 ± 0.96 (κ = 0.401 ± 0.017) found by Cardoso et al. (1989). The value of B depends on the value of A + through R. If A + = 26 then R = 2.67. The above values are valid for subcritical flows. For supercritical flows, the value of B can be smaller (Tominaga & Nezu 1992, Prinos & Zeris 1995). This implies smaller values of A + and R. The graphs of R and B as a function of A + are shown on figure 2(a) and (b) respectively. The dashed lines give the case A + = 26 used for subcritical flows.

Figure 2 :

 2 Figure 2: Graphs of R (a) and of the integration constant B (b) as a function of the parameter A + of Van Driest's damping factor. The dashed lines show the case A + = 26.

  [START_REF] Richard | Consistent equations for open-channel flows in the smooth turbulent regime with shearing effects[END_REF] and α 1 = 0.680. In the approximation of weakly sheared flows due to[START_REF] Teshukov | Gas-dynamic analogy for vortex free-boundary flows[END_REF], all terms of O(µ 3 ) are neglected (see Richard et al. 2017 for a complete discussion). The quantity µ 2 κ 2 /C 2 (µ) is important because it is equal to the friction coefficient (Richard et al. 2017), defined by

Figure 3 :

 3 Figure 3: Variations of the normalized enstrophy ϕ/ϕ max (black curve) and of the normalized boundary layer thickness δ/h n (red curve) according to law (114) with the normalized abscissa along the channel x/h n , in the case Re = 10 5 and F = 0.5.

Figure 4 :

 4 Figure 4: (a) Variation of the ratio of the length of the flow developing zone L over the normal depth h n with the Froude number: Re = 10 4 (black •); Re = 5 × 10 4 (blue ); Re = 10 5 (green ); Re = 5 × 10 5 (red ); Re = 10 6 (black ). (b) Variations of L/h n with the Reynolds number for F = 0.5 (dots); dashed line: L/h n = 25.8Re 1/10 .

Figure 5 :

 5 Figure 5: Ratio of the length of the flow developing zone over the normal depth L/h n as a function of 4Re/F : calculated with the model (red •); measured values of Kırkgöz & Ardıçhoglu (1997) (black ); values calculated with the empirical law of Kırkgöz & Ardıçhoglu (1997) (blue ).

Figure 6 :

 6 Figure 6: (a) Normalized depth h/h max as a function of the normalized velocity U/U max : cases C7M24 (black curve), C7M30 (blue curve) and C4M24 (red curve) (see table 1). (b) (c) (d) Evolution of the depth (black curve), velocity (red curve) and enstrophy (blue curve) normalized by their value for the normal flow as a function of normalized time: cases C7M24 (b), C4M24 (c) and C7M30 (d).

Figure 7 :

 7 Figure 7: (a) Comparison between the cases C3S30 (solid curves) and C3M30 (dashed curves): depth (black curves), velocity (red curves) and enstrophy (blue curves) normalized by their normal value as a function of normalized time. (b) Reconstruction of the bottom shear stress normalized by its value in the base flow (solid curves) for the cases C7M24 (black), C7M30 (blue) and C4M24 (red); dashed curves: normalized depth (same colours).

Figure 8 :

 8 Figure 8: reconstruction of the velocity profile. (a) Case C7M24: steady profile (black) and profile at the peak depth (red); (b) Case C7M30: steady profile (black) and profiles at the early part of the rising stage (green), at the peak depth (red) and at the end of the falling stage (blue).

Figure 9 :

 9 Figure 9: Apparent von Kármán constant (a) (b), apparent integration constant (c) (d) and apparent wake-strength parameter (e) (f) as a function of normalized time. (a) (c) (e) Cases C7M24 (black ), C3M24 (blue ) and C4M24 (red •); (b) (d) (f) Cases C7M30 (black ), C3M30 (blue ) and C4M30 (red •).

Figure 10 :

 10 Figure 10: Ratio u + /z + as a function of normalized time for the cases C7M24 (black ), C7M30 (blue ) and C4M24 (red ).

Figure 11 :

 11 Figure 11: Comparison of the apparent von Kármán constant calculated with the friction velocity obtained with (102) (black) or by assuming u + = z + (red): (a) case C7M24; (b) case C7M30.

and τ 2

 2 xz + τ 2 yz = 1 + z2 (1 -s) 1 -e -ξ/A 2 ∂

× 1 +

 1 2z 2 (1 -s) 1 -e -ξ/A 2 ∂ ũ0

2 2F 2 cos θ +grad h 2 2F 2

 222 )Forming U ⊗(138) + (138)⊗U yields ∂ ∂t (h U ⊗ U ) + div h U ⊗ U ⊗ U + h 3 U ⊗ ϕ +div h 3 ϕ ⊗U -grad U •h 3 ϕ +U ⊗grad h cos θ ⊗U = κ 2 ε h U ⊗ λ + h λ ⊗ U -U ⊗ τ sh (0) -τ sh (0) ⊗ U + O(ε) (139)

Table 1 :

 1 Parameters used for the numerical simulations in the unsteady case.

		1/2162	30	20	5	17
	C3M30 10 5 0.5 1/2162	30	30	5	17
	C3M24 10 5 0.5 1/2162	240	30	5	17
	C4M30 10 5 0.8	1/844	30	30	5	17
	C4M24 10 5 0.8	1/844	240	30	5	17
	C7M30 10 5 0.18 1/16680	30	40	5	30
	C7M24 10 5 0.18 1/16680 240	40	5	30

  This work was supported by the AQUA department of INRAE (project Aquanum). Declaration of Interests. The authors report no conflict of interest. A Integration of the constitutive law in the outer layer Using (13), the constitutive relation in the outer layer yields τ sh = z 2 (1 -s) ∂u ∂z

	and finally
	∂u 1 ∂z	=	1 2κ 2	λ 0 h	λ λ	λ λ	• gradh	T 0 (s) + µT 1 (s) + µ 2 T 2 (s) µs (1 -s)	+	cos θ λ 0 F 2	1 µs
												-	1 κ 2	λ 0 h	gradh	cos θ λ 0 F 2	1 µs	(123)
	Integrating this equation gives (60).
												2	+	∂v ∂z	2 ∂u ∂z	(117)
	From this relation, we obtain
									τ 2 xz + τ 2 yz = z 4 (1 -s) 2 ∂u ∂z	4	(118)
	At the order 1, this gives
			2 τ (0) xz τ (1) xz + τ (0) yz τ (1) yz	= 4z 4 (1 -s) 2 ∂u 0 ∂z	2	∂u 0 ∂z	•	∂u 1 ∂z	(119)
	Using the expressions found at order zero and those of τ xz and τ (1)	(1) yz , this
	leads to the relation
	λ λ	•	∂u 1 ∂z	=	1 2κ 2	λ 0 h	λ λ	• gradh	T 0 (s) + µT 1 (s) + µ 2 T 2 (s) µs (1 -s)	-	cos θ µλ 0 F 2 s
												(120)
	Then we can write at order 1
		τ	(1) sh = z 2 (1 -s)	∂u 0 ∂z	-1 ∂u 0 ∂z	•	∂u 1 ∂z	∂u 0 ∂z	+	∂u 0 ∂z	∂u 1 ∂z	(121)
	which gives
						τ	(1)				√	λh	λ λ	λ λ	•	∂u 1 ∂z	+	∂u 1 ∂z	(122)

sh = s (1 -s) h

and

The right-hand side of (137) can be written

Using the asymptotic expansions found above, we can write

Using ( 81) and (82), this expression enables to write the right-hand side of (90) as a sum of relaxation terms.

D Expressions of the zero-order and first-order velocity in the inner layer

In the inner layer, the expressions (32) and (63) lead to

and