ENAC Capucine Amielh

ENAC Alexandre Chabory

ENAC Christophe Macabiau

Laurent Azoulai

Airbus Biographies

GNSS Multipath Error Model for Aircraft Surface Navigation Based on Canonical Scenarios for Class F Airports

In 2016, she worked as a research engineer in the ElectroMagnetic and Antenna (EMA

INTRODUCTION

For surface operations, aircraft ground position and ground velocity can be estimated using a large variety of sensors, such as Global Navigation Satellite System, Inertial Navigation System, Wheel Speed Sensors and others. The sensor bringing global position information with a few meters of error is the GNSS, but measurements from this sensor are affected by signal blockage and multipath. The use of GNSS to navigate with guaranteed integrity on the airport surface relies on a valid and efficient GNSS surface multipath error model for airport surface navigation but such kind of multipath error model has not yet been standardized.

In [START_REF] Chen | Development of a GPS Deterministic Multipath Simulator for an Efficient Computation of the Positioning Errors[END_REF], a GPS multipath simulator for airport navigation has been proposed. By means of Geometrical Optics and Physical Optics, this electromagnetic model permits to simulate up to second order interactions between obstacles. The complex airport environment is rendered with uncertainties in the simulated scene such as the building positions and materials. This deterministic approach is therefore enhanced by MonteCarlo draws on these model parameters as presented in [START_REF] Chen | Hybrid deterministicstatistical GPS multipath simulator for airport navigation[END_REF]. A multipath ranging error model based on the outputs of the Fortran Multipath Range Error simulator (FMRE) has been proposed in [START_REF] Montloin | GNSS multipath error model for airport surface operations[END_REF]. Overbouding laws are obtained in [START_REF] Montloin | GNSS multipath error model for airport surface operations[END_REF] to get the mean and the variance of bounding distributions on the pseudorange errors due to multipath. Moreover, the multipath simulator has been recently enhanced by an antenna+aircraft model specific to ground airport navigation. This electromagnetic model accounts for aircraft structure multipath as well as antenna gain, phase and group delay variations for all direction of arrival in the presence of the close fuselage. Details on this radiating antenna+aircraft model can be found in [START_REF] Amielh | Importance of the Antenna Model to Assess GNSS Multipath in Airport Environments[END_REF] and [START_REF] Amielh | Including the Aircraft and the Antenna in a Wide Band Aeronautical LMS Channel Model[END_REF].

Because of the growth of the number of passengers carried by air transport and the dimensions of new aircraft models, a standardized multipath pseudorange error model seems more and more required to ensure a safe ground navigation and to maintain high airport capacity. To assess the performance of the FMRE multipath pseudorange error model in these conditions, this study is focused on reference code F airports in ICAO classification [START_REF]International Civil Aviation Volume I, Aerodrome Design and Operations[END_REF]. The following study is made for GPS L1 C/A signals but is applicable to L5, E1C and E5a.

The first basic objective of this paper is to propose a validation of the multipath pseudorange error simulator through the numerical validation of the electromagnetic prediction of multipath parameters.

Then, the main objective is to propose a multipath pseudorange error model based on the one proposed in [START_REF] Montloin | GNSS multipath error model for airport surface operations[END_REF] adapted to canonical taxi operation scenarios. An example of establishment of this multipath pseudorange error model is proposed for a typical scenario of taxi on a stand taxilane.

Pseudorange error multipath simulator for aircraft ground surface navigation

In this first section, the principle and the validation of the GNSS multipath pseudorange error simulator is described.

Pseudorange error multipath simulator principle

For a given aircraft trajectory and satellite position, the computation of the multipath error is based on a deterministic elec tromagnetic method applied to an airport scene (Fig. 1). The elements in the scene, e.g. airport buildings or surrounding aircraft, are described by means of a polygonal mesh. The simulator emulates first and secondorder interactions as follows.

For firstorder interactions, each facet yields a scattered field modeled by means of Physical Optics (PO). Currents on the facet induce an electromagnetic (EM) field by means of StrattonChu integrals. This way, each illuminated facet is seen as the origin of one multipath. The characteristics of the multipath are computed taking into account the antenna model which is detailed in [START_REF] Amielh | Importance of the Antenna Model to Assess GNSS Multipath in Airport Environments[END_REF] and [START_REF] Amielh | Including the Aircraft and the Antenna in a Wide Band Aeronautical LMS Channel Model[END_REF].

For secondorder interactions, Geometrical Optics (GO) is used to determine the illuminated facets and the associated incident fields of the first interaction. GO is a ray launching method which accounts for the reflection of rays by surfaces e.g. facets here. When an incident ray impacts a surface, it yields reflected rays, whose direction derives from SnellDescartes laws, and whose amplitude and divergence factor are modified. The divergence factor characterizes the evolution of the field amplitude due to the ray tube divergence with propagation. GO is limited as for instance it predicts a discontinuity of the field at shadow boundaries. The second scattering is obtained with PO similarly to firstorder interactions. The reflection over the ground is accounted via the image theorem in both cases.

After passing through the receiver antenna, each scattered field yields one multipath characterized by an amplitude, a phase, a delay and a Doppler shift. Amplitude and phase are directly related to the scattered field while the delay and Doppler shift are assessed geometrically from the satellite position and the aircraft trajectory. Finally, the pseudorange is determined from multipath parameters via a GNSS receiver model. The receiver model simulates correlator outputs in the presence of multipath.

Since the accuracy of the 3D environment model and aircraft trajectory are limited, the deterministic EM multipath model is enhanced by statistical draws associated with several environment uncertainties. This MonteCarlo analysis considers randomness for the building positions and orientations in the horizontal plane, their height and materials, as well as the material thicknesses, the ground material and the receiver altitude. At each point of a trajectory, this affects the mean and the variance of the pseudorange error.

Numerical validation of the EM multipath prediction of the multipath simulator

In this subsection, a numerical validation of the multipath pseudorange error simulator is proposed. This presented validation addresses the validation of the EM multipath model. The receiver simulator has been validated in various efforts not related in this paper.

The EM multipath model provides multipath parameters as inputs to the receiver simulator. To validate this multipath param eter prediction, one way is to compare the channel transfer function obtained with the multipath simulator to a reference transfer function derived from a different tool than FMRE. In this paragraph, the reference channel transfer function is obtained using FEKO, an EM software. Under FEKO, the Method of Moments is used. This exact method derives an analytical resolution of Maxwell's equations. The FEKO channel transfer function is computed through Eq. 1.

H F EKO (f ) = E RHCP (θ, ϕ, f ) E RHCP (θ = 0, ϕ = 0, f = f 0 ) (1) with
• θ and ϕ, the spherical coordinates,

• f , the frequency with f 0 the antenna central frequency equal to f L1 = 1575.42 MHz

• E RHCP , the righ handed cicrularly polarized electric field,

• and H F EKO , the normalized channel transfer function at the frequency f in the direction of arrival (θ, ϕ).

This reference solution is compared to the FMRE predicted transfer function of Eq. 2.

H F M RE (f ) = N ∑ n=0 a n e j(2πτn(f +f Dopp n )+ϕn) (2)
where • a n , the amplitude of the n th multipath,

• τ n , the delay of the n th multipath,

• ϕ n , the phase of the n th multipath,

• f Dopp n , the Doppler frequency of the n th multipath,

• and H F M RE , the FMRE channel transfer function at the frequency f in the direction of arrival (θ, ϕ).

N.B.:

To simplify the notation, the transfer functions are written H(f ) even if they depend on the direction of arrival and the Doppler frequency is omitted since we consider a static configuration.

The comparison of the transfer functions is executed in the context of a patch antenna mounted on an infinite ground plane and facing a metallic wall of 5 m by 5 m, 60 m away as presented is Fig. 2. To compare the transfer functions, similarly to [START_REF] Chen | Development of a GPS Deterministic Multipath Simulator for an Efficient Computation of the Positioning Errors[END_REF], the mean square (RMS) difference ϵ RM S is computed after applying the baseband power spectrum density of the GPS L1 C/A signal (cf. Eq. 3).

ϵ RM S = ∫ B/2 -B/2 |S GP S L1CA .(H F M RE (f ) -H F EKO (f ))| 2 df ∫ B/2 -B/2 |S GP S L1CA | 2 df (3)
with the normalized baseband power spectrum density of the GPS L1 C/A signal:

S GP S L1CA = sinc 2 (T c (f -f L1 ))
where T c the chip spacing is equal to 1 fc = 1 1.023e 6 = 0.9775 µs. If the RMS difference is lower than -20 dB then it corresponds to a power difference of 1% and therefore the two predictions fit.

The antenna model under FMRE does not render the frequency dependentbehavior of the antenna on the complete L1 C/A bandwidth (refer to [START_REF] Amielh | Importance of the Antenna Model to Assess GNSS Multipath in Airport Environments[END_REF]). In fact, the real antenna gain, phase and group delay are varying with the frequency but under the multipath simulator we only consider their value at the antenna central frequency i.e. f L1 for complexity reasons. Therefore, this antenna difference is removed from the FEKO transfer function as follows. Due to the theorem of reciprocity and superposition, we can write

H F EKO (f ) = 1 + |H W all+Ant Bias F EKO (f )| -|H Ant Bias F EKO (f )| (4)
H W all+Ant Bias F EKO (f )
is the transfer function obtained with FEKO for the patch antenna mounted on an infinite ground plane and a 5 m by 5 m metallic wall at 60 m. H Ant Bias F EKO (f ) is the transfer function of the same antenna radiating in free space. In the following, H F EKO (f ) refers to the transfer function resulting from Eq. 4 and is free of frequencydependent antenna biases. These transfer functions are plotted on Fig. 3. As it can be seen on Fig. 3, the major part of the antenna bias has been successfully removed. A decrease of the amplitude of the transfer function is still visible along with the frequency. The FMRE and FEKO transfer functions present the same oscillations for a similar magnitude and fit quite well. These oscillations are characteristics to the presence of a dominant multipath.

The transfer functions for an azimuth from 0 to 360 o and an elevation from 0 to 90 o are shown in Fig. 4 and5. The period of the oscillations is comparable for FEKO and FMRE transfer functions. The magnitude is of the same order even if we observe a frequencydependent amplitude residue up to 1 dB for some directions of arrival. The difference between the two transfer functions according to Eq. 3 is presented Fig. 6. On Fig. 6, we can see that the criteria ϵ RM S < -20 dB is mainly respected except for some azimuth at high elevation. In fact, for a satellite close to the zenith, the multipath simulator is still able to predict some multipath while in reality (see Fig. 7), as predicted by FEKO, the major part of the scattered electric field is not redirected towards the antenna. These directions with a RMS difference larger than to -20 dB only represent 13 cases on 5429 meaning that the RMS criteria is respected for 99,76% of the directions of arrival and then that the EM prediction of multipath is correct.

The next section illustrates the use of the validated multipath pseudorange error simulator for measurements comparison.

Comparison to measurements performed at ENAC

A measurement campaign has been performed at ENAC on a rooftop to limit the influence of non predictable obstacles such as pedestrians, bikes or cars and facing a metallic wall to maximize the chance to get a dominant multipath. The antenna used is an Antcom G8 Arinc mounted a ground plane. This antenna is a dualband GPS L1 C/A and L2 patch. We used a Septentrio PolarRX3eTR receiver. The position of the receiver antenna is measured using a Velodyne Lidar. The Lidar image of the scene is shown Fig. 8 and its representation under the multipath simulator is visible Fig. 9. We measured the code and the phase on L1 C/A and L2 during more than 6 hours. The multipath pseudorange error is extracted by means of a Code Minus Carrier process presented in details in [START_REF] Amielh | Validation of Existing GNSS Multipath Model[END_REF]. The resulting CMC are plotted Fig. 10. The existence of multipath is directly linked to the geometry of the scene. Therefore, for a sufficiently long period of obser vation, the satellite movement is significant enough to see a variation of the multipath parameters. This phenomenon implies that multipath variations are a low frequency component. Higher frequencies of the PSD mainly represent noise and diffuse multipath.

5 satellites seem to have significant low frequency PSD peaks. In the following, it has been chosen to present results for PRN 29 which shows the highest PSD amplitude at low frequency. The azimuth, elevation and CMC of PRN 29 can be seen Fig. 13 and the signaltonoise ratio Fig. 14. The mean and the standard deviation of the predicted multipath pseudorange error at each epoch is computed and plotted Fig. 16. We have to keep in mind that the CMC values are not noise free contrary to our multipath simulator predictions. To remove this noise term, one could model the signaltonoise ratio as a function of the elevation and link this to a noise versus elevation model. This would allow us to compare multipath pseudorange error variances free of noise. But since the 6 hours of measurement campaign are not sufficient to establish this model on this day.

Even if the CMC still includes the noise, the comparison provides good results as seen Fig. 16 and 17. Two phases can be distinguished: a transient state where the satellite is not locked yet and a permanent state after the last tracking loss. The transient state occurs at low elevation. The loss of tracking can be due to the presence of a tree not modeled in the simulator scene. Trees are not modeled since this multipath simulator is designed for airport navigation only. The prediction oscillations confirm the presence of a dominant multipath. With Tab. 2, we can compare the mean and the standard deviation of the observed CMC and the predictions of the multipath error variations for both transitory and permanent state. A good agreement of the means and the standard deviations can be observed except for the standard deviation of the transitory state. In fact the peaks occurring right before the loss of tracking is not present in the simulations. This is typically a result of a not modeled obstacle such as a tree or a building detail. Nevertheless, for the permanent state, the means perfectly match and the standard deviations are quite similar taking into account that the measured standard deviation contains a noise term. The comparison between simulated multipath pseudorange errors and CMC values provides an interesting illustration of the use of the multipath simulator. Even if the comparison with CMC is limited due to the noise presence and the small number of MonteCarlo draws, we get a good overall accordance.

Quantity

In this section, a numerical validation of the multipath simulator has been proposed. To illustrate the results which can be made with the multipath simulator, a comparison to measurements performed at ENAC has been developed. These two conclusive steps allow us to go further with this multipath simulator and to propose a multipath ranging error model for airport navigation.

Multipath pseudorange error model for canonical taxi operation scenarios

The purpose of this paper is to propose a new GNSS multipath pseudorange error model adapted to taxi operations on a class F airport. In [START_REF] Montloin | GNSS multipath error model for airport surface operations[END_REF], a multipath ranging error model has been proposed but this model was not airport independent. The next subsections refer to the model presented in [START_REF] Montloin | GNSS multipath error model for airport surface operations[END_REF] and goes further to adapt it to any class F airport by identifying classical situations.

Previous multipath pseudorange error model for taxi operations from [2]

First of all, the model of [START_REF] Montloin | GNSS multipath error model for airport surface operations[END_REF] starts by identifying an impact zone. This impact zone is defined as a location on the airport surface where multipath has a significant contribution to the pseudorange error with respect to the other sources of ranging errors (troposphere, ionosphere, satellite clock and ephemeris inaccuracies and the receiver thermal noise) in steady state. Once this location is delimited, the multipath ranging error model distinguishes a static and a dynamic situation.

Static configuration

When the aircraft is not moving in the impact zone, the multipath pseudorange error ε M P is decomposed as follows:

ε M P = b A+G (Θ, Φ) + b obs (5)
where

• Θ and Φ, the satellite elevation and the aircraft azimuth,

• b A+G , the aircraft+ground component induced by the ground reflection, which is considered as a deterministic bias,

• b obs , the oscillating component induced by the 1 st and 2 nd order interactions with obstacles.

The oscillating component depends on the GNSS receiver antenna radiation pattern model, the GNSS receiver settings, the satellite position, the aircraft azimuth, the aircraft model itself, the 3D model of the airport environment, and the type of ground material. The Probability Density Function (PDF) of the oscillating component b obs is then overbounded by a centered Gaussian distribution.

Dynamic configuration

In the context of a dynamic situation, the multipath ranging error is estimated via Eq. 6.

ε M P (t) = b A+C (Θ, Φ) + u obs (t) (6)
with

• b A+C (Θ, Φ), similar to the static one,

• u obs (t), the oscillating component induced by the 1 st and 2 nd order interactions with obstacles, a zeromean timedependent error.

The oscillating component u obs has the same dependency as b obs but it depends on a trajectory location and orientation dependency as well as an aircraft speed over the trajectory dependency. In [START_REF] Montloin | GNSS multipath error model for airport surface operations[END_REF], it has been chosen to provide a stochastic model of the oscillating component in the impact zone. The true trajectory is modeled by its own statistical distribution drawn from Monte Carlo simulations. In fact, in each MonteCarlo simulation, each point of a trajectory is randomly chosen within an interval around the true position. This drawn position ensures the true trajectory independence. The Power Spectral Density (PSD) of the oscillating component is overbounded by a zeromean 1 st order GaussMarkov process M(τ dyn , σ dyn ). τ dyn and σ dyn are respectively the correlation time and the standard deviation. The normalized PSD of the 1 st order GaussMarkov process is given in Eq. 7.

S GM (f ) = 2σ 2 dyn τ dyn ( (2πf ) 2 + 1 τ 2 dyn ) (7)
The overbounding PSD is chosen by fixing τ dyn graphically from simulation results. Then the minimal value of σ dyn is picked so that the PSD of the GaussMarkov process overbounds the estimated PSD. This multipath error model has been only tested on some specific trajectories on Blagnac airport and therefore, other conditions may lead to different results. The choice of a GaussMarkov process has also been made to accommodate the use of a further integrity algorithm.

The following multipath pseudorange error model is independent of the considered airport and the aircraft orientation and trajectory since it is based on canonical scenarios that can be found in ICAO documentation.

Canonical taxi operation scenario identification

Even if an airport is difficult to model due to all the uncertainties of the environment, the airport navigation is quite codified. In the following section, canonical scenarios are identified according to the ICAO airport navigation rules.

In [START_REF]International Civil Aviation Volume I, Aerodrome Design and Operations[END_REF], an ICAO airport classification is proposed with respect to the dimensions of the biggest aircraft which can use the platform. This classification is reminded Tab. 3.

Code letter

Wing span

Outer main gear wheel span A Up to but not including 15 m Up to but not including 4.5 m B 15 m up to but not including 24 m 4.5 m up to but not including 6 m C 24 m up to but not including 36 m 6 m up to but not including 9 m D 36 m up to but not including 52 m 9 m up to but not including 14 m E 52 m up to but not including 65 m 9 m up to but not including 14 m F 65 m up to but not including 80 m 14 m up to but not including 16 m Table 3: ICAO aerodrome reference code [START_REF]International Civil Aviation Volume I, Aerodrome Design and Operations[END_REF] For example, according to these dimensions an A380 can only navigate on a class F airport. In this same document, ICAO provides a list of the minimal distances between taxiway center lines and between taxiway center line and objects. A minimal wingtip to object distance is also given. These distances vary with the taxi operation subphases and the aerodrome reference code. In the following the surface operation areas are in accordance with the ICAO definitions (Fig. 18). The subphases of taxi operations which are considered are all the phases from taxi on the rapid exit taxiway to taxi on the aircraft stand leadin lines.

Figure 18: Surface operation areas [START_REF]International Civil Aviation Volume I, Aerodrome Design and Operations[END_REF] The surface clearances for class F airport used in this multipath ranging error model are extracetd from [START_REF]Aerodrome Design Manual Part 2 Taxiways, Aprons and Holding Bays 4th edition[END_REF] The object can be another taxilane center line or an air craft parked ... Apron taxiway centerline to object

d AT xw2o

The object can be an aircraft parked, a taxilane center line, an apron taxiway center line or a taxiway center line ... Taxiway center line to taxi way center line

d T xw2T xw

Table 4: Surface clearance for airport navigation on class F airport [START_REF]International Civil Aviation Volume I, Aerodrome Design and Operations[END_REF] Based on these surface navigation clearances and a study of airport environment, scenarios of ground navigation have been identified. They can be classified by obstacles in the scene:

• scenarios including airport buildings and surrounding parked aircraft in Tab. 5,

• scenarios including only surrounding aircraft in Tab. 6,

• and a scenario including only small obstacles such as cars, trucks, gateway, blast fence ... in Tab. 7.

The purpose of the small obstacles scenario is to determine which obstacles have a multipath pseudorange error contribution sufficient to be included into other scenarios. To ensure the worst case scenario, the surrounding parked aircraft are always A380 and the navigating aircraft is an A320. Considering an A320 navigating along with A380s may not be the most relevant scenario but it ensures the worst case. The chosen speeds are consistent with usual speeds in zerovisibility conditions.

The different conditions of distances and speed are summarize in the Comments of Tab. 5, 6 and 7. For example, when an aircraft passes by another aircraft, depending if the aircraft at stop is parked on the apron or at stop on another taxiway, the distances between both aircraft are different. All the possible cases are considered in the following tables. It means that for each presented configuration, several scenarios are considered. Table 7: Small obstacle scenario

The proposed multipath ranging error model aims at providing the mean and the variance of the multipath pseudorange error for each canonical scenario. The global error model for a given aircraft trajectory from leaving the gate to entering the runway or from taking the rapid exit taxiway to parking at the gate will be a juxtaposition of these canonical scenarios.

The following section presents the multipath ranging error model applied for each of these scenarios.

Establishment of the multipath pseudorange error model

For each scenario, hundreds to thousands MonteCarlo draws are computed with the multipath simulator. The parameters of the MonteCarlo draws are :

• the receiver antenna altitude rendering the uncertainty on the position of the antenna,

• the building position and orientation to account for the environment geometric uncertainties,

• the building material permittivities and thicknesses,

• the building heights,

• the building facade materials distribution to be independent of the type of building, depicted Fig. 19, • the satellite azimuth.

Figure 19: Building random material distribution for a ratio concrete/metal/glass of 50/10/40% Therefore, the obtained GNSS multipath pseudorange error model reflects the uncertainty on the previously listed parameters. The model is also made to be independent of the airport considered. The different aircraft orientations and trajectories are taken into account as different scenarios. Finally, the list of input parameters is as follows: the satellite elevation, the navigating aircraft type and the subphase of taxi operation via the assessed scenario. The navigating aircraft type is taken into account in the antenna+aircraft model where the airframe multipath is assessed and in the scene as the altitude of the receiver antenna depends on the dimensions of the aircraft.

Similarly to the dynamic configuration of the multipath ranging error in [START_REF] Montloin | GNSS multipath error model for airport surface operations[END_REF], we propose to overbound the Power Spectral Density (PSD) of the multipath pseudorange error using a first order GaussMarkov process.

Contrary to what was made in [START_REF] Montloin | GNSS multipath error model for airport surface operations[END_REF], the airframe multipath is gathered with the other obstacle multipath before being aggregated into delay bins. The deterministic bias b A+C cannot be retrieved in the final simulated multipath pseudorange errors. Therefore, it has been chosen to overbound the complete multipath pseudorange error ϵ M P by a GaussMarkov process.

ϵ M P (t) = U obs (t) (8)
with U obs a first order GaussMarkov process rendering the influence of the aircraft structure, the ground and the obstacles. The parameters of this GaussMarkov process are determined as follows. First, to ensure the independence of the ranging error overbound to the satellite azimuth, for one elevation, we aggregate the PSDs over all the azimuths. Then, for a range of correlation time, the GaussMarkov density larger in each point than the maximum of the simulated PSDs is computed. Among these PSDs, the GaussMarkov density with the minimal standard deviation is kept for each considered correlation time. The search is performed through a simple dichotomy. Then the optimal overbounding GaussMarkov density is chosen as the one minimizing the integral difference

∫ f S σmin (τ, f )df - ∫ f S + simu (f )df (9)
with S σmin the PSD with a minimal standard deviation for the given correlation time τ and S + simu the maximum of the simulated PSDs at each frequency. The way to proceed is presented in the next section with the example of taxi on a stand taxilane.

Example of multipath pseudorange prediction for taxi on stand taxilane

To illustrate the establishment of this multipath pseudorange error model, in this section we look at a scenario of taxi on a stand taxi lane between two airport buildings with 6 A380 parked.

Figure 1 :

 1 Figure 1: Multipath simulator principle

Figure 2 :

 2 Figure 2: Scene configuration for transfer function comparison

Figure 3 :

 3 Figure 3: Magnitude transfer function comparison at elevation 5 o

Figure 4 :Figure 5 :

 45 Figure 4: Magnitude of the FEKO transfer function with antenna bias removed for 0 to 360 o of azimuth and 0 to 90 o of elevation

Figure 6 :

 6 Figure 6: RMS difference of the two transfer functions for 0 to 360 o of azimuth and 0 to 90 o of elevation

Figure 7 :

 7 Figure 7: Magnitude transfer function comparison at elevation 89 o

Figure 8 :

 8 Figure 8: Lidar image of the measurement conditions Figure 9: Representation of the measurement environ ment under the multipath simulator

Figure 10 :

 10 Figure 10: Estimated multipath and noise pseudorange errors

Figure 11 :

 11 Figure 11: Power Spectral Density function of the CMC observables Figure 12: Zoom in low frequencies of CMC's PSD

Figure 13 :

 13 Figure 13: Azimuth, elevation and CMC over time for PRN 29

Figure 15 :

 15 Figure 15: Illumination of the scene by PRN29

Figure 16 :

 16 Figure 16: CMC versus mean ± the standard deviation of the predictions for PRN 29 Figure 17: Zoom in the early permanent state of Fig. 16

  stand If the aircraft is on a taxilane, D1 = d T xl2o , D2 = d M iniStndSp , D3 = d ClearReq and the speed is from 5 kt to stop. But if the aircraft is on an apron taxiway, D1 = d AT xw2o and the speed is from 5 kt to stop. Taxi on taxilane or apron taxiway af ter push out be tween two A380 from stand If the aircraft is on a taxilane, D1 = d T xl2o , D2 = d M iniStndSp , D3 = d ClearReq and the speed is from stop to 5 kt. But if the aircraft is on an apron taxiway, D1 = d AT xw2o and the speed is from stop to 5 kt. Taxi on taxilane or apron taxiway par allel to stand with 3 A380 parked If the aircraft is on a taxilane, D1 = d T xl2o , D2 = d M iniStndSp , D3 = d ClearReq and the speed is inferior to 5 kt. If the aircraft is on an apron taxi way, D1 = d AT xw2o and the speed is between 5 and 20 kt. And if the aircraft is on a taxiway, D1 = 2d AT xw2o and the speed is between 5 and 20 kt. Taxi on taxi lane/apron taxiway between two par allel stands with 3 A380 parked at each If the aircraft is on a taxilane, D1 = d T xl2o , D2 = d M iniStndSp , D3 = d ClearReq and the speed is inferior to 5 kt. If the aircraft is on an apron taxi way, D1 = d AT xw2o and the speed is between 5 and 20 kt. And if the aircraft is on a taxiway, D1 = 2d AT xw2o and the speed is between 5 and 20 kt.Table 5: Scenarios including buildings Scenario name Scenario illustration Comments Passing by a small obstacle D1 = d AT xw2o , 360 o trajectory at 5 kt max.

  Fig.20shows a representation of the scene under the multipath simulator. The navigating A320 is moving along a straight trajectory of 80 m at 3 kt. Due to the geometric symmetry of the scene, the satellite describes only half the space i.e. the satellite azimuth takes a range between 0 and 180 o . The azimuth is given relatively to the xaxis. The receiver RF filter has a ideal infinite bandwidth and the signal considered is the GPS L1 C/A. A total of 500 MonteCarlo draws are performed.

Figure 20 :

 20 Figure 20: Taxi on aircraft stand taxilane Fig. 21 presents the multipath pseudorange errors obtained at elevation 10 o and azimuth 60 o for some MonteCarlo draws. The associated PSD are plotted in Fig. 22.

Figure 21 :

 21 Figure 21: Multipath pseudorange errors obtained with the different MonteCarlo draws for an elevation of 10 o and an azimuth of 60 o

Figure 22 :

 22 Figure 22: Power spectral density of the multipath pseudorange errors obtained with the MonteCarlo draws for an elevation of 10 o and all the considered azimuths

Figure 23 :

 23 Figure 23: Optimal overbounding GaussMarkov deter mination Figure 24: Total simulated PSDs with max and over bounding PSDs

Table 1 :

 1 1. Number of trajectory per aircraft type

	Category	Parameter	Variation
	Receiver	Receiver altitude	10 cm
		Horizontal center translation	1 m
	Building geometry	Vertical height homothetie	1%
		Horizontal rotation	3 %
		Wall thickness	50 cm
	Material	Window thickness	3 cm
		Air layer thickness	1 cm %

Table 2 :

 2 Moments comparison

	Tansitory state Permanent state

  in Tab. 4.

	Distance name	Value (m)	Notation	Comments
	Clearance requirement	7.5	d ClearReq	Can be increased if required by safety but we focus on
				worst case
	Minimum stand spacing	87.5	d M iniStdSP Sum of the maximal wingspan authorized and the clear
				ance requirement
	Stand taxilane center line to	50.5	d T xl2o
	object

ACKNOWLEDGMENTS

The authors thank GSA and Airbus for founding this research.

The standard deviations are of the meter level and as expected decrease with the elevation. It can be noted that at elevation 45 o the standard deviation is higher than at 30 o . This can be explained by the fact that we considered only 6 different azimuths. Therefore the presented standard deviations don't represent the complete satellite geometry.

Besides, the correlation times are very small and decrease with the elevation. The evolution of the correlation time seems to present an anomaly. Indeed, for a satellite at high elevation, a small satellite position change is expected to have a small impact on the simulated multipath and therefore the multipath pseudorange errors should be correlated over a longer period. But the presented scenario corresponds to a 16 s trajectory at 50 Hz which prevents from getting access to low frequency PSD components. An extension of this study to lower frequency and more azimuths should enable to determine more accurately the correlation time and the standard deviation.

CONCLUSION

To conlude, the Fortran Multipath Range Error simulator principle has been presented. The outputs of the simulator have been illustrated through comparison to measurements performed at ENAC. Then, a numerical validation of the electromagnetic multipath prediction has been proposed based on transfer functions comparison to the FEKO reference. Then after reminding the work of [START_REF] Montloin | GNSS multipath error model for airport surface operations[END_REF], we adapted a multipath pseudorange error model for aircraft navigation on reference code F airport. This model is built for canonical scenarios. Its input parameters are the navigating aircraft type, the satellite elevation and the taxi operation subphase considered as among a list of canonical scenarios. In this model, the power spectral density of the multipath pseudo range error is overbounded by a first order GaussMarkov PSD which seems to fit quite well in the example of taxi on a stand taxilane.

To go further, a study of this scenario with a smaller azimuth grid step will be performed to get access to a better estimation of the correlation time and the standard deviation. Then this multipath ranging error model will be considered as the input of an hybridization algorithm to estimate a navigation solution for class F airport navigation.


J C González-Aguirre

M J Castro

Morales De Luna

A robust model for rapidly varying flows over movable bottom with suspended and bedload transport: modelling and numerical approach

Keywords: Sediment transport, Shallow water Exner models, bedload, suspended load, finite volume solvers, well-balanced schemes

We propose a coupled model for suspended and bedload sediment transport in the shallow water framework. The model is deduced under hydrostatic pressure assumptions and will not assume any Bossinesq hypothesis. The numerical resolution is carried out in a segregated way. First the underlying system of conservation laws is solved by using a first order path-conservative Riemann solver. Then, the source terms corresponding with the erosion and depositions rates are approximated in a semi-implicit way. The final scheme preserves the positivity of the density. Several numerical experiments were carried out in order to validate the model and the numerical scheme. The results obtained are in good agreement with the experimental data.

Introduction

In geophysics and earth sciences the term fluvial is used to refer to the processes associated with the rivers or streams and erosion or deposits and the morphology created by them. The sediment is transported by the river current ˚Dpto. Análisis Metemático, Universidad de Málaga, Spain

Email address: jcga@uma.es (González-Aguirre, J. C. )

as suspended load (finer fractions carried by the flow) and bed load (coarse fractions which move close to the bottom rolling, jumping and sliding), see [START_REF] Subhasish | Fluvial Hydrodynamics, Hydrodynamic and Sediment Transport Phenomena, GeoPlanet: Earth and Planetary Sciences[END_REF].

Knowledge of sediment transport can be applied in civil engineering, for example, to plan the extended life of a dam forming a reservoir. The sediment transport can be applied to solve environmental engineering problems, for example, it is important in providing habitat for fish in river or estuaries, sustaining a hygienic stream ecosystem. The sediment carried by a river is deposed downstream reducing the capacity of the river at that area, which can results

in frequent floods at outskirts of the river.

The study of sediment transport focuses on understanding the relationship that exists between the movement of water and the movement of sedimentary materials. A first approach to model bedload transport is to couple the shallow water equations with the so-called Exner equation (See [START_REF] Exner | Über die wechselwirkung zwischen wasser und geschiebe in flüsen[END_REF]). This equation depends on the empirical definition of the solid transport flux for the bedload transport. Several formulations have been proposed, see for instance [START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF][START_REF] Einstein | The bed load function for sediment transportation in open channel flows[END_REF][START_REF] Grass | Sediment transport by waves and currents[END_REF][START_REF] Van Rijn | Sediment transport, part I: Bed load transport[END_REF] among others. This approach has been extensively used to describe bedload transport, see [START_REF] Caleffi | High-order balanced cweno scheme for movable bed shallow water equations[END_REF][START_REF] Castro-Díaz | Sediment transport models in shallow water equations and numerical approach by high order finite volume methods[END_REF][START_REF] Liu | Two-dimensional scour simulations based on coupled model of shallow water equations and sediment transport on unstructured meshes[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] Murillo | An exner-based coupled model for twodimensional transient flow over erodible bed[END_REF][START_REF] Soares-Frazão | Hllc scheme with novel wave-speed estimators appropriate for two-dimensional shallow-water flow on erodible bed[END_REF][START_REF] Juez | A 2D weakly-coupled and efficient numerical model for transient shallow flow and movable bed[END_REF], among many others.

Suspended load accounts for sedimentary particles which have been eroded

from the bottom and remain floating on the current for some time until they are sedimented afterwards. In [START_REF] Pantin | Interaction between velocity and effective density in turbidity flow: Phase-plane analysis, with criteria for autosuspension[END_REF] was proposed a model with full considerations of the mass and momentum transfers due to the sediment exchange between the bottom and the water. This model includes the mass and momentum conservation equation for the water sediment mixture and the mass conservation of the sediment, but the bed deformation was not taken into account. In [START_REF] Cao | Computational dam-break hydraulics over erodible sediment bed[END_REF][START_REF] Hu | Fully coupled mathematical modeling of turbidity currents over erodible bed[END_REF] the mass conservation for the bed material was considered in order to incorporate the effects of the bed deformation on the fluid flux. Similar models have been used in [START_REF] Simpson | Coupled model of surface water flow, sediment transport and morphological evolution[END_REF][START_REF] Wu | One-dimensional modeling of dam-break flow over movable beds[END_REF][START_REF] Li | Fully coupled approach to modeling shallow water flow, sediment transport, and bed evolution in rivers[END_REF][START_REF] Benkhaldoun | A flux-limiter method for dam-break flows over erodible sediment beds[END_REF][START_REF] Canelas | Two-dimensional depth-averaged modelling of dam-break flows over mobile beds[END_REF][START_REF] Liu | A robust coupled 2-D model for rapidly varying flows over erodible bed using central-upwind method with wetting and drying[END_REF][START_REF] Liu | Wellbalanced central-upwind scheme for a fully coupled shallow water system modeling flows over erodible bed[END_REF]. All of these models are based on the assumption of hydrostatic pressure. A more recent work [START_REF] Cantero-Chinchilla | Nonhydrostatic dam break flows. II: one-dimensional depth-averaged modeling for movable bed flows[END_REF] introduces a non-hydrostatic model for sediment transport.

One of the assumptions that is usually done in sediment transport models is the Boussinesq hypothesis. This means that spatial variations of the density of the water sediment mixture are ignored in the convective part and are only taken into account into the pressure term. This assumption is valid when the variations of density are small. This can be inaccurate in many situations. One objective of this work is to deduce a general mathematical model for sediment transport without such assumption. This model will be proposed under the shallow water framework and aims at accurately describing both bedload and suspended sediment transport. Up to our knowledge, sediment transport models based on the shallow water framework usually use Boussinesq hypothesis which makes the numerical computations easier. Nevertheless this may result in inaccurate results in situations where we are far from this assumption.

A second objective of this paper is to described an efficient and robust numerical approximation of the model. This will be done by using a a two-step algorithm, following similar ideas as in [START_REF] Nord | Psem 2d: A physically based model of erosion processes at the plot scale[END_REF][START_REF] Kubatko | An unstructured grid morphodynamic model with a discontinuous galerkin method for bed evolution[END_REF][START_REF] Morales De Luna | On a shallow water model for the simulation of turbidity currents[END_REF][START_REF] Castro-Díaz | A HLLC scheme for nonconservative hyperbolic problems. application to turbidity currents with sediment transport[END_REF]. During the first step, the hyperbolic part of the system will be solved using a path-conservative wellbalanced numerical scheme (See [START_REF] Parés | Numerical methods for nonconservative hyperbolic systems: a theoretical framework[END_REF][START_REF] Xing | On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations[END_REF]). Moreover, both bedload and suspended load will be solved without decoupling the system. As shown in [START_REF] Cordier | Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help[END_REF], a splitting technique may result into instabilities. During the second step, erosiondeposition terms are solved using a novel semi-implicit approach. This will allow to guarantee the positivity preserving of density and water-thickness.

This work is structured as follows: First, the mathematical model is deduced.

In Section 3, a general compact formulation will be presented as well as a study of the hyperbolicity of the model and the characterization of some relevant stationary solutions. The numerical scheme to be used will be introduced in Section 4. In Section 5, the extension to the two-dimensional case will be briefly shown. Finally, in Section 6, several numerical experiments are carry out to validate the model and the numerical approach used.

Deduction of mathematical model

We consider a water flow over movable bottom bpx, tq that contains suspended particles of diameter d and density ρ s . The movement of water is characterized by the time evolution of both the water depth hpx, tq and the velocity upx, tq " pupx, tq, vpx, tq, wpx, tqq where x P R 3 . We shall denote by ηpx, tq " hpx, tq `bpx, tq the free surface of the fluid and by ρ 0 the density of water. The evolution of the bottom, assumed with porosity p, is the result of two processes: erosion-deposition of the bottom and bedload transport. The first one corresponds to a transfer of particles between the bottom and the suspended load in the fluid. This process gives a mixture of water and sediment with density ρpx, tq and volumetric sediment concentration cpx, tq. The erosion and deposition fluxes will be denoted by F e px, tq and F d px, tq respectively.

The second process, called bedload transport, takes into account sedimentary particles which move along the bottom by rolling, sliding and jumping. This transport is measured by the bedload discharge q b px, tq. A one dimensional sketch of this phenomenon is shown in Figure 1. The mathematical model to be used is based on the Navier-Stokes equations:

B t ρ `divpρuq " 0. (2.1) divpρuq " 0, (2.2)
B t pρuq `divpρu b uq " ´∇p `∇ ¨T `ρg, (2.3)
where ppx, tq denotes the instantaneous pressure. The density of the mixture is given by ρpx, tq " ρ 0 `cpx, tqpρ s ´ρ0 q, (2.4) the vector g " p0, 0, ´gq t is the gravity and

T " µp∇u `p∇uq t q (2.5)
is the stress tensor, whith µ the viscosity coefficient. We shall consider the kinematic boundary conditions

B t η `u| z"η B x η ´w| z"η " 0, (2.6)
B t b `u| z"b B x b ´w| z"b " ´r φ b , (2.7)
where r φ b " p1 `pB x bq 2 q 1{2 ϕ b . Here ϕ b is the normal erosion-deposition rate.

For the sake of simplicity, we shall neglect the variations on the horizontal direction y. Nevertheless everything said here is easily extended otherwise.

Vertical integration

In this section we shall consider that the horizontal effects are predominant with respect the vertical ones, and therefore, we suppose that the vertical variations of the velocity are small. Under this assumptions we proceed to do a vertical integration of the equations. We shall assume that the fluid is incom-

B t pρuq dz `ˆη b B x pρu 2 q dz `ˆη b B x p dz " ρu| z"b w| z"b ´ρu| z"η w| z"η `µpr τ | z"η ´r τ | z"b q. (2.14)
Taking into account kinematic boundary condition, (2.14) becomes

B t ˆη b ρu dz `Bx ˆη b ρu 2 dz `Bx ˆη b p dz " ρu| z"b r φ b `µpr τ | z"η ´r τ | z"b q.
(2.15) Definition 2.1. For any function f px, z, tq, we define its mean on the vertical direction as

f " 1 h ˆη b f dz. (2.16)
The fluctuating part p f of f relative to the mean is defined as

p f " f ´f . (2.17)
From definition (2.1) we can see that, given two functions f 1 , f 2 we get

hf 1 f 2 `ˆη b p f 1 p f 2 dz " ˆη b f 1 f 2 dz. ( 2
B t ˆη b pρuq dz `Bx ˆη b pρu 2 q dz `ghρB x η `g h 2 2 B x ρ `g ˆη b ˆBx ˆη z p ρdz 1 ˙dz. " µpr τ | z"η ´r τ | z"b q `ρu| z"b r φ b . (2.23)
Using definition 2.1 and equation (2.17 (2.24)

We shall assume that the fluctuating parts are negligible and we consider that the tensor stress can be expressed in terms of u, so that, we write

µpr τ | z"η ´r τ | z"b q " τ puq. (2.25)
We shall denote

φ b " c r φ b . (2.26)
It is usual to express the velocity at surface, u| z"η and bottom u| z"b , in terms of mean velocity u. Entropy considerations leads us to assume u| z"η " u and u| z"b " 1 2 u, see [START_REF] Parker | Self-accelerating turbidity currents[END_REF][START_REF] Bouchut | On new erosion models of savage-hutter type for avalanches[END_REF]. Therefore, equations (2.19) and (2.24) take the 90 form B t phρq `Bx phρ uq " ρ 0 r φ b `φb pρ s ´ρ0 q, (2.27)

B t phρ uq `Bx phρ u 2 q `Bx ˆg h 2 2 ρ ˙" ´ghρB x b `τ puq `ρ u 2 φ b .
(2.28)

Henceforth we neglect the symbol ´. Following [START_REF] Hu | Fully coupled mathematical modeling of turbidity currents over erodible bed[END_REF], in order to take into account the bed deformation, the density of sediment ρ s in the source term is usually written in terms of the density of saturated bottom ρ b given by ρ b " ρ 0 p `ρs p1 ´pq, (2.29) therefore (2.27) becomes

B t ph ρq `Bx phρ uq " φ b ρ b 1 ´p `ρ0 ˆr φ b ´φb 1 ´p ˙. (2.30)
The dynamic of bedload sediment transport is modeled using the Exner equation, see [START_REF] Exner | Zur physik der dünen[END_REF][START_REF] Exner | Über die wechselwirkung zwischen wasser und geschiebe in flüsen[END_REF]. This equation states that the time variation of the sediment layer over a certain volume is due to of the net variation of the solid transport through the boundaries of the volume. Exner equation reads

B t b `1 1 ´p B x q b " ´1 1 ´p φ b , (2.31)
where q b " q b ph, h uq denotes the solid transport flux.

Let us consider the simple case where the sediment is confined in a closed domain, so that there is not exchange through boundaries. In particular, assume for simplicity the case where u " 0. Then, one should get that the total sediment mass in the interior of the domain is preserved, that is

B t phρ `ρb bq " 0. ( 2
B t phρuq `Bx phρu 2 q `Bx ˆg h 2 2 ρ ˙" ´ghρB x b `τ puq `ρ u 2 φ b , (2.37)
B t b `1 1 ´p B x q b " ´1 1 ´p φ b , (2.38)

Erosion and deposition effects

To close the system the sediment exchange between the bottom and the water (φ b ) is determined from the rates of deposition (F d ) and erosion (F e ), namely φ b " F e ´Fd .

(2.39)

Deposition rate

The deposition rate is described as the product of the settling velocity of sediment (drop velocity of the particle under water at rest) v s , and the fractional concentration of suspension near by the bed c b ,

F d " v s c b . (2.40)
In order to compute c b we shall use the formula proposed in [START_REF] Bradford | Hydrodynamics of turbid underflows. i: Formulation and numerical analysis[END_REF], which states that

c b c " 0.4 ˆd D sg ˙1.64 `1.64, (2.41)
where D sg denotes the geometric mean size of the suspended sediment mixture and d is the sediment diameter.

The settling velocity of sediment may be obtained experimentally or using the formulation proposed by Zhang and Xie [START_REF] Zhang | Sedimentation research in China: Systematic selection[END_REF], which has been used in several works like [START_REF] Benkhaldoun | A flux-limiter method for dam-break flows over erodible sediment beds[END_REF][START_REF] Liu | Wellbalanced central-upwind scheme for a fully coupled shallow water system modeling flows over erodible bed[END_REF]

v s " d ˆ13.95 ν d ˙2 `1.09 s g d ´13.95 ν d , (2.42)
where ν denotes the kinematic viscosity of the water, s " ρ s ρ 0 ´1 is the relative density and s g is the specific gravity of the submerged sediment.

Erosion rate

The erosion rate from the bottom is computed using the following expression F e " v s P E s .

(2.43)

The sediment entrainment coefficient E s , is computed using the formulation developed by Garcia and Parker [START_REF] Garcia | Experiments on the entrainment of sediment into suspension by a dense bottom current[END_REF],

E s " 1.3 ˆ10 ´7Z 5 1 `4.3 ˆ10 ´7Z 5 , (2.44)
where

Z " α 1 ? c D |u| v s R α2 j , (2.45)
and

R p " ? s g d d ν (2.46)
is the particle Reynold number. Several values have been proposed for the parameters pα 1 , α 2 q as it can be seen in [START_REF] Parker | Self-accelerating turbidity currents[END_REF][START_REF] Bradford | Hydrodynamics of turbid underflows. i: Formulation and numerical analysis[END_REF][START_REF] Kubo | Laboratory experiments and numerical simulation of sediment-wave formation by turbidity currents[END_REF][START_REF] Kubo | Experimental and numerical study of topographic effects on deposition from two-dimensional, particle-driven density currents[END_REF]. In this work we shall employ the values (1,0.6) if R p ą 2.36 and (0.586,1.23) when R p ď 2.36. The constant P represents the volume fraction of the sediment in the bottom. In [START_REF] Goñi | Determination of the volumetric solids fraction of saturated polydisperse ore tailing sediments[END_REF] is it said that if P represents the volumetric solids fraction of the sediment, then 1 ´P denotes the porosity of the control volume.

Friction term

Following [START_REF] Hu | Fully coupled mathematical modeling of turbidity currents over erodible bed[END_REF], the friction term is computed as

τ puq " ρu 2 ˚p1 `rw q, (2.47)
where u ˚is the bed shear velocity and it is computed as [START_REF] Cao | Computational dam-break hydraulics over erodible sediment bed[END_REF] by

u 2 ˚" gµ 2 f u 2 h 1{3 , (2.48)
where µ f is the dimensionless Manning coefficient. The term r w denotes the 120 ratio of upper-interface resistance to bed resistance.

Bed load discharge

The computation of the bed load discharge q b is based on empirical laws.

Some of the most used ones were proposed Meyer-Peter and Müller [START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF] , Grass [START_REF] Grass | Sediment transport by waves and currents[END_REF], Van Rijn [START_REF] Van Rijn | Sediment transport, part I: Bed load transport[END_REF], among others. In this work we use the formulation proposed by Meyer-Peter and Müller, namely

q b Q " 8 sgnpuq pθ ´θc q 3{2 `, (2.49)
where

Q " d a g s d
is the characteristic discharge, θ denotes the Shields parameter (ratio between agitating and stabilizing forces on a sediment grain in the bottom)

θ " u 2 s g d (2.50)
and θ c is the critical Shields stress for incipient motion.

General formulation of the model

The system of partial differential equations (2.35)-(2.38) can be written in vectorial form as

B t W `Bx FpWq `BpWqB x W " GpWq `Sτ pWq, (3.1)
where the vector of conservative variables takes the form W " ph, hρ, hρ u, bq t .

The physical flux F, the non-conservative product matrix B and the source terms, G (erosion and deposition) and S τ (friction term), have the followings form FpWq "

¨hu hρu hρu 2 q b 1 ´p ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , BpWq " ¨0 0 0 0 0 0 0 0 g 2 hρ g 2 h 0 ghρ 0 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , (3.2)
GpWq "

¨φb 1 ´p φ b ρ b 1 ´p ρ u 2 φ b ´φb 1 ´p ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , S τ " ¨0 0 ρu 2 ˚p1 `rw q 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' . ( 3.3)
System (3.1) may also be written in quasi-linear form as

B t W `ApWqB x W " GpWq `Sτ pWq, (3.4)
where ApWq " BFpWq BW `BpWq.

The problem to be solved can be formulated as follows. Given the initial condition W 0 , we need to compute the state W t , t ą 0 by solving the system (3.4) on r0, Ls, where L is the length of channel, subject to some boundary conditions.

Hyperbolic character of the system

Given that the bed load discharge flux does not depend on the depth of bottom, we have B b q b " 0. Let us be consider the matrix

A " ¨u ´u ρ 1 ρ 0 0 0 1 0 g 2 hρ g 2 h ´u2 2u ghρ α β δ 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , (3.5)
where

α " B h q b , β " B hρ q b , δ " B hρu q b .
Following [START_REF] Cordier | Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help[END_REF], bedload transport fluxes may be written as a function of the Shields parameter, that is,

q b " sgnpuqΦpθq,
for some function Φ. In particular, if the Meyer-Peter&Müller formula is used, we have Φpθq " 8Qpθ ´θc q 3{2 `. Nevertheless everything said here may be generalized for any other solid transport flux formula.

Therefore, we get

α " sgnpuqΦ 1 pθqB h θ, β " sgnpuqΦ 1 pθqB hρ θ, δ " sgnpuqΦ 1 pθqB hρu θ. (3.6)
Remark that from (2.50) we have

B h θ " ´1 6 ρu B hρu θ, B hρ θ " ´u B hρu θ. (3.7) Therefore α " ´1 6 ρuδ, β " ´uδ. (3.8)
Some straightforward calculations show that the characteristic polynomial of A is ppλq " ´pu ´λqpλ 3 `a2 λ 2 `a1 λ `a0 q, (3.9) with a 2 " ´2u,

a 1 " ´ˆghρ ˆδ `1 2ρ ˙`´g 2 h ´u2 ¯˙, a 0 " ´ghρ ˆβ `α ρ ˙.
Thus, u is an eigenvalue of the matrix A. In order to prove the hyperbolicity of the system (3.4) we need to ensure that the cubic polynomial on the right hand side of (3.9) has three real roots. Remark that this polynomial may be written as

P A pλq " f pλq ´dpλq, (3.10)
where

f pλq " λ `pu ´λq 2 ´gh ˘, (3.11)
dpλq " ghρ ˆλ ´7 6 u ˙δ.

This means that we have three real roots whenever the line dpλq intersects f pλq in three points. In other words, we need to ensure that the root of dpλq, that More explicitly, the lines d ˘pλq take the form

d ˘pλq " f pλ ˘q ´ghρδpλ ´λ˘q , (3.13)
where

λ ˘" 2u `au 2 `3ghp1 ´ρδq 3 . (3.14)
The roots of the straight lines d ˘pλq are given by

α ˘" f pλ ˘q ghρδ `λ˘. (3.15)
System (3.4) will be hyperbolic whenever the next relation is fulfilled: Remark that a sufficient condition for the system to be hyperbolic is that |u| ă 6 ? gh, which is in general true for physical situations. Indeed, one can check that in this case 150 α ´ă u ´agh ă 7 6 u ă u `agh ă α `.

α ´ă 7 6 u ă α `. (3.16)
(3.17)

In order to compute an approximation for characteristics velocities, we shall employ the following procedure, which is depicted in Figure 3.

• Let us consider

S ´" u ´agh, S mid " 0, and S `" u `agh,

the roots of (3.11).

• We compute the tangent straight line l ˘pλq to f pλq at point pS ˘, 0q

l ˘pλq " f 1 pλq pλ ´S˘q . (3.19)
• We solve the equation l ˘pλq " ´dpλq, whose solution is • We compute the tangent straight line l mid to f pλq at point (0, 0) and we solve the equation l mid pλq " ´dpλq, then we get

λ ˚" f 1 pS ˘qS ˘´7 6 ghρuδ f 1 pS ˘q ´ghρδ . ( 3
λ m " ´7 6 ghρuδ u 2 ´ghp1 `ρδq . (3.21)
Let us consider λ 1 ă λ 2 ă λ 3 the roots of the polynomial (3.10). We shall use the approximations given by λ ˚ă λ m ă λ ˚.

155

The approximation λ ˚, is an upper bound for λ 3 , while λ ˚is a lower bound of λ 1 .

Stationary solutions

In this section we briefly consider some relevant stationary solutions of the system (3.4), when the source term related to de sediment erosion and depo-160 sition processes and friction term are neglected. In that case, the momentum conservation equations writes

u pB t phρq `Bx phρ uqq `hρ pB t u `uB x uq `g 2 h 2 B x ρ " ´ghρB x pηq. (3.22)
Taking into account the mass conservation equation in (3.4) and that stationary solutions are independent of time, we get In particular we are interested in those corresponding to u " 0. In such case, two different families of stationary solutions appear, those corresponding to lake at rest given by u " 0; ρ " cte; η " cte.

B x u 2 2 `gh 2ρ B x ρ `gB x η " 0. ( 3
(3.26)

In the case that the density is not a constant function, we get It is expected that the suspended sediment under water at rest must be deposed at the bottom. At the same time, due to this deposit, the thickness of the bottom increases, water thickness decreases. This will modify the density of the water-sediment mixture as follows:

gh 2 2 B x ρ `ghρB x η " 0. ( 3
• The density of mixture is a decreasing function through the time if c ă 180 1 ´p.

• The density of mixture is an increasing function through the time if c ą 1 ´p.

• The density of mixture is a constant function if c " 1 ´p.

Numerical scheme

In this section we describe the first order numerical finite volume solver that we used to approximate the solutions of system (3.4).

The spatial discretization is done by means of a finite volume mesh composed by cells C i " rx i´1{2 , x i`1{2 s, i " 1, ¨¨¨, N . For the sake of simplicity we shall assume that these cells have a constant size ∆x and denote by x i`1{2 " i ∆x the intercells and x i " `i `1 2 ˘∆x the center of the cells, see Figure 4. Let ∆t be the time step and t n " n∆t.

0 L x1 x 1{2 xi´1 xi xi`1 x i´1 2 x i`1 2
Ci ' E We use a two-step algorithm to compute the solution of the system (3.4).

Incorporating erosion and depositions effects in the second step as it is done in [START_REF] Nord | Psem 2d: A physically based model of erosion processes at the plot scale[END_REF][START_REF] Kubatko | An unstructured grid morphodynamic model with a discontinuous galerkin method for bed evolution[END_REF][START_REF] Morales De Luna | On a shallow water model for the simulation of turbidity currents[END_REF][START_REF] Castro-Díaz | A HLLC scheme for nonconservative hyperbolic problems. application to turbidity currents with sediment transport[END_REF].

On the first step, we solve the system

B t W `ApWqB x W " 0, (4.1)
using a finite volume well-balanced method and on the second step we solve the ordinary differential equation induced by the erosion deposition and friction source terms

B t W " GpWq `Sτ pWq. (4.2)
In order to solve (4.1) we shall use the path-conservative framework. Remark that the system (4.1) contains non conservative products ApWq W x , which in general do not make sense within the framework of the theory of distributions.

The theory described in [START_REF] Maso | Definition and weak stability of nonconservative products[END_REF] allows to give a rigorous definition of weak solutions using a family of Lipschithz continuous paths. This theory has been used to developed the so-called path conservative schemes, see [START_REF] Parés | Numerical methods for nonconservative hyperbolic systems: a theoretical framework[END_REF].

Let W n i be the approximation of the cell averages of the exact solution at the cell C i at time t n

W n i - 1 ∆x ˆxi`1{2 x i´1{2 Wpx, t n qdx, (4.3)
One can write a numerical scheme in the form

W n`1 i " W n`1 i ´∆t ∆x ´Dì ´1 2 `Dí `1 2 ¯, (4.4)
with

D ȋ`1 2 " D ˘pW i , W i`1 q " 1 2
´FpW n i`1 q ´FpW n i q `Bi`1 2 pW n i`1 ´Wn i q 1 2 Q i`1 2 `Wn i`1 ´Wn i ˘, (4.5)
where the matrix Q i`1 2 represents a viscosity matrix that approximates

|A i`1 2 |,
the absolute value of the Roe linearization for A at the intercell. In particular, we will use a PVM method as introduced in [START_REF] Castro Díaz | A class of computationally fast first order finite volume solvers: Pvm methods[END_REF].

Q i`1 2 " P i`1{2 r ´Ai`1 2 ¯(4.6)
where

P i`1{2 r
pxq is a polynomial of degree r, which approximates the absolute 205 value function.

In this work we employ a parabolic approximation in order to compute the viscosity matrix Q i`1 2 . Further details can be found in [START_REF] Degond | Polynomial upwind schemes for hyperbolic systems[END_REF][START_REF] Fernández-Nieto | On an intermediate field capturing riemann solver based on a parabolic viscosity matrix for the two-layer shallow water system[END_REF],

P i`1{2 2 ´Ai`1 2 ¯" α 0 Id `α1 A i`1 2 `α2 A 2 i`1 2 . ( 4.7)
where α 1 , α 2 , α 3 are defined in terms of the eigenvalues λ 1 " λ min , λ 2 " λ mid , λ 3 " λ max and they are computed by solving the linear system of equations

¨1 λ 1 λ 2 1 1 λ 3 λ 2 3 1 λ 2 λ 2 2 ‹ ‹ ‹ ' ¨α0 α 1 α 2 ‹ ‹ ‹ ' " ¨|λ 1 | |λ 3 | |λ 2 | ‹ ‹ ‹ ' . ( 4.8)
Let us considerer

γ 1 " ´`|λ 3 |pλ 1 λ 2 2 ´λ2 1 λ 2 q `|λ 2 |pλ 3 λ 2 1 ´λ2 3 λ 1 q `|λ 1 |pλ 2 λ 2 3 ´λ2 2 λ 3 q ˘, γ 2 " |λ 2 |pλ 2 1 ´λ2 3 q `|λ 1 |pλ 2 3 ´λ2 2 q `|λ 3 |pλ 2 2 ´λ2 1 q, γ 2 " ´p|λ 2 |pλ 1 ´λ3 q `|λ 1 |pλ 3 ´λ2 q `|λ 3 |pλ 2 ´λ1 qq, γ 4 " pλ 1 ´λ2 qpλ 3 pλ 3 ´λ1 ´λ2 q `λ1 λ 2 q,
then the solution of the linear system (4.8) is given by

α 0 " γ 1 γ 4 , α 1 " γ 2 γ 4 , α 2 " γ 3 γ 4 . (4.9)
Remark that here the approximation of the eigenvalues introduced in Section 3.1 shall be used.

In order to ensure the positivity of the numerical scheme, we proceed as follows: if h i`1{2 ă ε, where ε is threshold value, the viscosity matrix is computed by

P i`1{2 2 ´Ai`1 2 ¯" α 0 Id `α1 A i`1 2 (4.10)
where

α 0 " λ 3 |λ 1 | ´λ1 |λ 3 | λ 3 ´λ1 and α 1 " |λ 3 | ´|λ 1 | λ 3 ´λ1 .
This means that in a given interface where a wet-dry front is found, the selected PVM at such interface coincides with the HLL scheme. The choice of HLL scheme in such situations will grant the positivity of water and sediment thickness.

Numerical treatment of the erosion-deposition source term

The second step of the method employed is carried out as follow: Let W i be the approximation computed using the numerical scheme (4.4). We shall compute W n`1 i in a semi-implicit way as fallows, first the erosion and deposition process are incorporated in the solution

The sediment flux is computed as φ b,i " F e,i ´F d,i , being F e,i " K er,i p1 ´pq, and

F d,i " K dep,i ρ i , (4.11)
where

K er,i " v s E s,i , and K dep,i " v s ρ i c i ˜0.4 ˆd Dsg ˙1.64 `1.64 ¸. (4.
12)

The values phρq n`1 i and p1 ´pqb n`1 i are defined by the following equations

225 phρq n`1 i " phρq i `∆t 1 ´p ρ b ˆKe r,i p1 ´pq b n`1 i b i ´phρq n`1 i h i K dep,i ˙, (4.13) p1 ´pqb n`1 i " p1 ´pqb i ´∆t ˆKe r,i p1 ´pq b n`1 i b i ´phρq n`1 i h i K dep,i ˙. (4.14)
Let us be introduce the following notation

a " ∆tρ b 1 ´p K er,i b i , e " ∆tρ b 1 ´p K dep,i h i , c " ∆tK er,i b i , d " ∆tK dep,i h i ,
x " phρq n`1 i and y " p1 ´pqb n`1 i , and

x 0 " phρq i and y 0 " p1 ´pqb i , whose solution is given by:

then
x " p1 `cq x 0 `a y 0 1 `c `e , (4.16)
y " d x 0 `p1 `eq y 0 1 `c `e . (4.17)

Remark 4.1. Notice that if phρq i and p1 ´pqb i are positive then phρq n`1 i and p1 ´pqb n`1 i will be positive. This is an improvement with respect a explicit approximation of the erosion-deposition terms which does not grant positivity, which would be a major drawback.

230

From equation (4.13) we get

∆ 1 ´p φ n`1 b " 1 ρ b `phρq n`1 i ´phρq i ˘, (4.18)
this relation is employed to update the values h n`1 and phρuq n`1 . Finally the friction term is computed in a semi-implicit form.

Remark 4.2. Remark that the choice of PVM scheme results in the postivity preserving property of the first step of the numerical scheme. Moreover, the semi-implicit approach for the numerical treatment of the erosion and deposition 235 source terms is also positiviy preserving. Therefore this will assure that the global scheme introduced here for the numerical resolution of system (3.4) is positivity preserving.

Extension to bidimensional case

The mathematical model (3.1) can be easily extended to the bidimensional case. For the sake of completeness we shall describe here the model and numerical scheme in the 2D case.

The model may be written in compact form as

B t W `Bx F 1 pWq `By F 2 pWq " B 1 pWqB x W `B2 pWqB y W `GpWq `Sτ pWq, (5.1)
where W " ph, hρ, hρ u, hρ v, bq t and

F 1 pWq " ¨hu hρ u hρ u 2 `1 2 gh 2 ρ hρ u v 1 1 ´p q bx ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' ; F 2 pWq " ¨hv hρ v hρ u v hρ v 2 `1 2 gh 2 ρ 1 1 ´p q by ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' .
(5.2)

B 1 pWq " ´pg hρqe 3,5 , B 2 pWq " ´pg hρqe 4,5 , (

GpWq "

¨φb 1 ´p ρ b φ b 1 ´p ρ u 2 φ b ρ v 2 φ b ´φb 1 ´p ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' ; S τ pWq " ¨0 0 τ u τ v 0 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , (5.4)
where e r,s is the canonical basis of square matrix of order 5.

The computation of φ b is made essentially using (2.39) but replacing the absolute value of the velocity |u| by the norm of the vector velocity |u|. The 245 components of friction term τ f " pτ u , τ v q are computed using (2.47) by taking into account the direction of velocity when the bed shear velocity (2.48) is computed.

The Meyer-Peter & Müller formulation for the two-dimensional solid transport flux is given by

q b " 8Q sgnpuqpθ ´θc q 3{2 `, (5.5)
where sgnpuq corresponds to u{|u|. The Shields parameter, in this case, takes the form:

θ " µ 2 f dpρ s ´ρ0 qh 1{3 |u|. (5.6)
Straightforward computations allow us to prove that the system (5.1) is invariant under ratations (see [START_REF] Fernández-Nieto | Modelling and numerical simulation of submarine sediment shallow flows: transport and avalanches[END_REF]), that is to say, given a rotation matrix

T n " ¨1 0 0 0 0 0 1 0 0 0 0 0 n x n y 0 0 0 ´ny n x 0 0 0 0 0 1 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , (5.7)
where n " pn x , n y q t is any unitary vector, if we considerer F n " F 1 pWqn x F2 pWqn y and BpWq " pB 1 pWq, B 2 pWqq then, the following relations are ful-

250 filled F n pWq " T ´1 n F 1 pT n Wq, (5.8)
BpWq ¨n " T ´1 n B 1 pT n Wq.

(5.9)

Following [START_REF] Fernández-Nieto | Modelling and numerical simulation of submarine sediment shallow flows: transport and avalanches[END_REF][START_REF] De La Asunción | Efficient gpu implementation of a two waves tvd-waf method for the two-dimensional one layer shallow water system on structured meshes[END_REF] we can see that T n W verifies the system

B t pT n Wq `Bn F 1 pT n Wq " B 1 pT n WqB n W `Qn K , (5.10)
where Q n K " T n p´B n K F n K pWq `pBpWq ¨nK qB n K Wq is the tangential term.

This will allow us to easily extend the 1D numerical scheme to this case.

Numerical solution for the bidimensional model

The numerical scheme is built by using (5.10) at each edge the volumes of mesh. In order to solve (5.1), the computational domain is divided into finite volumes V i Ă R 2 . In this work we use structured meshes made up by rectangular cells whose edges are parallel to the Cartesian axes. Given a finite volume V i , N i is the set of indexes j such that V j is a neighbour of V i ; E ij is the common edge to two neighbour cells V i and V j and |E ij | represents its length; n ij is the normal unit vector of the edge E ij pointing towards to the cell V j ; and |V i | is the area of the finite volume V i .

The path conservative well balanced numerical scheme for the two dimensional model writes

W n`1 i " W n i ´∆t |V i | ÿ jPNi |E ij |D íj `Wn i , W n j ˘, (5.11)
where D íj `Wn i , W n j ˘is an approximation of the normal flux associated to the edge E ij . By using the invariant under rotation property (5.8) we can define

D íj `Wn i , W n j ˘" T ´1 nij D ´`T nij W i , T nij W j ˘, (5.12)
where D ´`T nij W i , T nij W j ˘is the numerical flux function defined Section 4.

In this case the one dimensional system is defined from Equation (5.10) by neglecting the tangential term Q n K , further details can be seen in [START_REF] Fernández-Nieto | Modelling and numerical simulation of submarine sediment shallow flows: transport and avalanches[END_REF].

Erosion deposition and friction source terms are treated as it was done in Section 4.1.

Numerical simulations

In this section, several numerical experiments have been carried out in order to validate the model and the numerical scheme used.

Lake at rest

In this test we try to reproduce numerically the characteristics of lake at rest state, as described in Section 3.2.1, for the system (5.1). We consider a closed flume with flat bottom and we consider that the flume is initially filled with a mixture of water and sediment. The density of sediment is 1.580g{cm 3 , the diameter of the particle is 3.9 mm and the porosity is 0.47. The initial concentration of sediment is 0.2, the initial depth of bottom is 1.05 m and the initial free surface of water is 3.0 m. The total time simulated is 30 s. The domain Ω " r´1, 1s ˆr´1, 1s was discretized by using 100 points in the x and y directions, respectively. Free boundary conditions are considered in all boundaries of the domain. Figure 5 shows, as expected, that the suspended sediment is deposed at the bottom, therefore the depth of bottom increases. We can notice that the free surface of water remains constant through the time. Remark that the models presented in [START_REF] Bradford | Hydrodynamics of turbid underflows. i: Formulation and numerical analysis[END_REF][START_REF] Khan | Numerical modeling of hyperpycnal plume[END_REF] do not verify this property. In those models, under lake at rest hypothesis, the freshwater mass is increased due to lack of a of a source term in the first equation, which makes that the free surface level does not remains constant through time. Figure 6 shows that the density of water sediment mixture decreases because we have selected a initial value of c such that the relation c ă 1 ´p is fulfilled. We can conclude that the numerical scheme reproduces correctly the conditions of the lake at rest state for the system (5.1).

Dune evolution

In this test we shall consider bedload transport for the morphological evolution of a dune. We shall consider clear water at initial time and we shall impose an entrance of suspended sediment concentration at left boundary. In this experiment the erosion and deposition terms are neglected. We consider a sedimentary material with the same characteristics described in Section 6.1 and the same domain with the same discretization. The initial conditions are given by bpx, 0q "

$ & % 1.1 `0.05 cosp10|x|q, if |x| ă π 10 , 1.05, Otherwise. (6.1)
hpx, 0q `bpx, 0q " 2.0, hpx, 0q upx, 0q " 1, hpx, 0q vpx, 0q " 0, ρpx, 0q " ρ 0 .

(6.2)

We have employed the Meyer-Peter & Mülller formulation to compute the bedload discharge (5.5), but in this case we have used

Q " 2d c g ρ s ´ρ0 ρ o d
in order to see movement in the bottom at short times. We have imposed the following boundary conditions, hpx, tq upx, tq " 1, for x P r´1, 1s ˆt´1, 1u, t ą 0, hpx, tq vpx, tq " 0, for x P t´1u ˆr´1, 1s, t ą 0.

(6.3)

We shall impose a homogeneous boundary condition on suspended sediment concentration given by ρpx, tq " ρ 0 `0.019 pρ s ´ρ0 q hpx, tq , for x P t´1u ˆr´1, 1s, t ą 0. (6.4)

The remaining boundaries are considered with open boundary condition.

The final time is t " 40 s. The numerical results obtained at time t " 40 s can be seen in Figure 7.

We can see that the numerical scheme preserves the symmetry, with respect of x axe, of the water free surface, bottom, and discharge in x direction. Remark that the provided numerical solution is smooth. We can see that the bottom at the final time presents a star-shape form which agrees to what is expected (see [START_REF] Vriend | 2dh mathematical modelling of morphological evolutions in shallow water[END_REF]).

Simulations of one-dimensional dam breaks flows

In this section one dimensional experimental dam breaks are simulated in order to do a comparison between the experimental data and the numerical results.

Experiment 1

This dam break experiment was performed in the University of Taiwan, Taipe and it was reported in [START_REF] Capart | Formation a jump by the dam-break wave over a granular bed[END_REF]. In this experiment the sediment particles considered were artificial spherical pearls covered with a shiny coating, with a diameter of d " 6.1 mm, density ρ s " 1.048g{cm 3 , porosity p " 0.4, roughness manning coefficient µ f " 0.03, critical Shields parameter θ c " 0.045 and final time T " 5 t 0 , where t 0 " a h 0 {g is the hydrodynamic time scale, h 0 being the initial upstream depth of water. The domain Ω " r´0.6, 0.6s is divided in 500 points and the initial conditions are given by hpx, 0q "

(a) (b) (c) (d)
$ & % 0.1, if x ď 0, 0, otherwise,
upx, 0q " 0, ρpx, 0q " 1, bpx, 0q " 0. (6.5)

Open boundary conditions are imposed in both ends of the channel. This is done by simply using the ghost cell technique.

A comparison between the numerical results obtained and the experimental data can be seen in Figures 8 to 10. From these figures we can see that the front of the wave is well captured. We notice that at time 3 t 0 the profile of free surface of water is under-estimated while for times 4 t 0 and 5 t 0 the profile computed depicts successfully the shape of free surface provided by the experimental data. On the other hand, for the times 3 t 0 and 4 t 0 the erosion rate computed under-estimate the observed one in the numerical data, while for the time 5 t 0 the erosion rate agrees with the experimental one. As it was pointed out by [START_REF] Capart | Formation a jump by the dam-break wave over a granular bed[END_REF], the non-hydrostatic effects, which are not taken into account by this model, can explain the discrepancies observed between the numerical results and the experimental data. The differences observed in the free surface may be explained by the rapid drop in the hydrostatic pressure. For instance, better approximations of the free surface for this case are achieved in [START_REF] Cantero-Chinchilla | Nonhydrostatic dam break flows. II: one-dimensional depth-averaged modeling for movable bed flows[END_REF] where a nonhydrostatic model is employed. The counterpart being that a non-hydrostatic model is computationally more expensive.

Experiment 2

A small scale dam break flux was performed in the Université catholique de Louvain. This experiment was reported in [START_REF] Fraccarrollo | Riemann wave description of erosional dambreak flows[END_REF]. The sediment considered was composed by cylindrical PVC pellet with diameter " 3.2 mm and height " 2.8 mm which is equivalent to a spherical diameter " 3.5 mm, the density ρ s " 1.540g{cm 3 , the porosity is p " 0.4 , the roughness coefficient µ f " 0.03, the critical Shields parameter θ c " 0.045 and the total time T " 10 t 0 , where t 0 " a h 0 {g is the hydrodynamic time scale, h 0 being the initial upstream depth of water . The domain Ω " r0, 2.5s is discretized using 500 points and the initial conditions are given by hpx, 0q "

$ & % 0.1 if x ď 1.25, 0 Otherwise, upx, 0q " 0, ρpx, 0q " 1, bpx, 0q " 0. (6.6)
Open boundary conditions are considered at upstream and downstream ends of 335 the channel.

Figures 11 to 13 show that the numerical results capture successfully the advancing front of the water. The time evolution of bottom computed for times t " 5 t 0 and t " 7.5 t 0 has a high degree of accuracy. The free surface for these computed times shows a good level of accuracy as well when compared with the 340 reported experimental free surface. For time t " 10 t 0 we can see that the upstream rate of sediment entrainment computed is greater than observed one in the experimental data,while the computed downstream rate of sediment entrainment approximates with a good level of accuracy the experimental data. Again, differences between the shape of the computed free surface and the reported 345 one are associated with the hydrostatic pressure considered in the mathematical model. We remark that the numerical results obtained in this experiment improve the results obtained for instance in [START_REF] Wu | One-dimensional modeling of dam-break flow over movable beds[END_REF][START_REF] Liu | A robust coupled 2-D model for rapidly varying flows over erodible bed using central-upwind method with wetting and drying[END_REF].

Experiment 3

This small scale dam break experiment was reported in [START_REF] Spinewine | Small-scale laboratory dam-break waves on movable beds[END_REF]. In this experiment the sediment taken into account is characterised by a diameter of particle d " 3.9 mm, density ρ s " 1.580 g{cm 3 , porosity p " 0.47, the roughness coefficient µ f " 0.0165, the critical Shields parameter θ c " 0.047 and the total time of compute T " 1.5 s. The domain for this experiment is Ω " r´3.0, 3.0s which has been divided into 500 points. The initial condition are give by hpx, 0q "

$ & % 0.35 if x ď 0, 0 Otherwise, upx, 0q " 0, ρpx, 0q " 1, bpx, 0q " 0. (6.7)

Open boundary conditions are considered at the upstream and downstream ends of the channel.

In Figures 14 to 16 a high level of precision of the numerical results is observed. The shape of computed free surface recreates outstandingly the free surface given by the experimental data at all computed times. The computed shape of bottom approximates successfully the bottom reported by the data, for all computational times.

Simulations of biimensional dam-breaks flows

In this section we recreate numerically some bidimensional dam breaks experiments which have been reported in the literature.

Experiment 4

This experiment was carried out at the Hydraulics Laboratory of Université catholique de Louvain and it was reported in [START_REF] Goutiere | Dam-break flow on mobile bed in abruptly widening channel: experimental data[END_REF]. The sediment used was coarse, almost uniform sand with a median diameter of 1.72 mm, density ρ s "

2.63 g{cm 3 and porosity of 39%, the roughness coefficient is µ f " 0.025. A fully saturated and compacted sand layer of thickness 0.10 m over whole flume is considered as initial condition for the bottom. The upstream condition for the clear water layer is 0.25 m and 0 m for the downstream. In Figure 17 we can see a sketch of the shape of the domain and the considered initial conditions. A free boundary condition is considered at right end of the channel while a wall boundary condition is assumed at remaining boundaries. Figure 22 shows the comparison of temporal evolution of free surface of water, we can see that the the evolution computed by the numerical results follows the shape of the temporal evolution given by the experimental data, little 390 difference can be observed after the eighth second. In Figure 23 we can notice that before t " 3 s the numerical results depict correctly the time evolution of the experimental data, between t " 3 s and t " 5 s the numerical results under-estimate the values reported by the data, after t " 6 s we can see an over-estimation gives by the numerical results. It can be observed that the 395 approach presented here improves the numerical results obtained in [START_REF] Swartenbroekx | Two-dimensional two-layer shallow water model for dam break flows with significant bed load transport[END_REF] for the bed evolution and they are comparable to the ones obtained in [START_REF] Liu | A robust coupled 2-D model for rapidly varying flows over erodible bed using central-upwind method with wetting and drying[END_REF].

Experiment 5

This experiment has been conducted at the Hydraulics Unit of Mechanical and Civil Engineering Laboratory, Université de Louvain. The set-up of the 400 flume and its dimensions can be seen in Figure 24. This experiment was reported in [START_REF] Soares-Frazão | Dam-break flows over mobile beds: experiments and benchmark tests for numerical models[END_REF]. The sediment used was coarse sand characterized by a diameter d "

1.61 mm with density ρ s " 2.63g{cm 3 and porosity p " 0.42, the manning roughness coefficient reported was µ f " 0.0165 s{m 1{3 for the movable bottom and µ f " 0.010 s{m 1{3 for the fixed bed. The boundary conditions consisted of 405 a closed wall at the upstream of the flume and free boundary condition at the downstream of the flume. The remaining boundary conditions correspond to wall condition. The initial water level in the upstream is denoted by η 0 , while the initial water level in the interval 0 ă x ă 9 is defined by η 1 , and initial is h s " 0.085 m, these initial conditions can be seen in Table 1. During the experiment test repeatability was cheked to measure its level of repeatability, further details can be found in [START_REF] Soares-Frazão | Dam-break flows over mobile beds: experiments and benchmark tests for numerical models[END_REF]. between the numerical results and the experimental data for the final bottom at different longitudinal sections. In Figure 26 we can see that the morphological evolution is well recreated by the numerical results. Notice that in the interval 0 ă x ă 2 the rate of erosion in the numerical results is greater than the observed one in the experimental data, but at the remaining comparison section 420 the numerical results capture successfully the shape of experimental data. Fig-

ure 27 shows that the general shape of the experimental data is recovered by the numerical results, little differences can be observed at 1 ă x ă 2, but the approximations computed have a good level of accuracy. In Figure 28 we can observe that the numerical results recreate successfully the experimental data.

425

The deposition rate is well computed in the section 0 ă x ă 2, after this interval the deposition rate is under-estimated (or the erosion rate is over-estimated).

During the experiment, the temporal evolution of water level was measured and the gauge positions are stated in Table 2. In Figures 29 to 32 we can see in a good way the experimental data during the five second, after that time the numerical results over-estimated the experimental data. In [START_REF] Soares-Frazão | Dam-break flows over mobile beds: experiments and benchmark tests for numerical models[END_REF] the numerical results provided by several models were compared with the experimental data.

These numerical results gave a good approximation for the bed evolution and the water level. In comparison with these models, the approach presented here improves the results at the sections y " 0.6 and y " 1.45 for the bed evolution.

The results for the water level are similar.

Conclusions

We have deduced a new model for suspended and bedload sediment transport under the shallow water hypothesis. The model does not make a Boussinesq type hypothesis. The numerical scheme proposed is based on a two-step algorithm. First the underlying hyperbolic system is solved using a PVM path conservative scheme, which is a fast and robust scheme. Then the source terms corresponding to erosion and deposition are applied. This is done by using a semi-implicit approach which allow us to guarantee that the computed densities remain positive, so that it results in a robust way of treating such terms.

From the comparison done between the numerical results and the experimental data we can conclude that the mathematical model and the numerical scheme proposed in this work provided successful results. Notice from the one dimensional experiments that a high level of accuracy in the results is obtained.

The model captures correctly the advance of the front wave and the shape the of free surface. Regarding to the morphological evolution, a high level of accuracy is observed as well. Some differences may be observed, specially for small times, between the numerical and experimental results. These discrepancies may be explained by the presence of non hydrostatic effects which are not considered in the model. Doing so would make it more cumbersome and less efficient.

For the two dimensional experiments, on the one hand we can conclude that the morphological evolution is computed successfully, the rate of erosion and deposition are according with the observed ones in the experimental data, but the amplitudes were underestimated. On the other hand, the free surface of water was good enough computed in some gauge while for other a overestimation of that is obtained.

In general, the mathematical model and the numerical scheme proposed here is able to recreate transient flows over movable bottom in an accurate and efficient way.

F

  d px, tq ηpx, tq " hpx, tq `bpx, tq, free surface of water

Figure 1 :

 1 Figure 1: Sketch of bed load and suspended-load sediment transport.

r φ b ´Bt ˆη b p ρp udz ´Bx ˆη b p ρp u 2 dz ´g ˆη b ˆBx ˆη z p ρdz 1

 1 ) we get the averaged momentum equationB t phρ uq `Bx hpρ u 2 q `Bx ˆg h 2 2 ρ ˙" ´ghρB x b`µpr τ | z"η ´r τ | z"b q `ρu| z"b ˙dz.

is 7 6 u

 6 Figure 2.

Figure 2 :

 2 Figure 2: Sketch of the proof for the hyperbolicity.

Figure 3 :

 3 Figure 3: Sketch of the procedure to compute the characteristics velocities of (3.4).

. 23 )

 23 Moreover the volume and mass conservation equation in (3.4) results in 165 B x phuq " 0, (3.24) B x phρuq " 0. (3.25)

Figure 4 :

 4 Figure 4: Finite volume mesh.

  the system of equations (4.13)-(4.14) reduces to

Figure 5 :

 5 Figure 5: Free surface and bottom at time t " 0 s (left) and free surface and bottom at t " 30 s (right).

Figure 6 :

 6 Figure 6: Density at time t " 0 s (left) and density at t " 30 s (right).

Figure 7 :

 7 Figure 7: Numerical results at time t " 40 s; (a) Free surface, (b) Bottom, (c) Discharge on x direction and (d) Discharge on y direction.

Figure 8 :

 8 Figure 8: Comparison between numerical results and experimental data at time t " 3 t 0 s.

Figure 9 :Figure 10 :

 910 Figure 9: Comparison between numerical results and experimental data at time t " 4 t 0 s.

Figure 11 :

 11 Figure 11: Comparison between numerical results and experimental data at time t " 5 t 0 .

Figure 14 :Figure 15 :

 1415 Figure 14: Comparison between numerical results and experimental data at time t " 1.0 s.

Figure 16 :

 16 Figure 16: Comparison between numerical results and experimental data at time t " 1.5 s.

Figure 17 :

 17 Figure 17: Experimental set-up and initial conditions: (a) Longitudinal profile and (b)plan view.

Figure 18

 18 Figure 18 shows the free surface of the water and the depth of bottom at final simulation time t " 12 s. In Figures 19 to 21 we can see a comparison between the numerical results and the experimental data for the depth of bot-

Figure 20

 20 shows that the numerical results under-estimate the experimental data for y " 0.25 m, after that point the numerical results recreate good enough the tendency of the experimental data. The comparison at cross section x " 1.6 m is showed in Figure21, we can see that the numerical results385adjust successfully the experimental data. Therefore we can conclude that the morphological evolution for this experiment is computed successfully.

Figure 18 :

 18 Figure 18: Free surface at t " 12 s (left) and bottom at t " 12 s (right).

  0

Figure 22 :Figure 23 :

 2223 Figure 22: Temporal evolution of free surface at point p4.2, 0.125q.

Figure 24 :

 24 Figure 24: Experimental set-up, flume dimensions: (a) plan view, (b) elevation and (c) crosssections.

Figure 25 showsTable 1 :Figure 25 :

 25125 Figure 25 shows the final state for the free surface of water and the depth of bottom after 20 s of simulation. In Figures 26 to 28 we can see the comparison415

Figure 26 :Figure 27 :

 2627 Figure 26: Bottom surface at section y " 0.2 at time t " 20 s.

Figure 28 :

 28 Figure 28: Bottom surface at section y " 1.45 at time t " 20 s.

Figure 32 :

 32 Figure 32: Temporal evolution of free surface at point US6.

  Comparison between numerical results and experimental data at time t " 7.5 t 0 . Comparison between numerical results and experimental data at time t " 10 t 0 .

		0.1
								Numerical free surface
		0.08						Numerical bottom Experimental data
	b(m),FS(m)	0.04 0.06
		0.02
		0
	-0.02 Figure 12: 0.6 Figure 13:	0.8	1	1.2	1.4 x(m)	1.6	Numerical free surface Numerical bottom Experimental data 1.8 2	2.2

  Cross section at x " 1.3 m at time t " 12 s.

		0.125		Numerical results
			0.12		Experimental data
	Bottom (m)	0.1 0.105 0.11 0.115
		0.095
			0.09
		0.085
				0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5
								Y (m)
	Bottom (m)	0.09 0.095 0.1 0.105 0.115 0.12 0.11	0.05 Figure 20: 0 0.1 0.15 Numerical results 0.2 Experimental data	0.25	0.3	0.35	0.4	0.45	0.5
								Y (m)
		0.11		Numerical results
					Experimental data
	(m)		0.1
	Bottom	0.08 0.09
		0.07
				0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5
								Y (m)
					Figure 19: Cross section at x " 1.15 m at time t " 12 s.
								36

Figure 21: Cross section at x " 1.6 m at time t " 12 s.

Table 2 :

 2 Gauge locations. the experimental data and the numerical results for the temporal evolution of the water level. We can see in Figures29 and 30that at gauge US1 and US3 that the numerical results depicts successfully the temporal evolution of the water level provided by the experimental data. On the other hand from Figures 31 and 32 we notice that the numerical results approximate

	Gauge Id xpmq ypmq
	US1	0.640 -0.500
	US3	0.640 0.165
	US5	1.940 -0.990
	US6	1.940 -0.330
	the comparison between

  Figure 30: Temporal evolution of free surface at point US3.

			0.28
			0.26
			0.24
			0.22
	FS (m)	0.16 0.18 0.2
			0.14
			0.12													Numerical results
			0.1													Experimental data
			0.08	0	2	4	6	8	10	12	14	16	18	20
															t (s)
			0.2
		0.18
	(m)	0.16
	FS	0.14
		0.12
																	Numerical results
			0.1														Experimental data
		0.08	0		2		4		6		8		10	12	14	16	18	20
															t (s)
									Figure 31: Temporal evolution of free surface at point US5.
			0.22
	0.22
			0.2
		0.2
			0.18
	0.14 0.16 0.18 0.14 FS FS (m) 0.16 (m)
	0.1 0.12 0.12 0.1												Numerical results Numerical results Experimental data Experimental data
	0 0.08 0.08	0	2	2	4	4	6	6	8	8	10 t (s) t (s) 10	12 12	14 14	16 16	18 18	20 20
									Figure 29: Temporal evolution of free surface at point US1.

iment is progressively deposed. Under water at rest supposition pu " 0q, the