Henrique Guilherme

Fabiana Fiorot
email: gfiorot@insa-rennes.fr

Oliveira De

Pascal Ferreira

Geraldo Dupont

De Freitas Maciel

Guilherme Henrique Fiorot

Fabiana De Oliveira Ferreira

Pascal Dupont

Geraldo De Freitas Maciel

Roll waves study on Acquabona watershed

Keywords: Open channel flow, Mudflow, Natural risks, Roll waves

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Natural hazards are a subject each day more current and that has been broadly disclosed on news and media. As results from global warming, landscape usage and bad urban management, the climate change emerges as worldwide events, thus leading to natural catastrophes all across the globe [START_REF] Burton | The social construction of natural disasters: An evolutionary perspective[END_REF]. Those catastrophes are result of the impact of natural phenomena over vulnerable environments. In fact, they are in fact the conjunction of meteorological, geological/geomorphological effects with anthropic actions. Also affected by climate changes, Brazil has confronted successive natural events such as floods and landslides during the rainy season [START_REF] Pinho | Analise espacial dos movimentos de massa em Nova Friburgo/RJ: o caso do desastre natural de janeiro de 2011[END_REF]. On Brazilian South and Southeast regions, mudflows are especially prominent during the rainy season over Serra do Mar, the mountains chains that extends along the shore [START_REF] Guidicini | A Herschel-Bulkley model for mud flow down a slope[END_REF]. According to Brazilian Ministry of Science and Technology, landslides and mudflows are the greatest cause of death among all Brazilian natural hazards.

Apart from the great damage incurred during Rio de Janeiro's event, the scenario established in the event (topography, mud constitution, etc.) was favorable to roll waves instabilities generation. However, properties such as front velocity propagation, mudflow/landslide depth from events of such nature haven't been recorded by Brazilian authorities and/or researchers which add up difficulty to any further investigation/conclusion. To analyze such kind of events based on mathematical models, the work here presented will explore an analogous debris/mud flow event recorded by [START_REF] Berti | Debris flow monitoring in the Acquabona watershed on the Dolomites (Italian Alps)[END_REF] on Acquabona watershed, Cortina d'Ampezzo, Italy. The aim here is to identify how good a simplified mathematical model is on predicting roll waves formation on a natural environment as the presence of such phenomenon could be translated as a different level of risk to the environment.

From a mechanical point of view, the representation of mass-movement can be done through many physical approaches [START_REF] Iverson | The physics of debris flows[END_REF]. To choose a specific model over another relies on the representativeness of that model to the target phenomenon. Among them, mudflows and mud floods are natural events well represent by hydraulic approaches. They are known by their ability to erode soil and entrain sediment, sometimes even changing its rheological behaviour. From this last perspective, fluid in mudflows/mud floods can be Newtonian (clean water, low concentration of sediment, less than 1% in volume) or non-Newtonian (hyperconcentrated, high concentration of sediment, greater than 8% in volume). The importance of the correct approach to risk management is clear when prediction of values such as flow height, front velocity, flow runout are at stake. Such task has been well developed, but still lack of application.

Besides all those aspects, secondary effects a like bank erosion, front propagation and roll waves instabilities could make natural hazard events even more dangerous as they can change flow features and consequently interact more with the environment. The last, roll waves, are periodic free surface instabilities with considerable energy level (i.e. high amplitude and propagation velocity) which might appear on open channel flow configuration such as those we perceive on mountain regions and cliffs. Natural landscape can constitute slopes and preferable path for flows. Slope and fluid rheology could then reach (or surpass) threshold values for roll waves generation. In the past 10 years, several mathematical models have been developed on mud and debris flows, adding up complex features such as arbitrary topography, non-Newtonian fluid rheology and even assuming unstable conditions of the free surface [START_REF] Van Asch | Techniques, issues and advances in numerical modelling of landslide hazard[END_REF][START_REF] Balmforth | Dynamics of roll waves[END_REF][START_REF] Zanuttigh | Instability and surge development in debris flows[END_REF].

Roll waves appearance on artificial channel sections can lead to periodic overflows what induce debris entrainment and banks destabilization. From this perspective, they could not only elevate probability of mudflows events but also increase social and economic risks inherent to those hazards. Although no substantial proof of roll waves appearance have been taken on Brazil's territory, countries like New Zealand, China, Switzerland and Italy have already confirmed presence of surges and wave trains on mud flows over mountain regions [START_REF] Zanuttigh | Instability and surge development in debris flows[END_REF].

Especially on Italy, an event at Cortina d'Ampezzo, on August 17th, 1998, draws attention regarding wave train detection, formed by almost periodic surges of mud and debris flow. This event, reported by [START_REF] Berti | Debris flow monitoring in the Acquabona watershed on the Dolomites (Italian Alps)[END_REF], will be explored on this paper to promote a realistic data confrontation and a validation attempt to mathematical model described by [START_REF] Maciel | Control of instabilities in non-Newtonian free surface fluid flows[END_REF], which will be recalled on this work. Identify formation conditions and determine how roll waves instabilities will develop are task to be performed to better understand the risks dependence on such instabilities.

It's important to remark that the correct representation of fluids in this kind of phenomena requires attention. The variability found on solid-liquid fraction that composes the mixture contributes to the complexity of the fluid which could acquire cohesive or non-cohesive behaviour, depending on the material in suspension. Thus, exactly quantify shear ratio suffered from fluid due to stress applied is not usually possible. Mixtures where sediment concentration is high will be characterized mainly by the presence of a yield stress and viscous interactions, usually called mud flows. On the other hand, when dispersive turbulent stress are remarkably noticed and become the main mechanism of energy dispersion characterize mud floods and debris flow events.

The variability of physical properties found in situ during or even after events (mud floods, mud flows, debris flow) gives most of the difficulty that arises when categorizing those events. Conversely, even though mathematical approximations and modelling are usually employed, often they don't properly represent realistic scenarios. Usually data employed on those models come from experimental data and/or in situ measurements sometimes poorly collected, i.e. their representation of wholeness of event, environment, etc. are not complete. This paper will show the applicability of a 1-D model, based on shallow water equations and viscoplastic rheology, to that kind of event. No turbulent dispersion should be considered. Roll waves properties estimative (amplitude, wavelength, propagation velocity) will be made. The event used in this paper was reported by [START_REF] Berti | Debris flow monitoring in the Acquabona watershed on the Dolomites (Italian Alps)[END_REF] and make reference to a 3 debris flow event from 1998 in Acquabona watershed, Cortina d'Ampezzo, Italy. However, authors did not have quantified the contribution of turbulent dispersions to the stress tensor.

MATHEMATICAL AND NUMERICAL MODELING OF ROLL WAVES

The modeling departs from a mechanical point of view, considering a constant discharge rate (𝑞 [m 2 /s]) of mud and debris aqueous mixture with volumetric concentration 𝜙 𝑣 , composing an incompressible viscous fluid with density 𝜌 [kg/m 3 ] flowing down an inclined open channel (𝜃 degrees), where the driven force is gravity, represented by its acceleration (𝑔 [m 2 /s]). For the compounded mixture (water + clay or water + fine sand + clay) investigated in this work, rheometric studies from [START_REF] Coussot | Steady, laminar, flow of concentrated mud suspensions in open channel[END_REF], [START_REF] Piau | Flow of a yield stress fluid in a long domain -Application to flow on an inclined plane[END_REF], Huang & Garcia (1998), [START_REF] Maciel | Rheological analysis of water clay compositions in order to investigate mudflows developing in canals[END_REF] and many others proved that when in absence of sedimentation effects, fluid rheology can be described by the Herschel-Bulkley non-linear rheological model with permanent simple shearing conditions, as shown by the following expression:

𝜏 𝑥𝑧 = 𝜏 𝑐 + 𝐾 𝑛 ( 𝜕𝑢 𝜕𝑧 ) 𝑛 , for 𝜏 𝑥𝑧 > 𝜏 𝑐 , (1)
elsewhere 𝜕𝑢 𝜕𝑧 ⁄ = 0. On Eq. 1, 𝜏 𝑥𝑧 is the shear stress magnitude acting on longitudinal coordinate 𝑥 due to a gradient in vertical coordinate 𝑧; 𝑢 is the longitudinal velocity, 𝜕𝑢 𝜕𝑧 ⁄ is the shear rate magnitude. This model is then inserted into the viscous part of the stress tensor on Cauchy's equation.

The flow is then considered laminar, where the Reynolds number 𝑅𝑒 = 𝜌(𝑢 0 2-𝑛 ℎ 0 𝑛 )/𝐾 𝑛 should be not greater than 500, and in shallow water conditions, where mathematical approximations lead to a Saint-Venant alike system. Equations are then vertically averaged and mean values for flow height (ℎ 0 [m]) and velocity (𝑢 0 [m/s]) are found. Through these values, a dimensionless system is constituted where three main dimensionless numbers are identified: the Froude number 𝐹𝑟 = 𝑢 0 /√𝑔ℎ 0 cos 𝜃 ; the dimensionless yield stress 𝐶 * = 𝜏 𝑐 /(𝜌𝑔ℎ 0 sin 𝜃); and the flow index 𝑛. The roll wave equation is given in a mobile coordinate system (𝑥 ′ = 𝑥 * -𝑈 * 𝑡 * ), as follows:

where 𝑈 is roll wave celerity, computed as:

𝑈 * = 𝑈 𝑢 0 = 𝛼 + √𝛼 2 -𝛼 + 𝐹𝑟 -2 , ( 3
)
And 𝛼 is a momentum correction factor that assumes values between 1 and 1.2 (for complete plug flow and pure laminar, respectively). This equation should have continuous solution inside the wavelength if sufficient conditions are reached. The domain of possible solution depends on flow Froude number, which should be greater than the threshold 𝐹𝑟 𝑚𝑖𝑛 , found through linear stability analysis (Trowbridge, 1987;[START_REF] Coussot | Steady, laminar, flow of concentrated mud suspensions in open channel[END_REF][START_REF] Maciel | Control of instabilities in non-Newtonian free surface fluid flows[END_REF] and given by the Equation (4):

𝐹𝑟 > 𝐹𝑟 𝑚𝑖𝑛 = 𝜗 √𝜑 2 -2𝛼𝜑𝜗 + 𝛼𝜗 2 , ( 4
)
where 𝜑(𝑛, 𝐶 * ) = (𝑛 + 1)(2𝑛 + 1)/(𝑛 + 1 + 𝑛𝐶 * ) and 𝜗(𝑛, 𝐶 * ) = 𝑛(1 -𝐶 * ).

𝜕ℎ * 𝜕𝑥′ = ℎ * -𝐶 * -(1 -𝐶 * ) [(1 + 𝑈 * (ℎ * -1)) ( 1 -𝐶 * ℎ * -𝐶 * ) ( 𝑛 + 1 + 𝑛𝐶 * (𝑛 + 1)ℎ * + 𝑛𝐶 * )] 𝑛 (𝛼 -1)𝑈 * 2 - 𝛼(1 -𝑈 * ) 2 ℎ * 2 + ℎ * 𝐹𝑟 2 . ( 2
)
Wavelength 𝜆 is calculated through defined integration of Equation ( 2) inverted, between wave trough and crest (ℎ 1 and ℎ 2 , respectively). Then, ℎ * (𝑥′) is given by defined integration of Equation ( 2) inside wavelength 𝜆. Finally, dimensionless results are brought to scale through mean flow parameters in steady and permanent configuration. All steps and details of this mathematical model can be found at [START_REF] Maciel | Control of instabilities in non-Newtonian free surface fluid flows[END_REF]. Schematically, roll waves properties can be described as shows Figure 1: Finally, the model utilization can be summarized as shows Figure 2. This methodology will be adopted in this article to estimate roll waves formation on a realistic case scenario, based on the event detailed by [START_REF] Berti | Debris flow monitoring in the Acquabona watershed on the Dolomites (Italian Alps)[END_REF]. The model here applied was only calibrated to fluid materials having silt, clay and sand in its matrix, corresponding mostly to the granulometry variability of downstream material related during the chosen Acquabona event.

Figure 2. Diagram showing input and output data of roll waves prediction methodology. 𝑇 𝑐 and 𝜆 𝑐 are characteristic time and length used for scaling process, following [START_REF] Maciel | Control of instabilities in non-Newtonian free surface fluid flows[END_REF].

MODELING ACQUABONA EVENT

The purpose of this work is to explore in situ measurements collected by [START_REF] Berti | Debris flow monitoring in the Acquabona watershed on the Dolomites (Italian Alps)[END_REF] and find out if roll waves were prone to develop under those circumstances and what would have been their main features. Based on the mathematical model recalled before, this work will estimate roll waves properties that could have happened on Acquabona site based on basic information given by those authors such as precipitation, landscape, sediment material displaced, etc.

Site description

The Acquabona channel is localized on the Dolomites Mountains, Northeast Italy, Southeast from Cortina d'Ampezzo city, having downstream the Boite river valley. Morphopedologic properties of region are prone to mass movement natural events that occur almost every year during rainy season and for this reason is explored by researchers for studies. The upstream region is composed by a limestone rock basin, presenting 30° steepness. Along the channel, the superficial soil is composed by rock and gravel from Dolomites and present mean inclination of 18°. Downstream region reaches Boite river, crossing road Roma, and presents 7° steepness. The total length of the channel is 1632 meters.

The fluid formed from the mass movement changes as it flows down the channel. The viscous matrix is the result of and aqueous solution of gravel (with mean diameter between 2 and 20 mm) and sand (diameter less than 2 mm). The flow is capable of dislodge rocks and boulders up to 2 meter diameter. Fluid fraction of rock with diameter greater than 20 mm can reach 28% of total mixture mass.

During flow evolution from upstream to downstream, a densification process is remarked by fluid samples taken from past events (Berti et al;[START_REF] Berti | Field observations of a debris flow event in the Dolomites[END_REF][START_REF] Berti | Debris flow monitoring in the Acquabona watershed on the Dolomites (Italian Alps)[END_REF] which can be explained by soil composition of intermediate part of the channel that is more susceptible to be eroded. Upstream, roughly 60% of fluid matrix is gravel while rough sands (from diameter between 65 µm and 2 mm) and fines (silt and clay with diameter smaller than 65 µm) have 35% and 5% of its composition. At intermediate length, fines quantities are greater, reaching 20% in volume of fluid matrix. Downstream, fines reach 30% of mixture, pointing out channel erodibility.

Measurements stations are installed at strategic position along the channel, allowing data acquisition from events. Data collected comprise pluviosity, wind speed, surges velocity (geophones), flow height (ultrasonic devices), flow superficial velocity (image processing) and pressure gauges. Manipulation of these data will make possible estimative of fluid rheological properties and flow dynamics (discharge rate) which along site topography constitute all necessary data input for model.

Event sample

During summer of 1998, at August 17th, a large debris flow event occurred on Acquabona channel when data were collected by 3 measurement stations installed along the channel.

Stations were positioned at specific points of the channel to facilitate comprehension of event evolution. The first station, installed 300 meters from channel inlet (rock basin), allowed observations from cameras (image data, later processed into surface velocity information), pluviosity, surge velocity and pressure from the event. The second station, positioned on the intermediate part of the channel, just wind speed and surge velocity data were collected. The last station, installed 500 meters away from channel outlet close to Boite river, image data, flow height surges velocities were collected. All data, methods and descriptions are related on Berti et al work [START_REF] Berti | Debris flow monitoring in the Acquabona watershed on the Dolomites (Italian Alps)[END_REF]. Measurements from pluviosity done at inlet station showed that the rain that initiated the event reached 25 mm during 30 minutes and that the first debris surge happened after 45 minutes of rain, when pluviosity reached 20 mm. The total volume displaced estimated was 8000-9000 m 3 .

Startup mechanisms of bed mobilization were elucidated and discussed by the authors in details. However, the observation concerning the event of August 17th, 1998, was particularly interesting as more than 15 surges of mud and debris were identified during 35 minutes after the first surge, as shows figure 1. Nevertheless, the stage were the wave trains is detected is highly unpredictable and waves with non-constant amplitude, celerity and wavelength are measured. This could be explained as we deal with a natural phenomenon where a broad set of frequencies could have disturbed the system, thus altering waves on its transitional stage. Another important remark concerns the equipment used for measurements which has great sensibility depending on the roughness of the surface (ultrasound system versus shock wave). We should then identify a well-behaved set of waves which could be a good sample of a roll wave and gather from it all input parameters necessary for its numerical simulation and prediction.

Figure 3. Measurements from station 3, downstream channel, collected through ultrasonic equipment on August 17th, 1998, related by [START_REF] Berti | Debris flow monitoring in the Acquabona watershed on the Dolomites (Italian Alps)[END_REF]. Highlighted sequence of surges identifies possible roll wave formation used in this work.

The fronts identified by A, B, C, D and E on Figure 3 present relatively well-behaved amplitude and wavelength when compared others event fronts. As known, ultrasonic devices suffer from great losses of signal when sharp/steep properties are at stake. Then, considering this system flaw, measurements could have suffered from an intense scattering when larger boulder passed through, thus showing high peaks followed by low valleys as shows Figure 3 in the end of highlighted area where no measurement are registered. Therefore, we will identify the end of our wave train on 1209 seconds, where this elevated peak abruptly appears on data series. Characteristics of surges are then collected directly from data and showed on Table 1.

Parameters gathering

First of all, results collected from pressure gauges showed that fluid pressure was equal to the bottom shear stress, thus leading to the conclusion that flow was liquefied. This observation allows employment of rheology models for mixture representation and obtaining of others flow characteristics as will follow.

Secondly, qualitative observations and data collected from soil show that sediment material density (from rocks, gravel and sands) should be around 2000 kg/m 3 . Concerning mixture composition, according to [START_REF] Zanuttigh | Instability and surge development in debris flows[END_REF], the concentration was roughly estimated by 60% of sediment in water expressed in terms of mass. Then, from mass concentration, we can compute volume concentration which should be around 43%, and the mixture density on 1430 kg/m 3 . As explained by [START_REF] Ancey | Plasticity and geophysical flows: A review[END_REF], mixtures mainly composed by clay (colloidal mixtures) have strong non-Newtonian behavior for which the yield stress exponentially increases with material concentration. However, when sand takes part into the mixture, at high shear rate, inert particles of sand take part into the flow, configuring sediments inertial regimen [START_REF] Bagnold | Experiments on a gravity-free dispersion of large solid spheres in a Newtonian fluid under shear[END_REF], thus changing rheology behavior of mixture.

Experimental analysis shows that when concentration of colloidal material is kept and concentration of inert material is added up, non-Newtonian effects are reduced. In fact, this last mixture represents how the fluid should have behaved on the case-scenario of Acquabona. Based on experimental rheological studies [START_REF] Coussot | Steady, laminar, flow of concentrated mud suspensions in open channel[END_REF][START_REF] Ancey | Plasticity and geophysical flows: A review[END_REF][START_REF] Maciel | Rheological analysis of water clay compositions in order to investigate mudflows developing in canals[END_REF], fluid yield stress (𝜏 𝑐 ) is roughly estimated on 200 Pa. Concerning index flow, results from the same author show that it stabilizes around 0.30, which is not far from 1/3 calibrated for French torrential lava [START_REF] Coussot | Steady, laminar, flow of concentrated mud suspensions in open channel[END_REF]. Although difficulty to predict rheological parameters, consistency index 𝐾 𝑛 can be estimated and should assume values between 50 and 1000 Pas n based on the same works.

And finally, the two main flow parameters have to be found: mean flow velocity 𝑢 0 and mean flow height ℎ 0 . In order to guarantee model validity, we should assure that flow is at laminar regime, i.e. Reynolds number not greater than 500. As a value of reference, a first estimative of flow regime can be considered through the maximum flow velocity, which is the surface velocity 𝑢 𝑠 . Values obtained after the first surge and before the wave train (𝑡 < 900 s) suggests that mean surface velocity was about 3.4 m/s. As roll waves crests and troughs originated around mean flow height (Ng & Mei, 1994;Pascal, 1999), one way to estimate ℎ 0 is to calculate the mean value between all waves' crests and troughs. Based on the register from the station, mean flow height should be seen around 0.65 meters. However, the scoured depth of bed should also be considered for further calculations as it doesn't appear directly on the register. Thus, adding he scoured depth (0.25 meters) to the mean value calculated before (0.65 meters), ℎ 0 ≈ 0.90 meters. Flow Reynolds number should not be greater than 160 (for 𝐾 𝑛 = 50 Pas n ). Flow main parameters greatly influence model response to input parameters, so another method is used for sake of comparison.

Since we dispose from measured waves celerity 𝑈, we developed an iterative process to compute ℎ 0 through Equation (3) and 𝑢 0 , through surface velocity 𝑢 𝑠 as relates Equation ( 5):

𝑢 0 = ( 𝑛 + 1 + 𝑛𝐶 * 2𝑛 + 1 ) 𝑢 𝑠 . ( 5
)
Assuming that the momentum coefficient distribution is equal to 1 in the beginning and that we have a mean value for waves propagation velocity equal to 6.21 m/s (mean value calculated from Table 1), variables are solved and rapidly converge to constant values for ℎ 0 = 0.91 meters and 𝑢 0 = 2.85 m/s. Values show good representation of flow properties, in agreement with first assumptions. We consider that the precision is greater to front propagation velocity estimative and then the last value for ℎ 0 and 𝑢 0 are used as input parameter in the model. Fluid and mean properties of flow are summarized on Table 2 as well as dimensionless input parameters to 1D model.

Prediction results

Based on the input parameters set from Table 2, roll waves estimate is performed and results are shown on Figure 4. All waves represented here are the greatest roll waves possible, i.e. the ones with greater amplitude and longest period, given a set of input parameters. This assumption is possible because when flow disturbance is composed by many frequencies, the lowest disturbance frequency will be responsible for permanent waves generated [START_REF] Kranenburg | On the evolution of roll waves[END_REF]. Sensitivity analysis is also performed supposing that we had an uncertainty of 10% in the wave celerity, which would then change ℎ 0 and 𝑢 0 , following the reasoning explained before. All waves characteristics are summarized in Table 3, along with measured waves characteristics from [START_REF] Berti | Debris flow monitoring in the Acquabona watershed on the Dolomites (Italian Alps)[END_REF] for sake of comparison. Not only flow characteristics are important but also rheological properties of fluid play an important role. Based on typical results given by the model from mean values from measures, we can evaluate how waves' properties would change as rheological parameters vary. The assumption that sediment concentration would change along the channel as bed scouring happened is completely plausible. Two other tests cases are simulated considering that fluid could assume a less non-Newtonian behavior (𝑛 = 0.5; 𝜏 𝑐 = 50 Pa) and a more 9 non-Newtonian behavior (𝑛 = 0.2; 𝜏 𝑐 = 500 Pa). These results are presented in comparison with the reference test case for mean flow properties on Figure 5 are summarized in Table 4.

FINAL REMARKS

The first evidence is that the simulation made from the mean value for ℎ 0 have good agreement to the mean amplitude measured, although showing weak correspondence in wave period. The poor agreement between period/wavelength calculated and measured could be a systematic flaw of the model the mathematical model which neglects effects such as surface tension, fluid densification, variations in fluids properties, and others. However, it's also important to note that as we suppose variations on measured properties (thus slightly changing model input parameters) a considerable influence on model results is perceived, reaching higher wave periods and amplitudes. Simulated waves properties can significantly vary depending on rheological properties of fluid. As it can be seen, when non-Newtonian effects are greater, wave amplitude increase and period decreases, although the proportionality to respect of its mean flow does not change significantly. Conversely, when rheological effects are smaller, period is greater, and get closer to measured values, but on the contrary, amplitude decrease.

From an Engineering perspective, results show that the model is capable to estimate maximum roll wave in mudflows events. The results presented show that model is capable to reasonably identify the most important features (amplitude, period and propagation velocity) of possible roll wave formation in mudflow events. However, both flow and rheological properties of fluid significantly affect results on wave properties. This points out that the correct and detailed measurement of properties play an important role to estimative of possible events.

Figure 1 .

 1 Figure 1. Dimensionless schematics on roll waves properties for permanent flow. Main physical properties are remarked. Mean flow height ℎ 0 is a singularity point of roll waves equation.

Figure 4 .

 4 Figure 4. Waves' profiles obtained for Acquabona event from August 17th, 1998. Wave expected for mean flow configuration is simulated and represented by dark line. Two other cases are showed for ± 10% variation on mean waves celerity estimative.

Figure 5 .

 5 Figure 5. Influence of rheological properties (𝑛 and 𝜏 𝑐 ) on waves profiles ℎ(𝑡). Profiles ℎ(𝑡) are showed as function of time 𝑡, scaled to the each mean flow height and to period of waves on mean flow prediction 𝑇 𝑚 , respectively.

Table 1

 1

	. Characteristics of identified surges. Values directly collected from Berti et al (2000)
	measurements.
	Surge	Time [s]	𝑈 [m/s]	Δℎ [m]	𝑇 [s]
	A	895	6.67	1.58	87
	B	982	5.80	1.28	93
	C	1075	7.69	1.07	67
	D	1142	5.88	0.50	33
	E	1175	5.00	1.06	34
	Mean value	6.21	1.10	63

Table 2 .

 2 Fluid and flow mean properties and input parameters for 1D model simulation.

	Fluid properties	Flow properties	Dimensionless parameters
	𝑛	0.3	ℎ 0	0.91 m	𝐹𝑟 𝑚𝑖𝑛 (Eq. 4)	0.21
	𝜏 𝑐 𝐾 𝑛 50-1000 Pas n 200 Pa 𝜌 1430 kg/m 3	𝑢 0 𝜃	2.85 m/s 7°	𝐹𝑟 𝐶 * 𝑅𝑒	0.96 0.13 8-160

Table 3 .

 3 Output data from mathematical/numerical simulation from 1D model to respect of Acquabona event from August 17th, 1998 in comparison with direct results from event measurements.

		Measured	Model prediction
		Berti et al (2000)	-10% Mean + 10%
	𝑈	6.21	5.59	6.22	6.84
	Δℎ	1.10	0.74	1.02	1.27
	𝑇	63	10	17	25

Table 4

 4

	𝑈	6.21	6.83	6.22	6.83
	Δℎ	1.10	0.89	1.02	1.40
	𝑇	63	54	17	16

. Output data from mathematical/numerical simulation from 1D model to respect of Acquabona event from August 17th, 1998 in comparison with direct results from event measurements, when assuming rheological properties variations.

Measured Model prediction

[START_REF] Berti | Debris flow monitoring in the Acquabona watershed on the Dolomites (Italian Alps)[END_REF]

Less non-Newtonian Mean More non-Newtonian 𝑛 = 0.5; 𝜏 𝑐 = 50 Pa 𝑛 = 0.2; 𝜏 𝑐 = 500 Pa

ACKNOWLEDGEMENTS

Authors would like to express their most sincere gratitude to: CNPq (Conselho Nacional de Desenvolvimento Cientifico e Tecnologico do Brasil) for scholarship funding; and Organizing Comitee from CHANGES 2014 for financial support on author attendance in the conference.


Joris Sansen
email: jsansen@labri.fr

Frédéric Lalanne
email: lalanne@labri.fr

David Auber
email: auber@labri.fr

Romain Bourqui
email: bourqui@labri.fr

Adjasankey: Visualization of huge hierarchical weighted and directed graphs

come

Introduction

In various fields, analysing how data flow through a system is central to monitor, understand and predict variations of this system. For example determining in social network, how rumors spread and fade, who starts or stops rumors and if they are transmitted to different communities are important questions. In the same way, in web analytics, users navigation information is valuable and its analysis brings plenty of information on buying habits, fashionable articles, what leads a customer to a purchase or on the contrary, what hinders him/her. Due to improvement in data acquisition techniques, the size and complexity of such data prohibit manual representation. This is particulary true in the era of Big Data as dataset size can exceed Mega/Giga/Tera bytes. Information visualization therefore focuses on the design of automatic, fast and efficient techniques.

A common approach consists in modeling such data with a directed weighted graph where the nodes represent the network entities, the edges and their associated weight represent the information flowing between entities. Various techniques have been proposed during the last decades to visualize such directed weighted graphs. Among these techniques, the most popular are node-link diagrams and adjacency matrices. In node-link diagrams (e.g. [START_REF] Kornaropoulos | DAGView: An Approach for Visualizing Large Graphs[END_REF][START_REF] Selassie | Divided Edge Bundling for Directional Network Data[END_REF]), nodes are visually encoded with a shape (usually a square or a circle) and an edge between two nodes is represented as a poly-line between the corresponding shapes. In adjacency matrices (e.g. [START_REF] Elmqvist | ZAME: Interactive Large-Scale Graph Visualization[END_REF][START_REF] Shen | Path Visualization for Adjacency Matrices[END_REF]), nodes are displayed twice, on the abscissa and on the ordinate, and an edge between two nodes is represented by a shape in the corresponding row and column. These visualization techniques have been experimentally evaluated [START_REF] Ghoniem | On the readability of graphs using node-link and matrix-based representations: a controlled experiment and statistical analysis[END_REF][START_REF] Keller | Matrices or Node-link Diagrams: Which Visual Representation is Better for Visualising Connectivity Models?[END_REF] in order to determine which depiction performs better in different cases. These studies showed that adjacency matrix outperforms node-link diagram when the considered graph becomes large and dense. Other researches defined hybrid encodings that try to take advantage of several methods (e.g. [START_REF] Dinkla | Compressed Adjacency Matrices: Untangling Gene Regulatory Networks[END_REF][START_REF] Hlawatsch | Visual Adjacency Lists for Dynamic Graphs[END_REF][START_REF] Watson | Visualizing very large layered graphs with quilts[END_REF]).

To highlight interesting parts of the network, a hierarchical partition of the graph can be used. Such a hierarchical partition can either be inherent to the data or computed with a clustering algorithm (e.g [START_REF] Vincent D Blondel | Fast unfolding of communities in large networks[END_REF][START_REF] Marinus | Graph clustering by flow simulation[END_REF]). In that case, the main idea is to emphasize each set of the partition as well as the relationships within and between them. Again, both node-link diagram and adjacency matrix have been used to show this hierarchical partition (e.g. [START_REF] Eades | Multilevel visualization of clustered graphs[END_REF][START_REF] Holten | Hierarchical edge bundles: Visualization of adjacency relations in hierarchical data[END_REF][START_REF] Noel | Understanding complex network attack graphs through clustered adjacency matrices[END_REF]).

Even if the techniques mentioned above can scale up to millions of elements, three importants issues still remain: i) clutter ii) lack of efficiency for weighting representation, and iii) slow rendering. The clutter (node-node overlaps, node-edge overlaps and edge-edge crossings) increases with the size and the complexity of the graph, resulting in unreadable visualizations. Moreover, depicting weighting of the graph increases this effect. Also, a large number of elements to be displayed may lead to slow rendering and therefore hinder interactive exploration.

To reduce visual clutter and speed up the rendering, a common approach is to build an abstraction of the original graph. It allows to guide the user in his/her analysis by identifying interesting parts of the network and reduces the focus of his/her study. This is usually done by replacing each subset of a partition of the nodes by a single node (called metanode); and replacing all edges linking two of these subsets by one single edge (called meta-edge), linking the corresponding metanodes. This operation reduces the number of displayed nodes and edges and therefore re-D R A F T duces visual clutter and increases rendering speed. This abstraction is usually called a compound graph, associated with the original graph and the partition. If repeated iteratively, this process allows to build abstractions of higher and higher levels. This use of such method have already been widely studied during the last decade and applied to various depiction techniques [START_REF] Abello | Askgraphview: A large scale graph visualization sys-D R A F T tem[END_REF][START_REF] Archambault | GrouseFlocks: Steerable Exploration of Graph Hierarchy Space[END_REF][START_REF] Balzer | Level-of-detail visualization of clustered graph layouts[END_REF][START_REF] Elmqvist | ZAME: Interactive Large-Scale Graph Visualization[END_REF].

Depicting weighting of relations is possible for most of the above techniques but they suffer from important drawbacks (e.g. lack of efficiency, clutter, readability). A common and effective weighting depiction relates on representation inspired from flow diagrams. Since the first diagrams to modern works [START_REF] Alemasoom | Interactive Visualization of Energy System[END_REF][START_REF] Riehmann | Interactive sankey diagrams[END_REF][START_REF] Vehlow | Visualizing the Evolution of Communities in Dynamic Graphs[END_REF][START_REF] Wongsuphasawat | Exploring Flow, Factors, and Outcomes of Temporal Event Sequences with the Outflow Visualization[END_REF] its visual design kept its efficiency. Such technique emphasizes major trends while depicting quantities and thus allows to locate dominant contribution to the overall flow. However, it faces the problem of scalability and clutter when the number of elements to represent increases. In our knowledge of the visualization domain, there is a lot of publications focusing either on directed graphs, or on weighting depiction, or on hierarchical depiction and usability as a tool but very little on efficient depiction of the three all together.

In this paper, we present Adjasankey, a new visualization technique for representing weighted directed graphs together with a hierarchical partition. To do so, Adjasankey takes advantages of compound graph visualization, adjacency matrix layout and Sankey's diagram that supports the visualization of flows between entities. The aim of Adjasankey is to help the user to answer the following questions: i) Which entities are connected to important flows? ii) Does information mainly flow between or within clusters ? iii) How many flows start from or end to a particular cluster/entity?

The remainder of this paper is structured as follows. In Section 2, we present Adjasankey: design and technical aspects. Then we describe it's implementation in Section 3. We provide some results with a case study based on the analysis of a website in Section 4. Finally, we present conclusions and discuss possible future work.

Adjacency matrix layout for Sankey diagram : Adjasankey

This section presents the design of Adjasankey which uses a node positionning method similar to adjacency matrices combined to Sankey's diagram design for edges. As mentioned in Section 1, we consider the data as a directed and weighted graph together with a hierarchical tree. One can see in Fig. 2(a) an example of a clustered directed graph and the associated hierarchy tree as used as input for our representation.

Positioning and sizing

In Adjasankey, such as in an adjacency matrix, nodes are displayed on rows and columns. As for matrices that represent a directed graph, source nodes are laid out on ordinate and target nodes on abscissa. A well-known matrix issue lies in nodes ordering which can reveal or hide various information. Such node ordering problem is known to be NP-hard and is not the aim of this study. In our method, the height of a node row (resp. width of a node column) is set according to its weight (sum of weights of its incident edges). It helps to identify nodes (and clusters) that are incident to edges with high weights. Since node height (resp. width) is set according to its weight, a node having an out-degree (resp. in-degree) equal to 0 is not displayed on abscissa (resp. on ordinate) which reduces the matrix size. In Fig. 3(a), one can see that the height of E on ordinate is set to 1 while the width of E on abscissa is set to 3.

One can also notice that node B is not drawn on ordinate D R A F T as its out-degree equal to 0. The same observation can be done for nodes A and D on abscissa.

In adjacency matrix, node column and row contains multiple connections to depict all its incident edges. As a node height or width has been set according to its out-and in-degree, it can be split to provide a complete row/column to each of its incident edge and to set their size proportionally to their weight. In Fig. 3(a), the column of node E have been split into two row of width 1 and 2 as node E has two in-coming edges of weight 1 and 2.

Such node and edge positioning allows to bundle together edges having the same source or target without misrepresentation of weights. In order to strengthen this bundling, edges are ordered according first to the in-degree nodes order and then to the out-degree nodes order. In order to emphasize the hierarchical partition of the graph, we designed an algorithm which tries to affect similar colors to nested clusters. The main idea of that algorithm is to assign to each node of the hierarchy tree, a Hrange of value from the HSV colorspace. First an H-range of [0, 360] is assigned to the root node of the hierarchy tree. Then, for each internal node, its H-range is split according to the weight of its sons in the hierarchy tree and the resulting ranges are assigned to them. For instance in Fig. 4, the H-range [0, 360] of the root node is equally split and assign to C1,C2 and E as their total weight are all equal to 4. Finally, the hue value of each node of the hierarchy tree is set to the middle of its H-range. The positioning of nodes and edges allows to preserve for each edge an entire row and column. These row and column can therefore be colored to emphasize edges source and target. Edges color are set according to their source and target colors with a progressive and reversed interpolation. On the source (resp. target) side, edges are colored with target (resp. source) color. Such interpolation results in a cumulative histogram that shows the distribution of weights over the successors (resp. predecessors). For instance, in Fig. 3(b), we can see that in-coming flows of node E on abscissa are emitted from node C and D. Furthermore, the colors emphasize the contribution of node C and D (resp. one third and two third). So as to guarantee the visibility of edges and histogram, a spacing is kept near the entities.

Coloring

Rendering

Due to node placement, adjacent edges are bundled together. In order to emphasize these bundles of edges, they are depicted with a luminance texture that creates a tubular effect. Furthermore, it eases the distinction between edges crossings and edges bends. Indeed crossings identification without such tubular effect can be tedious (see Fig. 3(b)) even when edge borders are drawn. However, we are aware that such design have the side effect to modify color intensity which is a limitation of the technique.

Interaction

The first depiction of the graph is the highest level of abstraction provided by the hierarchical partition (see Fig. 1). We provide together with Adjasankey an interaction tool that allows to modify the diplayed level of detail. Using that interaction tool, the user can request more (resp. less) details on a given cluster. It induces the appearance of descendant nodes and the break up of meta-edges in the corresponding edges, potentially metanodes and metaedges themselves. In addition to that interaction tool, we also provide a classical zoom-and-pan as well as a neighbor identification interaction tools. The latter allows to highlight in the representation the neighborhood of a given node but also the source and target of a given edge/bundle of edges (see Fig. 7).

Implementation

The size of the data that Adjasankey have to handle makes impossible to store and manipulate them on a single computer. In this section, we present some technical aspects of Adjasankey, and explain how we combined a cloud architecture for the pre-processing of the data and a light-weight client that make possible to visualize an abstraction in real-time in web browser. Pre-processing on a cloud architecture To request more or less detail on a particular metanode in reasonable time, the first pre-processing step consists in duplicating the graph and hierarchy tree nodes. For instance in Fig. 5, one can see the result of that duplication step on the graph and hierarchy of Fig. 2. During this process, entities E and C are duplicated since they emit and receive flows while A, B and D are not (see Fig. 2(b)). The hierarchy tree is also updated and clusters containing both sources and targets are duplicated as well. As C1 and C2 contain both emiting and receiving entities (A and B for C1, C for C2), they are duplicated into C1'/C1" and C2'/C2". The second step consists in generating the meta-edge connecting each node of the duplicated hierarchy tree to all other possible neighbors. These two steps are achieved using Apache Spark (a distributed data processing framework) and then stored in an HBase database. A web service finally manages access requests and connection to HBase for data transmission. To improve query response time, our system also keep in memory nodes and metanodes that have already been requested. Using such an architecture, the system can answer in reasonable time when the user requests more (or less) details on the cluster of the hierarchy.

Graphic User Interface Once the duplicated graph and hierarchy tree have been computed, the visualization program is only responsible for the rendering of the coumpound graph. The rendering program is implemented in C++ associated to openGL shader programming. It is then transformed into Javascript and WebGL using the llvmbased project Emscripten [START_REF] Zakai | Emscripten: an LLVM-to-JavaScript compiler[END_REF]. The generated program is integrated within the visualization client written in html.

Case study: website user navigation

This section presents the results that we have obtained on a website user navigation dataset. More precisely, the dataset contains a set of 46, 500, 000 clicked hyperlinks in a Customer To Customer website containing approximatively 50, 000 webpages and 823, 000 hyperlinks. The dataset have been modeled as an oriented (hyperlinks direction) and weighted (number of clics per hyperlink) graph. The hierarchical tree associated to the graph is given by the hierarchical partition in categories (sub-categories and so on) of products for sale. After the duplication pre-process and the computation of the meta-edges, we obtain a set of 96, 300 nodes/metanodes and 1, 291, 300 edges/meta-edges. Fig. 1 shows the representation in Adjasankey of the highest level of abstaction of the data, i.e. the representation of the main categories of products for sale on the website. One can easily notice that the website contains five main categories (or webpages) and that the flows are mainly within category. Among these categories, the Professional equipment category (in blue in Fig. 1) represents approximatively 50% of the visited webpage while the House equipment represents approximatively 30%. Surprisingly, the home page of the website only represents 4% of the visited pages. Two factors can explain that phenomena, first many users enter on the website using search engines that redirect them directlty on particular webpages and second, the number of times the homepage is visited is low compared to the number of visited webpages per visit.

Within the Professional equipment category (see Fig. 6(b) and 7(a)), most of the flows are again within their sub-categories which means that, at this level of organization, the website seems also well-structured. One can also notice that the Agricultural equipment sub-category represents about 50% of the visited webpages of Professional equipment (in yellow in Fig. 7(a)). Other sub-categories either refer to other kinds of equipement or to territorial regions. Two type of behaviours can be observed: i) Users who navigate in territorial regions sub-categories, ii) Users who navigate in semantic sub-categories and then, at lower levels, in territorial regions categories. Such behaviours reveal a distinction between users who wants to promote either territorial proximity or accuracy of their search. As no offer can appear in two categories, it also indicates that users can obtain different results if they choose to follow a territorial strategy or not.

Considering the House equipment category (see Fig. 6(c) and 7(b)), one can observe two major subcategories (Gardening and Shoes), one medium category (Baby clothes) and two small categories (Luggage and Tools). While major part of the flows are within categories, two interesting behaviors can be noticed: i) almost all of the flows connected to Tools are emitted from and toward Gardening, ii) almost all of the flows connected to Luggage are emitted from and toward Shoes. the case of Tools and Gardening, when approximatively 90% of the flows from Tools are emitted toward Gardening.

We assume that this is the result of a miss-classification of the offers and that the Gardening category should be included in the Tools one. Indeed many offers of the Gardening category refers to gardening tools offers. The same reasoning applies to the flows between Luggage and Shoes flows case.

Conclusion and future work

We have presented Adjasankey, a novel visualization technique designed for weighted directed graphs associated to an hierarchical partition. Our visualization method combines Sankey's diagram visual design with an adjacency matrix layout and thus emphasizes edge weighting while optimizing overall readability. To scale to huge dataset, Adjasankey supports the multi-scale exploration of various levels of abstraction. We have implemented our solution in a Big Data infrastructure with cloud computing for data pre-processing and a light-weight client for visualization. This system allows a time-responsive interaction and visualization of various abstractions of large network. A trial of this system was made with a case study on website user navigation and was positively received by domain experts.

As future work, we plan to extend Adjasankey to path or sequence visualization as a combination of our depiction and of space-filling curves. We also plan to perform an experimental evaluation for comparing classical adjacency matrix diagram and Adjasankey for the depiction of directed and weighted graphs.

Figure 1 :

 1 Figure 1: Visualization of the highest level of abstraction with Adjasankey obtained on a dataset of 46, 500, 000 clicked hyperlinks in a C2C website.

Figure 2 :

 2 Figure 2: Example of directed weighted graph (a) and associated hierarchical partition (b).

Figure 3 :

 3 Figure 3: (a) Nodes and edges positionning on matrix layout according to the input graph Fig. 2; (b) resulting visualization with node and edge coloring.

Figure 4 :

 4 Figure 4: H-range of HSV colorspace distribution for node coloring.

Figure 5 :

 5 Figure 5: Nodes having both ingoing and outgoing edges in the graph of Fig. 2 are duplicated in both the graph (a) and the hierarchy (b).

Figure 6 :

 6 Figure 6: Focus on the home page (a) and on the Professional equipment (b) and House equipment (c) categories.

Figure 7 :

 7 Figure 7: Different levels of detail and zoom levels. (a) Expansion of Professional equipment metanode; (b) Expansion of and zoom on the House equipment metanode; (c) Flows emitted from the Gardening category ; (d) Details of Tools (in green) and Gardening (in yellow) flows: most of the flow emitted from Tools is emitted toward Gardening.

Acknowledgement

This work has been carried out as part of SpeedData project supported by the French Investissement dAvenir Program (Big Data -Cloud Computing topic) managed by Direction Gnrale des Entreprises (DGE) and BPIFrance.