Florent Taccone

Germain Antoine

Olivier Delestre

Nicole Goutal

A new criterion for the evaluation of the velocity field for rainfall-runoff modelling using a shallow-water model

. On this test case, measured velocities have been compared with the simulated ones. Thanks to an infiltration source term added to the mass conservation equation, hydrographs at the outlet of the plot have been well reproduced, as well as the main part (including the highest values) of the measured velocities. In particular, the velocities are well represented when the criterion given in the first test case is satisfied. Finally, the model is applied to a real watershed with an area of 1 km 2 on which hydrographs at the outlet have been successfully reproduced for two extreme rainy events [Le Bouteiller et al., 2015]. It appeared that both physical and numerical error were low in the hydrographic network and that the local velocities are correctly reproduced there, unlike in the hillslopes where the velocities are not well reproduced even if the mass fluxes are correct.

Introduction

In mountain areas, civil and industrial safety issues require an assessment of the sediment exportations over some very productive watersheds. As a first step, the comprehension and prediction of the sediment dynamic require a model capable of reproducing hydraulic transfers at the catchment scale, and more particularly to be able to reproduce local velocity field in runoff modelling. Thanks to the increasing computer capacities and the availability of high resolution DEM (Digital Elevation Models), it is becoming more and more relevant to use physically based models. The Shallow Water Equations (SWE) [START_REF] De Saint-Venant | Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit[END_REF] are widely used to represent the overland flow ( [START_REF] Zhang | Modeling of two-dimensional overland flow[END_REF]], [START_REF] Gottardi | A control-volume finite-element model for two-dimensional overland flow[END_REF], [START_REF] Nord | PSEM2D: A physically based model of erosion processes at the plot scale[END_REF], [START_REF] Burguete | Friction term discretization and limitation to preserve stability and conservation in the 1d shallow-water model: Application to unsteady irrigation and river flow[END_REF], [START_REF] Delestre | FullSWOF: A software for overland flow simulation[END_REF], [START_REF] Cea | Experimental validation of a 2D overland flow model using high resolution water depth and velocity data[END_REF]). To solve these equations, two main challenges must be fixed. The first one concerns the numerical aspect of the problem, in particular the treatment of the wet/dry transitions with the presence of steep slopes over a large domain. The need to define a criterion to characterize the accuracy of the numerical resolution is then highlighted. The second challenge is the modelling of an essential physical process: the friction. Regarding the first challenge, many authors have dealt with different aspects of the numerical solution of the 2D Shallow Water models for runoff applications. Since the nineties, the finite volume method ( [Bouchut, 2004], [Toro, 2009], [Vázquez-Cendón, 2015]) is mainly used to solve the 2D Shallow Water Equations. Thus, it appears that the key point of this numerical resolution is to preserve the positivity of the water depths, to deal with dry zones in presence of steep slope, and to be well-balanced in the sense of [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF], meaning preserving the hydrostatic balance of the lake at rest. In the literature, some numerical schemes have been developed to verify all these properties ( [START_REF] Kurganov | A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system[END_REF], [START_REF] Gallardo | On a well-balanced high-order finite volume scheme for shallow water equations with topography and dry areas[END_REF], [START_REF] Canestrelli | Well-balanced high-order centred schemes for non-conservative hyperbolic systems. applications to shallow water equations with fixed and mobile bed[END_REF], [START_REF] Hou | An improved hydrostatic reconstruction method for shallow water model[END_REF], [START_REF] Cea | A simple and efficient unstructured finite volume scheme for solving the shallow water equations in overland flow applications[END_REF], [START_REF] Chertock | Well-balanced positivity preserving central-upwind scheme for the shallow water system with friction terms[END_REF], [START_REF] Xia | An efficient and stable hydrodynamic model with novel source term discretization schemes for overland flow and flood simulations[END_REF]). For the purpose of having a robust resolution of the Shallow Water Equations for runoff applications, the best compromise needs to be found between accuracy of the numerical scheme and the good properties above-mentioned. One common numerical method to respect states of equilibrium from a classical numerical scheme such as the HLLC scheme introduced by [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF] or the kinetic scheme is the hydrostatic reconstruction introduced in [START_REF] Audusse | A fast stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. Nevertheless, in case of steep slopes and low water depths, the hydrostatic reconstruction needs very fine mesh (see [START_REF] Delestre | A limitation of the hydrostatic reconstruction technique for shallow water equations[END_REF]) and this is not compatible with large scale applications like in entire watershed. Several authors propose improvements of the hydrostatic reconstruction ( [START_REF] Caleffi | Well-balanced bottom discontinuities treatment for high-order shallow water equations WENO scheme[END_REF], [Bouchut and De Luna, 2010], [START_REF] Díaz | High order exactly well-balanced numerical methods for shallow water systems[END_REF], [Morales de Luna et al., 2013], [START_REF] Hou | An improved hydrostatic reconstruction method for shallow water model[END_REF], [START_REF] Xia | An efficient and stable hydrodynamic model with novel source term discretization schemes for overland flow and flood simulations[END_REF]). In this paper, a modification of the hydrostatic reconstruction proposed by [START_REF] Chen | A new hydrostatic reconstruction scheme based on subcell reconstructions[END_REF] will be compared to experimental data at different scales from representative test cases of runoff applications and a criterion to qualify the numerical results will be defined. The second challenge concerns the modelling of the friction terms which is one of the essential physical processes for runoff applications. Indeed, in the presence of steep slopes, the solution of the Saint-Venant equations is provided by a slope-friction balance. Some limitations are encountered for the physical representation of important phenomena such as the effect of bottom friction on flow. The diversity of flows from rain-induced runoff leads to varied inundation regimes and the laws of friction are often specific to each regime. To represent the friction of water on the bottom, the model of [Lawrence, 1997] is chosen for its ability to represent different flow regimes. However, since this model is defined for turbulent flows, a Poiseuille friction law is used as in [START_REF] Kirstetter | Modeling rain-driven overland flow: empirical versus analytical friction terms in the shallow water approximation[END_REF] for laminar flows.

The key point of this study is to evaluate the relevance of a numerical model to represent the velocity field in a watershed with rain-induced runoff, using a friction law adapted to each turbulence and inundation regime. More generally, under which conditions the velocity field can be well represented and how much confidence one can have in numerical model for the velocity field representation.

The objectives of this paper are:

• to consider the different aspects of the problem: the numerical treatment of the slope source term and the physical representation of the friction term for application on a wide range of scales;

• to propose a double criterion to qualify the simulation results, and more particularly the relevance of the velocities. This double criterion is defined from both the theoretical error on the discretization of the source term [START_REF] Chen | A new hydrostatic reconstruction scheme based on subcell reconstructions[END_REF] and a heuristic analysis of the root mean square error on the laboratory test cases. The first criterion ensures the mesh convergence of the numerical scheme. The second criterion shows that in the case of low water depths and large slope gradients, it is possible to have well represented velocities in the model;

• to validate and deepen the whole methodology with representative test cases of increasing difficulty from laboratory experiments to real applications.

In this work, the model is confronted with velocity measurements at different scales and the velocity field is evaluated on a watershed application. Many works are interested in the representation of the velocity fields on a specific case ( [START_REF] Tatard | Measurement and modelling of high-resolution flowvelocity data under simulated rainfall on a low-slope sandy soil[END_REF], [START_REF] Cea | Experimental validation of a 2D overland flow model using high resolution water depth and velocity data[END_REF], [START_REF] Kirstetter | Modeling rain-driven overland flow: empirical versus analytical friction terms in the shallow water approximation[END_REF]) and some works that have addressed numerical problems on steep slopes at different scales with a focus on discharges at the domain outlet [START_REF] Kim | Coupled modeling of hydrologic and hydrodynamic processes including overland channel flow[END_REF]. The key contribution of this work is to propose a method for evaluating velocities in the computational domain:

• at different scales, ranging from a laboratory channel of a few meters to a mountain watershed;

• with steep slopes and friction laws adapted to each case;

• with a recent, simple and massively parallel numerical method.

The paper will be organized as follows: in the first section, the governing equations and the finite volume scheme will be described and we will focus on the improvement of the hydrostatic reconstruction by the Chen and Noelle's scheme. The theoretical error due to the discretization of the slope source term, given by [START_REF] Chen | A new hydrostatic reconstruction scheme based on subcell reconstructions[END_REF], is following a dimensionless number depending on the slope, the mesh size and the water depth. The friction modelling will be adressed in details. The second section will be dedicated to laboratory test cases used on one hand for the definition and validation of a heuristic criterion to qualify the simulation results and, on the other hand, the validation of the friction modelling. The first experiment is a laboratory case [START_REF] Kirstetter | Modeling rain-driven overland flow: empirical versus analytical friction terms in the shallow water approximation[END_REF] which simulates a rain induced runoff in a straight channel with a constant slope. It will be used to evaluate the error made on the numerical resolution of the Shallow Water Equations and its effects on the velocities. The second experiment is also a laboratory case, but with a more complex geometry [START_REF] Cea | Experimental validation of a 2D overland flow model using high resolution water depth and velocity data[END_REF]. It is a sinusoidal mould with a rough fixed bottom on which several types of flow are developing. This allows the observation of the velocities for different inundation regimes and consequently the evaluation of the effect of the friction law. Finally, we present the results of the methodology on a more complex test case. Data from an instrumented plot of 40 m 2 in Thies (Senegal) have been selected [START_REF] Tatard | Measurement and modelling of high-resolution flowvelocity data under simulated rainfall on a low-slope sandy soil[END_REF]. They are used to qualify velocities on a real field with rain induced runoff on a larger scale than laboratory experiments. The data of the Draix-Blone observatory in France [Le Bouteiller et al., 2015] are used for an application on the Laval watershed, a well instrumented watershed of about 1 km 2 . Even if the velocities are not measured on this basin, their impact on the outlet discharge and catchment response will be evaluated.

Material and method

Overland flow simulation

The 2D Shallow Water Equations are derived from the Navier-Stokes equations with hydrostatic pressure and uniform vertical velocity assumptions [START_REF] Gerbeau | Derivation of viscous Saint Venant system for laminar shallow water; numerical validation[END_REF]. These assumptions are available if the flow is shallow and the variations of the slope are not large [START_REF] Gerbeau | Derivation of viscous Saint Venant system for laminar shallow water; numerical validation[END_REF]. They write:

             ∂h ∂t + ∂hu ∂x + ∂hv ∂y = R -I ∂hu ∂t + ∂ hu 2 + gh 2 2 ∂x + ∂huv ∂y = gh -∂z ∂x -S f x ∂hv ∂t + ∂huv ∂x + ∂   hv 2 + gh 2 2   ∂y = gh -∂z ∂y -S f y , ( 1
)
where h is the water height, t the time, u the flow velocity in the x direction, v the flow velocity in the y direction, R the rain intensity, I the infiltration rate given by an infiltration model (see section 2.1.3), g the gravity constant, z the bottom elevation and S f x and S f y the friction slope respectively in the x and y directions.

Numerical resolution

An unstructured mesh is used to discretize the domain in space. This allows better adaptation to the complex topography of a basin by better representing its topography and especially the shape of rivers and gullies in the river system. Using the formalism presented in [START_REF] Audusse | A well-balanced positivity "secondorder" scheme for shallow water flows on unstructured meshes[END_REF], a control volume is constructed at each node of the mesh, creating a boundary between each neighbor node of the mesh, as shown in figure 1.

Figure 1: Representation of the formalism used for the 2D unstructured mesh resolution.

The Shallow Water Equations (1) can be solved with a finite volume scheme which writes:

U t+1 i = U t i - j∈Vi ∆t ∆x ij F t ij + ∆t ∆x ij S t ij , (2)
where

U t i = (h t i , h t i u t i , h t i v t i
) the states at time t in the node P i , V i all the nodes surrounding P i , ∆t the time step, ∆x ij = Ci Lij the space step with C i the area of the cell constructed around the node P i and L ij the length of the boundary edge between the cells C i and C j , F t ij is the numerical flux at the interface between the cells C i and C j along n ij , the normal vector to the boundary edge and S t ij = (R i -I i , s t ijx , s t ijy ) are the source terms of mass and of momentum along x and y directions. The main difficulty in solving these equations lies in how to deal with the slope source terms s t ijx and s t ijy to respect the positivity of the water depth and the hydrostatic balance and more particularly in case of wet/dry transitions and steep slopes applications.

The Shallow Water Equations are solved with the scheme (2) and the interface fluxes are calculated with the HLLC method [START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF]. The intermediate states involved in the flux calculation are defined as

U ij = (h ij , h ij u i , h ij v i ).
The modified hydrostatic reconstruction method finally gives:

z ij = min(max(z i , z j ), min(h i + z i , h j + z j )), h ij = min(h i + z i -z ij , h i ), S ij = g 2 (h i + h ij )(z i -z ij )n ij L ij .
Based on the definition given by [START_REF] Chen | A new hydrostatic reconstruction scheme based on subcell reconstructions[END_REF], a fully wet state when the water depth is higher than the bottom difference between two cells, and a partially wet regime otherwise, are distinguished at the interface of two cells (Figure 2 illustrates these regimes).

Partially wet Fully wet

Figure 2: Illustration of the partially wet and fully wet regimes at the interface of two cells.

If we consider the non dimensional number β = h ∆x∂xz , if β is inferior to 1 the regime is partially wet and fully wet otherwise. The exact error made on the source term then can be estimated. Indeed, [START_REF] Chen | A new hydrostatic reconstruction scheme based on subcell reconstructions[END_REF] give the theoretical slope source term error like:

S ij S ex = 1 - 1 2β , (3)
for the fully wet regime (β ≥ 1) and

S ij S ex = 1 - β 2 , (4)
for the partially wet case (β < 1). The simulations are performed in two dimensions with the open source computational code TELEMAC2D [Hervouet, 1999] mainly used in river applications. The computation code is parallelized with the MPI method. The computation times given in the next sections correspond to simulations on 56 cpu cores on Intel Xeon CPUs E5-2697 v3 -2.60GHz. The presented scheme has been added the TELEMAC2D computation code and is available since the V8P0 version.

Friction term

Whatever the formulation considered, the friction term is computed with a semi-implicit treatment [START_REF] Kirstetter | Modeling rain-driven overland flow: empirical versus analytical friction terms in the shallow water approximation[END_REF]. [START_REF] Bristeau | Boundary conditions for the shallow water equations solved by kinetic schemes[END_REF]] demonstrated that the friction term S f defined as:

S f = 3νu gh 2 , ( 5
)
where ν is the water kinematic viscosity, is the approach that fits best with the data for the flume test case (section 2.2.1). More generally, this formula is used to represent friction in cases where the flow is laminar. For turbulent flows, the friction term defined by [Lawrence, 1997] is introduced to spatialize the friction term according to a constant size of roughness k s . This roughness size is linked to the Darcy-Weisbach coefficient with the following relation:

f =                      1 1.64 + 0.803ln(Λ) 2 if Λ ≥ 10 10 Λ 2 if 1 ≤ Λ ≤ 10 8 π C D min π 4 , Λ else , (6)
with Λ = h ks the inundation rate and C D the drag force of the rough particles fixed to 1 in [Lawrence, 1997]. With this formulation, the friction coefficient varies continuously for different inundation ratios and is related with the friction slope with the relation:

S f = f 8g u|u| h . (7)
It has been demonstrated in [Lawrence, 1997] that this friction model can fail for Reynolds numbers below 1000 and inundation ratios above 10 (Λ ≥ 10 in the equation ( 6)). In these cases, as the Reynolds number is low, the flow is laminar. By replacing the Lawrence's friction model by the Poiseuille model, presented in the equation ( 5), a mixed friction model is introduced.

All these friction models are not yet available in the official version of TELEMAC2D and have been implemented for this work.

Infiltration model

An improvement of the [START_REF] Green | Studies on soil physics[END_REF] model, presented in [START_REF] Brakensiek | Agricultural management effects on soil water processes part II: Green and Ampt parameters for crusting soils[END_REF]] and [START_REF] Esteves | Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values[END_REF] is used to represent the infiltration in the model. It is a vertical 1D model computed at each cell of the domain. As this infiltration model is only 1D vertical, its resolution is much faster than a 3D model with groundwater exchange like [Richards, 1931] type models. Moreover, the conductivity of the layers can be variable in the domain. The Horton model [Horton, 1933] is also widely used, but it does not allow to represent the vertical structure of the soil like this model. The infiltration rate is described like:

I = K 1 + h f + h z f , (8)
with K the soil conductivity under less than 1 cm of hydraulic head, h f the capillarity head and z f the wetting front. Then, the wetting front is updated following the equation:

z f = I C θ s -θ i , (9)
with I C the cumulated height infiltrated since the beginning of the event, θ s the saturated soil moisture or the porosity and θ i the initial soil moisture. In the watershed application, the soil is vertically divided in two layers, a first one with an associated depth Z c and conductivity K c and a second one with a depth considered as infinite and a conductivity K s . The conductivity K varies with the wetting front evolution according to this equation:

   K = K c if z f ≤ Z c K = z f z f -Zc Ks + Zc Kc if z f > Z c . (10)
This model can be used with one or two layers depending on the nature of the soil in each case. Considering the thickness of the first layer as infinite allows to reduce it to a single-layer model. This single-layer model is used in particular in the case of the Thies instrumented plot. The use of the two-layer model is justified by the fact that the measurements carried out in the Laval watershed clearly show a different structure between the surface layer and the inner layer of the soil. Indeed, [Mathys, 2006] shows that for out of more than 160 profiles measured in the basin, there is a superficial structure, unstructured and porous of 8 cm deep on average. Deeper, there is a much more structured layer with less porosity. This model is not present in the TELEMAC2D official sources and has been added by the authors via users routine.

Experimental test cases

Four test cases are selected at different scales and degrees of complexity. First two laboratory test cases are selected with a simple geometry. The first laboratory case is used to evaluate the effect of the resolution method and the other one to evaluate the impact of the law of friction on local velocities. Then, the model is applied to a real case on a small plot to confirm the results of the two previous test cases. Finally, the model is evaluated over a watershed of about 1 km 2 on which it is possible to reproduce the discharge at the outlet.

Flume experiment

[ [START_REF] Kirstetter | Modeling rain-driven overland flow: empirical versus analytical friction terms in the shallow water approximation[END_REF] provide data for a test case which is a straight channel of dimensions 4.04 m by 0.115 m. With a rain simulator over the domain, a steady rain, with an intensity of 25 mm/h or 50 mm/h depending on the chosen scenario, is applied during 600 s. The temporal evolution of the discharge at the outlet of this channel is observed with a time step of 0.1 s. This temporal evolution of the discharge is obtained by weighting the water flowing out of the channel. The water depth and the velocities are also measured along the channel every 60 cm when the steady state is reached. They are measured with a salt-tracing method [START_REF] Planchon | An automated salt-tracing gauge for flow-velocity measurement[END_REF]. The dataset is available for three slope values: 2%, 5% and 25%. All slope and rain intensity configurations are used with the exception of the case with a slope of 5% and a rain intensity of 50 mm/h as no data are available on this case. The table 1

Sinusoidal mould experiment

The second test case is a laboratory test case representing a real field with several values of inundation ratio [START_REF] Cea | Experimental validation of a 2D overland flow model using high resolution water depth and velocity data[END_REF]. At the inlet, runoff with water depths from 1 to 3 mm is observed, with high velocities. Then the flow is directed to furrows where the water depths can reach 2 cm. This test case is representative of the hydraulic transition which can be observed on a real plot with hydraulic jump between the runoff and the flow in the furrows. The number of furrows varies from one to three depending on the choice of the slope and the inlet discharge. Three configurations are considered. The basic configuration has an inflow of 94.12 cm 3 /s, a mean slope of 11.54 %, so we have only the first furrow filled with water. In the second configuration, the mean slope is increased to 20.4 %, the inlet discharge stays close to the initial configuration 108.52 cm 3 /s and water is flowing in the three furrows. The last configuration takes the same slope as the first configuration 11.54 %, but the inflow is three times higher: 291.15 cm 3 /s, and two furrows are supplied by water. The data available on this test case are measured at the steady state so they do not evolve in time. The water depths and the bottom elevation are measured with a laser scanner as presented in [START_REF] Legout | High spatial resolution mapping of surface velocities and depths for shallow overland flow[END_REF] in the whole domain with a spatial resolution of 0.5 mm and the experimental uncertainty of the water depth is 1 mm. Concerning the velocities, they are measured with the LSPIV method [START_REF] Legout | High spatial resolution mapping of surface velocities and depths for shallow overland flow[END_REF], with a spatial resolution of 5 mm and an experimental uncertainty of 0.1 m/s. Figure 3 shows the bottom topography of the sinusoidal mould.

There is no infiltration in this test case.

Figure 3: Bottom elevation of the sinusoidal mould experiment.

Plot experiment

The third experiment is described in [START_REF] Tatard | Measurement and modelling of high-resolution flowvelocity data under simulated rainfall on a low-slope sandy soil[END_REF]. The measurement site is located in Thies, Senegal. The plot is 10 m long and 4 m wide. An uniform rain is applied to the domain using a rain simulator as described in [START_REF] Esteves | Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values[END_REF]. After several rain events that have eroded the soil, the sediment mobility of the soil was considered as negligible and the soil is considered fixed for the period of interest. The soil is mainly composed of sand with a roughness size that will be considered spatially uniform. The average rain intensity applied by the simulator is 67 mm/h for 2 h. We note that there is a significant infiltration process in the domain because the total mass of rainfall does not correspond to the mass that was collected at the outlet. Velocities measurements were taken at 62 different locations on the plot. The measurement procedure is salt-tracing as in the first experimental case. Unlike the previous test case, there is no water level measurement in this case. Figure 4 shows the bottom topography of the plot and the location of the 62 velocity measurements. Some irregularities can be observed on the bottom topography which can mainly be due to measurement errors. However, since the true elevation of the bottom at these locations is not known, and they are quite far from each measurement point, the topography is kept as it is for the simulations. Overland flow models have already been tested on this dataset [Mugler et al., 2011].

Figure 4: Bottom elevation and velocity measurement points of the plot experiment in Thies.

Watershed application

The model is then evaluated by comparing results of the model to measured outlet discharges from field measurements on a real catchment. The Laval watershed is a sub-catchment of the Bouinenc watershed, located on the Draix-Bleone site [Le [START_REF] Bouteiller | Observatoire hydrosédimentaire de montagne Draix-Bléone[END_REF] in the Southern French Alps. The climate is a mountainous climate influenced by the Mediterranean climate which is reflected in particular by heavy summer storms and an average annual rainfall of 900 mm. Its total area is 0.86 km 2 and the mean slope is 58 %. The soil is mostly constituted of black marls and 68% of the surface is a bare soil. At the outlet of the catchment, the data are available for many rainfall events with a time step of 60 s. The rainfall for every event is also measured every 60 s. The flow rate is evaluated by laser measurement of the water depth over a known fixed section. The size of the basin and the steep slopes cause the hydrological response to be fast. Consequently, the short and intense storm events lead to the highest flow rates at the outlet [Mathys, 2006]. Two different events are selected, one spring rain following three rainy days with a high soil moisture initial condition, and a summer storm with a dry initial soil state. The choice of these events is motivated by the fact that the initial state of the soil is very different from one to the other, but also because they are the two events that had the greatest impact on sediment transport in a year. Figure 5 shows the bottom topography of the Laval watershed. Different studies have been performed in this observatory, mainly related to sediment transport ( [START_REF] Lukey | Test of the SHETRAN technology for modelling the impact of reforestation on badlands runoff and sediment yield at Draix, France[END_REF], [START_REF] Bechet | Detection of seasonal cycles of erosion processes in a black marl gully from a time series of high-resolution digital elevation models (DEMs)[END_REF]), but also to hydrological processes [START_REF] Esteves | Soil hydraulic properties in a marly gully catchment (Draix, France)[END_REF]. 3 Results

Flume experiment

The domain has been discretized in three triangular meshes with a space step of respectively 10 cm, 1 cm and 1 mm. Thus, the outlet discharge considered is rescaled for a width of 0.115 m. On the 1 mm mesh, the number of elements is 1.3 million and the computation time is about 3 hours. The simulation duration is 1000 s, the first 600 s representing the rise of the hydrograph and the steady state and the last 400 s the emptying of the domain.

To evaluate the error between the simulation and the measurements, the root mean square error is calculated on the outlet discharge. Then, it is normalized by the mean values of the measured outlet discharge at the steady states. The error is first evaluated on the outlet discharge, because it is the most accurate data available for this test case. Indeed, velocities are measured at only seven points along the channel with uncertainties regularly exceeding 10%.

In this test case, the slope source terms treatment is crucial in the momentum conservation equation because the water depths are low and the slopes are steep. In addition, the equations ( 3) and (4) show that the numerical accuracy of the slope source terms treatment is strongly dependent on the mesh size (via the definition of β). The evolution of the water depth is regular in this test case, so the mean value in the domain at the steady state is considered as the representative depth for the calculation of β. In figure 6, the simulations error made on the outlet discharge is compared with the theoretical error made on the source term depending on the value of β.

For the fully wet regime (β > 1), the RMSE on the outlet discharge follows the theoretical source term error. This result implies that we have a convergence in mesh size because whatever the slope and the positive water height defined, a mesh size can always be adapted to be in the case where β is large enough to minimize the error. But due to CPU limitation β could be less than 1 so the partially wet regime must be studied. There is no clear relation between the simulation error and the theoretical source term error. It appears that an error is made on the homogeneous problem, coming from the reconstructed values used in the calculation of the numerical fluxes. Focusing on the points that are in the partially-wet regime, it shows that the six points corresponding to the smallest errors are the simulations with a 25% slope. In addition, the two simulations corresponding to the smallest errors are the simulations with a rain intensity of 50 mm/h, and the error is very close for two different mesh sizes. Starting from these observations, it seems that the error made in the partially-wet case is not strongly dependent on the mesh size. Moreover, it seems to decrease when the slope and the water depths increase. By plotting the simulation error and the value of the water depth multiplied by the slope, a clear relation is highlighted in figure 7. This result suggests that for the partially wet case, the error is decreasing when the slope source term gh∂ x z is increasing.

To illustrate the impact of the numerical resolution on the velocities behavior in the domain, three simulations in particular are selected. The first case is the case where the error is minimal, i.e. the fully-wet case for which β is maximal. This is the case with a 50 mm/h rain on a 2% slope with a mesh size of 1 mm. Figure 8 shows the simulated discharge compared to the experimental data at the channel outlet, but also the velocities along the channel at the steady state. This result shows that for a fully-wet case with a large enough value of β, the numerical resolution is accurate and the simulations fit well with the experimental data. The second case chosen is where the error is maximal. This corresponds to a rainfall of 25 mm/h on a slope of 2% with a mesh size of 10 cm. In this case, a loss of quantity of movement is observed and the outlet discharge is very underestimated, as shown on figure 8. The values of the mass fluxes at the steady part of the hydrograph show that the numerical scheme is mass conservative. This corresponds to the sum of the F ij terms in the equation (2) at the output boundary condition. The simulated discharge presented here is the product between the longitudinal component of the velocity and the water depth (hu in the shallow water equations ( 1)). Therefore, the difference between the measured discharge and the experimental signal is directly related to the treatment of the slope source term in the momentum conservation equation. This also has an effect on the velocities, which are largely underestimated. The third case chosen is the case where the error is minimal when β is less than 1, when the term h∂ x z is maximal. This is the case with a 25% slope, a rainfall of 50 mm/h and a mesh size of 10 cm. A discharge lower than the experimental discharge is observed at the outlet of the channel. However, the discharge obtained is greater than the one obtained in the previous case. Indeed, to compensate for the loss of motion quantity, the model slightly overestimates the velocities, as seen in figure 8.

To sum up, velocities in the domain are well calculated when β is large enough. The criterion chosen to consider the error acceptable (RMSE<25%) is β > 2 for the fully wet case. However, this is not always possible, especially for flows with low water depths on steep slopes to construct a mesh fine enough to always be in this case. When β is lower than 1, the velocities can be significantly underestimated. Nevertheless, this only happens when the slope source term gh∂ x z is small, and the velocities are a bit overestimated if it is large. The criteria chosen to consider the error acceptable (RMSE<25%) is gh∂ x z > 2 • 10 -4 m 2 /s 2 , and it is independent from the mesh size on the partially wet case.

Sinusoidal mould experiment

In the previous section, a double criterion has been defined on a simple test case to evaluate the runoff resolution error and its impact on velocities in the domain. By evaluating the velocities under these criteria, which are β > 2 or gh∂ x z > 2.10 -4 m 2 /s 2 , this second test case allows to study several types of inundation regime. Thanks to its sinusoidal geometry, there are furrows in which the flow develops and between these furrows there is runoff with small water depths. The model is calibrated with the water depths, because it is the most accurate and a finely spatialized data available in this experiment. This measurement at the steady state is used to calibrate the representative size of the roughness in the friction law.

To choose the parameter of the roughness size k s in the [Lawrence, 1997]'s friction model ( 6), several simulations have been made on each configuration. As the uncertainty measurement on the water depth is low, we choose to minimize the RMSE on this parameter. The minimal value of the RMSE is reached for a roughness size between 0.25 mm and 0.27 mm depending on the configuration. The value of k s is then fixed to 0.26 mm. This value is coherent as the soil is made of grains between 0.5 mm and 1 mm and the roughness visible in the topography plays a role in slowing the flow. Moreover, the mesh sizes and the topographic measurements are fine enough to represent a part of these grains.

Accordingly to the flume test case, we created three triangular meshes with different grid sizes. The mesh sizes chosen are 10 mm, 5 mm and 2 mm for each configuration. For the finest mesh configuration, the number of elements is 536 654 and the computational time for 1000 s of simulation is about 7 minutes. The water depth in the furrow is well represented by the simulations for each tested configuration, as it is shown in figure 9. The location of the hydraulic jump and the flow in the furrows are in good agreement with the experiments. Indeed, the maximum water depths, as well as the location of the hydraulic jumps at the entrance of the furrows are very well reproduced. However, the water level is a little overestimated in the farthest furrows, especially in the experiments with three furrows. As shown in [START_REF] Cea | Experimental validation of a 2D overland flow model using high resolution water depth and velocity data[END_REF], the velocities measurement method can give inaccurate results where the water depths are small. In this test case, the main interest is in the physical representation of friction, so a filter is applied to the velocities and only those for which the numerical resolution is good are considered. Consequently, only velocities measurement points for which β is greater than 2 are therefore taken into account. The table 2 shows that the RMSE is smaller with the finest mesh, but the degradation of the spatial discretization does not impact significantly the accuracy. As said previously, the velocity measurement points on the domain are chosen under the criteria given in 3.1, which means β is greater than 2, the slope being too gentle to be in the case where β < 1 and gh∂ x z > 2.10 -4 m 2 /s 2 . It is therefore considered that the resolution error is negligible and that friction formula plays a major role in the velocities description. Figure 10 shows that the smallest velocities are well represented in the model, but the highest velocities are underestimated with the [Lawrence, 1997]'s friction model. The trend of the curve shows that the modeled velocities are limited. Indeed, high velocities do not seem to exceed a threshold for each simulation. The second configuration is also presented in [START_REF] Cea | Experimental validation of a 2D overland flow model using high resolution water depth and velocity data[END_REF] and the results are very similar.

By studying the distribution of Reynolds numbers in the domain, it appears that a velocity limitation often occurs for relatively small Reynolds numbers. In addition, figure 10 shows that there is no more limitation for the highest velocities with the mixed friction model. The roughness coefficient k s of the Lawrence's model has been recalculated in function of the water depth by minizing the RMSE as before. The new value found is 0.3 mm. The highest velocities are better represented by the model, even some unexplained underestimations remain for the test cases 1 and 2. The Lawrence's friction model is more accurate to represent the smallest velocities in the domain, but it is necessary to use the mixed friction model to represent the highest values of the velocities. In the perspective of a coupling with a erosion and deposition sedimentary model, these velocities are the most important contributors.

Plot experiment

This experiment has been chosen to evaluate whether the conclusions drawn from the laboratory experiments could be exported to a real plot. It also acts as a transition between laboratory cases and a real watershed, especially with the introduction of the infiltration model described in 2.1.3. From a given average rainfall, the infiltration model is calibrated to reproduce the flow at the outlet of the plot. Then the velocities of the model are compared to measured velocities at several points. The error induced by the resolution of the shallow water is first evaluated on the velocities measurement points according to the conclusions drawn on the flume test case. Then, the goal is to see if the same limitation as on the sinusoidal mould test case is observed and if it is corrected by the mixed friction model.

On the plot, the finest mesh created is a 25 mm mesh, which corresponds to the resolution of the topographic data. This mesh has 137 921 elements and the computational time for the 8000 s simulation is about 7 minutes. The rain applied on the plot is variable with a mean value of 66.73 mm/h. The rain intensity is measured with rain gauges and fluctuations are observed around the mean value. These fluctuations are reflected in the outlet discharge observations. However, the salt-tracing velocity measurement technique requires a steady state in the plot. These fluctuations therefore lead to average several velocity measures spread over time. For the simulations, the choice made is to apply a constant rain to compare measured velocities with steady state simulated velocities. Thus, the infiltration model is calibrated to reproduce the total amount of water exiting the domain during the rain event. Figure 11 shows the measured and simulated total rain applied to the plot and the induced discharges at the outlet. To obtain this discharge with the Green-Ampt infiltration model, only one layer of soil is considered. The parameters used are:

• the layer conductivity: K = 7 mm/h,

• the effective porosity: θ s -θ i = 0.125.

The behavior of the hydrograph is well reproduced with this infiltration model. The raise of the hydrograph is very coherent with the observations. However, the hydrograph is drying up a little bit more rapidly than that observed, which can be explained by the fact that exfiltration is not represented by the model. The simulated velocities are compared with the 62 measurement points when the flow rate is stabilized (t=7000s). To set the representative size of the bottom roughness, k s in the friction model, the RMSE on the velocities is minimized. Unlike the previous experiment where the friction was set on the water depths, here only the measured velocities are available in the dataset. The parameter k s minimizing the RMSE is 0.37 mm for both Lawrence's friction and mixed friction model. Figure 12 represents the simulated velocities compared to the measured velocities in the domain. For both friction models, the highest velocities are well represented and an underestimation is observed for a part of the lowest velocities. The difference between the two friction models is observed only on the two points corresponding to the highest velocities.

Based on the results of the flume test case, the velocity can be underestimated for low values of β. This underestimation is observed for cases where β is less than 2 and the slope source term is not dominant (gh∂ x z < 2.10 -4 ). In this experiment, the slope is not steep enough to be in the case where β < 1 and gh∂ x z > 2.10 -4 so we can distinguish two cases: the case where β ≤ 2 and the velocity is underestimated because of the numerical scheme and the case where β > 2 and the velocity is well represented. These two cases are highlighted in figure 12 and explain the part where low velocities are underestimated. A mesh convergence analysis has been performed on this test case and the results are displayed in the table 3. It shows that the RMSE decreases with the mesh size. In addition, the topography is less well represented and the slopes are not strong enough for the source term to become dominant. To conclude, this experiment confirms the relevance of the use of the simplified Green-Ampt infiltration model to estimate the outlet discharge. The velocity behavior corroborates the results obtained from the two laboratory experiments. Velocities are underestimated when β is less than 2 due to the limitations of the numerical scheme highlighted in the flume test case. The two highest velocities are underestimated if the Lawrence model is located in the well inundated case with Reynolds numbers below 1000, which is corrected by the mixed friction model as shown in the sinusoidal mould test case. However, the larger is the mesh size, the smaller is the difference between the models. The choice of the mixed friction model is thus relevant for small scale experiment, especially to represent the highest velocities, but its impact is limited when the scale and the mesh size is increased.

Watershed application

The model is now applied to a real watershed on which outlet discharge measurements are available. At this scale, the finest mesh size used is 1 m and the Lawrence's friction model and the mixed friction model give the same results. The goal is to observe the effect of the numerical resolution on the representation of the outlet discharges.

On this domain, three mesh sizes are considered, with a space step of 1 m, 2 m and 5 m. The 1 m mesh has almost 2 millions elements and the computational time for a simulation is about 30 minutes. In the hydraulic network, the mesh is forced to follow the lowest point of the river. This leads to a greater representation of the channels geometry, in particular for the coarser mesh, and the water is flowing in the river direction. The model is calibrated to field data by comparing the measured and simulated outlet discharges on two events. The events are selected because they are the most erosive of the 2012 year, in terms of sediment volume exported at the outlet. The first event is recorded the 29th of May succeeding six rainy events from the 21st to the 27th of May. The maximal value of the rain intensity is 84 mm/h and the peak discharge is 3 m 3 /s. Concerning the other event, it is a summer storm of the 25th August with an initial dry state, the last recorded event being the 25th July. The maximal intensity of the rain is 156 mm/h and the peak discharge is 6.6 m 3 /s.

To simulate these events, we set the properties of the soil constant in time and space, modifying only the initial soil moisture (θ i in the infiltration model). The calibration has been performed on the finest mesh and gives the following soil parameters:

• first layer thickness: Z c = 80 mm,

• first layer porosity: θ 1 = 0.35,

• first layer conductivity: K c = 30 mm/h,

• second layer porosity: θ 2 = 0.25,

• second layer conductivity: K s = 1 mm/h,

• capillarity head: h f = 50 mm.

Inspired by the measurements of [Mathys, 2006], the parameters are chosen to have a porous surface layer and a more structured base layer. The value chosen for the base layer conductivity is clearly underestimated to balance the effects of the exfiltration, which is not represented by the model. the parameters representing the soil properties have been set to be as close as possible to the measurements. Then, the calibration process was carried out in two steps: a first step in which the initial water content was set to represent the peak intensity for each event, and a second step in which the representative size of the roughness was set according to the peak arrival time for the two events. The initial soil moisture is set to θ i =0.22 for the spring event and θ i =0.03 for the summer storm. Figure 13 presents the measured and simulated outlet discharges for the spring event and the summer event. The simulated results are in good agreement with the observations. The outlet hydrographs are well represented for the two events, especially for the 1 m and 2 m meshes. For the 5 m mesh, the flow is slowed and a dephasing is observed. The table 4 shows the RMSE on the outlet hydrographs normalized by each maximal value for the two events. The mixed friction model is used with a roughness size fixed to 5 cm, as the RMSE on the outlet discharge is minimized for this roughness size. Figure 14 shows the RMSE in function of the roughness size parameter of the Lawrence's friction model, and both events are minimized for the same value. The soil properties are considered as uniform in space, as well as the friction coefficient. The representation of the vertical structure of the soil is sufficient to have a satisfactory reproduction of the hydrographs. By isolating the areas where the error on the resolution is low according to the criteria given in 3.1, it appears that in the hydraulic network one has a good representation of the velocities by the numerical scheme. However, figure 15 shows that the dynamic is a little different for the two events. For the spring event where the soil is already wet and therefore the infiltration capacity is lower, the well resolved zone is better distributed upstream of the network. Concerning the summer event, the infiltration capacity is higher, but also is the discharge. The quantity of water concerned by the well resolved zone is more important downstream of the basin, but is less high in the network. The accuracy of the hydrographs seems to be in relation with the presence of large well resolved zones close to the outlet, where the discharge is measured. Indeed, the RMSE of the simulated outlet hydrograph of the summer event is much lower than the one for the spring event.

Conclusion

This paper focuses on the relevance of the Shallow Water Equations to model rain-induced runoff modelling and particularly on the evaluation of the local velocities. Four experimental test cases were used to analyze and assess the numerical model and to better understand its behavior. The first experiment is a flume with a constant slope on which drops a constant rain used to evaluate the error due to the numerical scheme. The second test case is a sinusoidal mould with furrows representing the agricultural terrain. Then the conclusions from these two laboratory experiments were used to interpret the simulated results obtained on the Thies test case (10 meters long real sandy plot) in comparison with measurements. Finally, the model was applied on a real catchment (Laval catchment, in the southern French Alps) where only outlet discharge measurements were available.

The flume test case demonstrates that the error of the numerical resolution with the [Chen and Noelle, 2017]'s scheme can be quantified. This test case shows explicitly that this numerical error can lead to very important deviations on the simulated velocities for well reproduced mass fluxes. Thus, two criteria were defined to reproduce conditions in which the error is acceptable. The first one is verified if the water depth is two times higher than the bottom elevation difference between two adjacent cells. The second criterion is verified if the slope and the runoff water level are high enough for the slope source term gh∂ x z to be greater than 2 • 10 4 m 2 /s 2 .

The [Lawrence, 1997]'s friction model was used for its ability to reproduce a wide range of inundation ratios. Comparisons with velocities of the second test case showed that Lawrence's friction model is well adapted but presents a limitation to reproduce the highest velocities. This limitation was corrected by adding a laminar friction law to represent the velocities in a well inundated regime combined with low Reynolds numbers. The difference between Lawrence's model and the mixed friction model became less visible when the mesh size was larger, which was the case for larger scale applications (Thies plot and Laval watershed test cases).

On the Thies plot test case, the numerical error criterion allowed to discriminate well reproduced velocities from the other ones. When the numerical criteria are respected, the mixed friction model gave very good results compared to the measured velocities. The error coming from the friction representation could not be compared with the error coming from the numerical resolution, but the numerical criterion found was able to frame the good representation of the velocities.

On the Laval watershed, two rainy events with two different seasons and initial conditions were selected. The outlet hydrographs were well reproduced on both events with the mixed friction model. The numerical criterion showed that the velocities are well reproduced mainly in the hydraulic network. Depending on the initial soil moisture and the rain intensity, the well resolved zone (the zone respecting the numerical criterion) spread more or less in the gullies upstream.

This work proposes a new methodology to evaluate the quality of the representation of the velocity field in a watershed for rain-induced runoff. It highlights well known problems in solving the Shallow Water Equations in the presence of steep slopes and low water depths [START_REF] Delestre | A limitation of the hydrostatic reconstruction technique for shallow water equations[END_REF], [START_REF] Kim | Coupled modeling of hydrologic and hydrodynamic processes including overland channel flow[END_REF]). This provides an objective criterion for evaluating the discretization of the slope source term in a real case. This was possible thanks to a friction law identification work adapted to each studied test case allowing to neglect the impact of the physical representation of this one compared to the numerical error committed on the discretization of the slope source term.

Moreover, this methodology is built using a first-order numerical scheme implemented in open source computational code. This, combined with the choice of the right friction law for each case, makes this methodology transferable to other numerical schemes, and can be used as a basis for an objective choice of the numerical method to be used for velocity modelling for rain-induced runoff cases.

Figure 5 :

 5 Figure 5: Bottom elevation of the Laval watershed.

Figure 6 :Figure 7 :

 67 Figure 6: Normalized Root Mean Square Error made in the simulations compared to the theoretical source term error.

  Figure 8: Temporal evolution of the outlet discharge (left) and velocities along the channel at the steady state (right) for the three chosen test cases, simulations vs measurements.

Figure 9 :

 9 Figure 9: Water depth (m) on the domain, experimental measure (left) versus simulation results (right) for the three configurations.

Figure 10 :

 10 Figure 10: Measured velocities compared to simulated velocities, comparison between the Lawrence's friction model and the mixed friction model for the sinusoidal mould test case.

Figure 11 :

 11 Figure 11: Rain and outlet discharge of the event in the plot for Thies case.

Figure 12 :

 12 Figure 12: Measured velocity compared to simulated velocities, comparison between the Lawrence's friction model (left) and the mixed friction model (right) depending on the value of β for Thies case.

Figure 13 :

 13 Figure 13: Measured outlet discharge compared to simulated results for the watershed case.

Figure 14 :

 14 Figure14: Normalized RMSE on the outlet discharge for the two events considered for the watershed case.

Figure 15 :

 15 Figure 15: Location of areas with minimal resolution error for the watershed case.

Table 1 :

 1 summarizes the simulated cases. The infiltration is not taken into account. Configurations for the flume test case simulations.

	Slope (%) Rain intensity (mm/h) Data available
	2	25	X
	5	25	X
	25	25	X
	2	50	X
	5	50
	25	50	X

Table 2 :

 2 Mesh convergence analysis on the velocities.

	Configuration ∆x (mm) RMSE (m/s)
	1	2	0.0658396
		5	0.0718715
		10	0.0781000
	2	2	0.0584047
		5	0.0626524
		10	0.071708
	3	2	0.0384667
		5	0.0410869
		10	0.0471743

Table 3 :

 3 Mesh convergence analysis on the velocities for Thies case.

	∆x (mm)	RMSE (m/s)	RMSE (m/s)
		Lawrence's model mixed model
	25	0.031507	0.030972
	50	0.078562	0.078575
	100	0.089981	0.089981

Table 4 :

 4 RMSE on the outlet hydrographs for each event and each mesh size for the watershed case.

	Event	∆x (m) RMSE (-)
	29/05/2012	1	0.041959
		2	0.052898
		5	0.239660
	28/08/2012	1	0.025281
		2	0.028930
		5	0.155558

Acknowledgements

This work is a part of EDF R&D thesis within the CIFRE agreement number 2015/1102. The authors would like to thank O. Planchon, G. Nord, C. Le Bouteiller and E. Audusse for the constructive discussions.


L Daudeville

P Ladeveze

A damage mechanics tool for laminate delamination

A simplified method based upon damage mechanics for the delamination analysis of carbon-resin composites is presented. In the neighbourhood of a laminate structure quasi-straight edge, damage is taken concentrated on the interface between layers. The finite element code EDA, acting as a post processor of an elastic laminate shell computation, allows the onset and propagation forecast of delamination. First numerical simulations of delamina tion are given and compared with experimental results from literature.

INTRODUCTION

laminate shell computation. They are associated to criteria based upon the average of normal The degradation modes of carbon-resin lamistresses along a characteristic distance from the nates (e.g. T300-914) can be split into two classes. free edge.5 Progressive degradations are not yet On the one hand, intralaminar damages: transtaken into account. In addition, delamination does verse matrix cracking, fibre-matrix debonding, not always occur where stresses are maxima and fibre ruptures. On the other hand, delamination, the intrinsic feature of the distance from the free i.e. degradation of the link between layers. Lamiedge can be discussed when geometry and stacknate structures being thin (plates or shells), the ing sequence vary. Linear fracture mechanics, i.e. stress state far away from the edge zones (free computation of an energy release rate G and its edges, zones of force application) is a plane stress comparison with a critical value Ge, is generally state, whereas the stress state can be three-dimenused for crack propagation study. 6 -8 Nevertheless, sional in the edge zones. Thus, delamination, fracture mechanics cannot be used for crack onset a phenomenon of layer debonding due to possible study. stresses that are normal to the shell surface,

Progressive degradation modelling allows one occurs principally in the edge zones of laminate to predict accurately the structure behaviour until structures.

rupture: that is damage mechanics of compo-Numerical rupture simulation of a laminate sites.9-11 Numerical methods based upon damage structure, taking into account all the progressive mechanics proposed by Ladeveze, 1 1 acting as degradation phenomena, leads to a time depenpost-processors of an elastic laminate shell comdent three-dimensional non-linear problem.

putation have been developed for the forecast of Nowadays, finite element calcualtions under a (i) rupture far away from the edge, (ii) rupture in classical laminate theory and combined with the vicinity of a circular hole, and (iii) delaminastress criteria, 1 are generally used for the rupture tion near a quasi-straight edge. forecast far away from the edge zones. Such an Method (iii) is presented. The general ideas of approach can be discussed because such criteria that computational tool were first presented in cannot well describe the different damage pheno-Ref. 12. We present now the development and the mena.

first results of this simplified method that allows For the delamination analysis, more especially the forecast of delamination onset and its propafor the onset forecast of a free edge crack, compugation on a short distance. Damage phenomena tation methods of elastic edge effects2-4 are are taken concentrated on the interfaces between generally used as post-processors of an elastic the layers. The edge is quasi-straight, thus the for the identification and the checking of the modelling.

DELAMINATION ANALYSIS BY DAMAGE MECHANICS

.1 Laminate modelling

The laminate is modelled as a stacking (i) of homogeneous layers, and (ii) of interfaces con necting layers (see Fig. 1 ).

All damage phenomena are taken concentrated on the interface. The principal features of the interface are as listed below. Interface.

• Its behaviour is assumed orthotropicthe influence of the adjacent layers on the inter face damage behaviour is taken into account through the fibre direction bisectors that are supposed to be the directions of orthotropy.

The constitutive law links the normal stresses to the jump of displacement between layers, its expression in elasticity is with • Damage is unilateral in the 3 directionfor the normal to the shell direction (mode I)

I [U] = U+ -U-= [u1]N1 + [u2]N2 + [u3]N3
there is no damage under compression.

• Relative variations of stiffness are the damage indicatorsthree internal variables d; are associated to the three stiffness moduli k7.

• The interface behaviour is assumed to be elastic with damagethat is at first, an inelastic behavour can also be introduced. 1 7

The thermodynamical potential is the strain energy:

= _!_ (<-033)� + (033)� ED 2 k� k�(l -d3) ( • )+ denotes the positive part.
The variables Ydi' that are similar to the energy release rate introduced in fracture mechanics, are conjugated to d;: Y d; =a E ofad;

1 < 033)� Y,13 = 2 k�( 1 -d3) 2 mode 3 =mode I 2 1 0 31 Yd i = 2 k7(1-d i ) 2 mode 1 2 1 032 }d 2 = 1 k�(l -d2 ) 2 mode2
Note that modes 1, 2 and 3 are here relative to the axes of orthotropy, contrary to the fracture mechanics classical appellation that refers to crack opening directions.

{d3 =w(Y) if d3< 1 d1 =y1w(Y) if d1<1 and d3< 1 d 2 = y2 w(Y) if d 2 < 1 and d3 < 1 with d3 = 1 otherwise d 1 = 1 otherwise d 2 = 1 otherwise
where Y0 is threshold energy, Ye is critical energy, n is characteristic parameters of the damage evolution law of the interface.

Link with fracture mechanics

Classical tests of fracture mechanics18 allow us to obtain the three inter-laminar fracture toughnesses G1c, Gue, Gm0 relative to the modes I, II and Ill. G1c, Gue, Gme are different. [START_REF] Sela | lnterlaminar fracture toughness and toughening of laminated composite materials: a review[END_REF] Under the assumption that Ge is J. constant parameter (there is no R-curve effect), Ge is a characteristic interface parameter. Ge can be interpreted as the necessary work per unit surface for the interface debonding. Then, it can be dependent (for instance under pure mode I= mode 3) on the characteristic parameters of the interface modelling.

J [ u3 , [

Glc = 0 033d[ U3] with [ u3e] = [ u3] Jd3=1: the interface is locally cracked
G1c is, therefore, the area under the curve 033 -[u3] (see Fig. 2). In pure mode I one finds

n G!c= Yo + --(¥;,-Yo) n+l
For the identification under pure modes 2 or 3, one assumes that the crack opening direction coincides with a bisector of adjacent fibre directions. In that case one finds (mode II= mode 1 and mode III = mode 2 ):

Gae =Yo (� -1) + (l� n+l) ( 1 -(�+!) (¥;,-Yo) ) Y; Y; Y; with la= II, III i= 1, 2
The above relations allow a first identification of damage evolution law parameters.

Note that delamination modelling through an interface with damage includes the notion of process zone given by some authors, because damage varies continually in the vicinity of the crack tip. The process zone length depends on the structure (the damaged zone at the crack tip is different for DCB beams with layers at 0° and 90°).

3 NUMERICAL STRATEGIES The laminate shell, of thickness 2h, occupies the domain Q = Q 1 U Q2 (see Fig. 3 ).

When the thickness 2h is weak compared with the transverse dimensions, the three-dimensional displacement-stress solution (U, o) of the elastic problem set on Q can be written as follows:20 r u =Up+ UE with r (Up, Op): solution of a shell problem lo= op+ oE VuE, oE): solution localised near the edge aQ2

Notably because of the layer constitutive law discontinuity, the following boundary condition is not respected on the edge 0Q2:

\lzE[-h, h] op(z)n=F(z)
where F and n denote the force density and the normal. Hence (UE, oE) is added to the valid solution far away from the edge aQ2(Up, Op) such that (U, o) satisfies the equations and boundary conditions of the three-dimensional problem. (UE, oE) is localised near aQ2, thus the edge effect solution (UE, oE) can be computed in post-processor of a shell computation set on Q2• If, in addition, the radius of curvature of the edge aQ2 is large compared with the thickness, the variations of (UE, oE) along the tangent direction ( y) are negligible, the problem becomes set into a band$ of which the length is in order of the thickness 2 h.

The problem (.9) to solve becomes: find (UE, oE)

The force R applied on the edge ao2 n $ of the structure is the residue of external forces and of the shell solution:

R=F-opn (i) means that (UE, oE) is localised and (iii) is the constitutive law.

Solving the boundary layer problem in the non-linear case

In the proposed modelling of the laminate, only the interfaces can be damaged due to edge effect normal stresses. The interfaces taking no part in the computation of (Up, ap ), it is possible to apply the principle of superimposed stress for the computation of (UE, aE).

The problem (.9J) to solve becomes a non-linear time dependent problem set into the band.%'.

We compared the damage states near the free edge obtained by two non-linear computations of a specimen submitted to tension along the x direction (see Fig. 4 ).

The first computation C 1 has been carried out by resolution of problem (.9J ). The second one C2 has been a reference computation because it allows us to have the exact three-dimensional solution by searching the displacement under the form proposed by Pipes and Pagano.22 {u(y, z) +t:x v(y, z) V(x, y, z):

( ) w y, z C1 and C2 results are rigorously equal. That confirms it is possible to apply the principle of superimposed stress in that particular case where non-linearities are taken concentrated on interfaces.

A Riks-like method

Due to damage on interfaces, an instability point may appear. That critical point cannot be passed with a Newton method that pilots the computation in terms of 'force' (see Fig. 5). A Riks-like method 1 5 • 1 6 allows to control the computation and pass such a limit point.

Current step:

By the Newton method: OAi = 0 { Xuo;= Ai+1F-f(U;) (.97 ) A i+i =A i+ oAi constraint: g( OAi) = 0 By the Riks method the load factor is released, it is then necessary to impose a constraint: g( OA i) = 0 Riks 1 5 and Crisfield 1 6 propose a constraint g( cUi) = 0 that concerns the norm of the global vector oi. Such a global constraint may lead to a non-convergence when damage localises. To ensure a good convergence we propose to use a local constraint that considers only the more significant degrees of freedom in the increase of damage. Note a is the closer node of the Gauss point where the increase of damage has been the more important at the initial step (/3 is the node located on the same interface but on the adjacent layer). n (1, 2 or 3) is the mode of principal damage (see Fig. 6).

The local constraint consists in imposing a constant value for the jump of displacement between a and /3 along the n direction. 3.3 A method for a large number of layers The operator Jr;, (see (9' )) can be the tangent structural stiffness matrix %1 or the secant one.�, (i.e. computed through an elastic with damage behaviour). Solving (9') can be carried out with the direct Crout method. In the case of a large number of layers, that method may be expensive because of .�1 matrix size. We propose then to use the conjugate gradient method with use of a conditioner .

. %a=X1:

In the iterative Riks algorithm, the linear systems to solve can be written under the form .Jt;, 0 0 andJtim are respectively the contributions of layers (Jt;,00 is a constant matrix) and of interfaces Ptini is a non-constant matrix because the interfaces can damage). The conjugate gradient method consists in a series of iterating resolutions of The matrix..X:: 0 0 (conditioner) being diagonal by units, the resolution is parallel on each layer. 4 The interface can be considered as a rich resin zone of weak thickness e1 compared with the layer thickness ec. By this analogy and by considering that e1/ec � 1, the interface constitutive law can be expressed as

NUMERICAL SIMULATIONS OF DELAMINATION

k o 2 --2G 23 with e 1 ( 1
)
where G13, G23, £3: elastic moduli of the rich resin zone, they can be chosen equal to the homo genised layer moduli. The identification of the other parameters can be carried out by the use of classical fracture mechanics tests21 results for the mode I (specimen DCB see Fig. 7), the mode II (specimen ENF) and the mode III (specimen SCB). The use of tough ness values issued from literature 1 9 is difficult because those values can present, for the same material (for instance, T300-5208), an important divergence. For the following simulations, we identified the parameters by comparison with experimental results. been carried out with no initial crack, whereas the experimental specimen had an initial crack of 28•5 mm. The authors 23 found Gic = 450 J /m 2 and the following elastic moduli (the elastic moduli relative to the normal direction 3 are supposed to be equal to those in the transverse direction 2 ).

n G 1c = 1';1+ --(� -Y o) n+l
In view of numerical results, one finds

(2)

• the numerical simulation predicts quite well the experimental results since the crack length (simulation) reaches the experimental initial crack length; • Y0, Ye, n values have an influence on the (F, U) curve (see Fig. 8), at the condition these values satisfy eqn ( 2 );

• by the knowledge of the load F associated to the crack length a, it is possible to compute the critical energy release rate G 1c (by the area method 1 8).

One finds the experimental value of G1c and one verifies that the critical energy release rate is independent of the crack length (see Fig. 9). This result has been verified experimentally 23 and theoretically established14 for an initial crack which is large compared with the thickness 2h. The important point is the numerical simulation pre dicts this result even if the crack length is short.

Delamination near the free edge of a specimen under tension or compression

We compared our simulations of delarnination near the free edge of a specimen under tension or compression (see Fig. 4) with first experimental results5• 2 4 for mode I delarnination (on the mid plane interface) of a T300-5208 material (see Table 1 ). Then, we also compared our results with experimental ones 2 5 for mixed mode of delamina tion of a T300-1034C (see Ta ble 2). The G1c-a curve. We have chosen the interface stiffnesses by eqn

5 1---+------+-------. . -t------t-•--t-• -+--t-
(1) and e1 = ec/5. First, we have identified under mode I (T300-5208) the parameters �. �1 and n with the (±45,0,90), laminate. Secondly, we have iden tified the parameters y 1 and y 2 under mixed mode (T 300-1034C) such that the results are quite good. So, we assumed that both T300-5208 and T 300-1034C damage behaviours were close.

The calculated strain ( Ecail in the tables corre sponds to the maximum of strain of the strain crack length curve (see Fig. 10).

S CONCLUSION

The analysis of delamination near an edge of a laminate structure by damage mechanics has been presented. The basic assumptions of this method the edge is quasi-straight, damage is taken concentrated on interfaces between layerslead to a simplified method for the delamination onset putation. Note there is no mesh dependency because the treated problem is plane and damage is concentrated on interfaces. The first compari sons of numerical simulations with experimental results are quite good.

  initial three-dimensional problem becomes a two dimensional one set into a band perpendicular to the edge. The laminate is modelled as a stacking of homogeneous elastic layers connected by elastic with damage interfaces. The interface is a two-dimensional entity which ensures stress and displacement transfers between layers. The degradation effect is taken into account through the relative variations of the interface elastic moduli. An interface with damage has also been used recently by Schellekens 13 for delamination simulation of a laminate under tension. The inter face modelling of this author is principally based upon numerical considerations. The inter face modelling we propose has already been pre sented. In the previous presentation, 1 4 for the particular example of a double cantilever beam (DCB), a link has been established between damage mechanics of the interface and fracture mechanics giving a mechanical interpretation to that modelling. Note that the crack length was considered large compared with the beam thick ness. In the present paper, the numerical simula tion of the DCB test is performed with no initial crack. The finite element method is used. Due to damage, the structure may present a critical state. A method 1 5 • 1 6 in which crack opening displace ment is controlled, allows such a limit point to be passed. Numerical simulations of delamination, com puted with the post-processor edge damage analy sis (EDA), are given in the framework of (i) a DCB beam: study of onset and propagation of delamination, and (ii) a specimen under tension: study of delamination onset near the free edge. The proposed numerical simulations allow a com parison with experimental results from literature

Fig

  Fig. 1.Interface.

  2.1.1 A first choice for the damage evolution lawThe evolution of the three damage variables is linked because the same micro-cracks participate to damage. Damage evolution is assumed to be governed by Y =sup lr,,/ Yd3 + y1 Yd, + y 2 Yd2) y1, y 2 coupling factors Damage evolution law is

3. 1

 1 The boundary layer problem 3.1.1 Problem formulation in elasticity

Fig. 2 .

 2 Fig. 2. Interface elastic with damage behaviour (mode I traction).

Fig. 3 .

 3 Fig. 3. Boundary layer.

  ( o ;)� -( o ;)� = 0 for i � 1

Fig. 4 .

 4 Fig. 4. Specimen under tension.

Fig. 5 .

 5 Fig. 5. Riks-like method.

4. 1

 1 Identification of model parametersThe characteristic parameters of the interface modelling are k�, kg, k�, }(J> ¥,, , n, Y 1 > y2•

4. 2

 2 TheDCBWe have simulated the DCB test with the EDA program for a comparison with the experimental results2 3 (the simulation of the DCB test is not a boundary layer problem). The studied structure is constituted of 48 layers of thickness ec = 0• 13 7 5 mm and of T 300 fibres with an MlO resin. The numerical simulation has

Fig. 6 .

 6 Fig. 6. Concerned nodes by the local constraint.

  2h

Fig. 7 .

 7 Fig. 7. DCB test.

  Yci. Ye, n values satisfy the following relationship:

Fig. 8 .

 8 Fig. 8. F-U curve.

  Fig. 9.

  ( Oz , ± 1 5 2 ), 15, -15 -0•51 (04, ± 154), 15, -15 -0•35 (04, ± 154), 15, -15 0•57 "Data taken from Ref. 25. propose the same elastic moduli for both materials T300-5208 and T300-1034C.

E 11 =

 11 138 000 MPaE22 = E33 = 9700 MPa G 1 2 = G13 = G23 = 5500 MPa V12 = v 1 3 = 0•3, v23 = 0•6

  Fig. 10. r-a curve for a (O, ±45, 90), laminate: delamination on the mid-plane interface.

Table 1 .

 1 Mode I delamination"

	Laminate	Cexp
	(±45,0,90t	0•53
	( ± 45 z , O z , 90 z ), ( ± 45 3, 0 3 , 90 3ls (0, ± 45, 90),	0-45 0•36 0•66
	(45, 0, -45, 90),	0•54
	( ± 30, 90), ( ± 30 2. 90 2) ,	0•39 0•36

"Data taken from Refs 5 and 23.

Table 2 .

 2 Mixed mode delamination"

	Laminate	Interface	Eexp