Sophie Cerf
email: sophie.cerf@gipsa-lab.fr

Bogdan Robu
email: bogdan.robu@gipsa-lab.fr

Nicolas Marchand
email: nicolas.marchand@gipsa-lab.fr

Antoine Boutet
email: antoine.boutet@insa-lyon.fr

Vincent Primault
email: vincent.primault@insa-lyon.fr

Sonia Ben Mokhtar
email: sonia.benmokhtar@insa-lyon.fr

Sara Bouchenak
email: sara.bouchenak@insa-lyon.fr

Toward an Easy Configuration of Location Privacy Protection Mechanisms

Keywords: CCS Concepts, Security and privacy → Pseudonymity, anonymity and untraceability, Privacy-preserving protocols, File system security

The widespread adoption of Location-Based Services (LBSs) has come with controversy about privacy. While leveraging location information leads to improving services through geo-contextualization, it rises privacy concerns as new knowledge can be inferred from location records, such as home/work places, habits or religious beliefs. To overcome this problem, several Location Privacy Protection Mechanisms (LPPMs) have been proposed in the literature these last years. However, every mechanism comes with its own configuration parameters that directly impact the privacy guarantees and the resulting utility of protected data. In this context, it can be difficult for a non-expert system designer to choose appropriate configuration parameters to use according to the expected privacy and utility.

In this paper, we present a framework enabling the easy configuration of LPPMs. To achieve that, our framework performs an offline, in-depth automated analysis of LPPMs to provide the formal relationship between their configuration parameters and both privacy and the utility metrics. This framework is modular: by using different metrics, a system designer is able to fine-tune her LPPM according to her expected privacy and utility guarantees (i.e., the guarantee itself and the level of this guarantee). To illustrate the capability of our framework, we analyse Geo-Indistinguishability (a well known differentially private LPPM) and we provide the formal relationship between its configuration parameter and two privacy and utility metrics.

INTRODUCTION

Location Based Services (LBSs) such as navigation or recommendation applications have been widely adopted by peo-Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s). ACM ISBN 978-1-4503-4666-5/16/12. DOI: http://dx.doi.org/10.1145/3007592.3007599 ple using mobile devices. Although location-aware systems have greatly improved the quality of many services by introducing geo-contextualization, such systems rise important concerns about privacy.

Middleware Posters and

Indeed, a collection of mobility traces (a set of timestamped locations reflecting the user's moving activity) can reveal many sensitive information about its user such as home and work places, favorite restaurants, religious leaning, or the individuals she met. To overcome this privacy issue, many Location Privacy Protection Mechanisms (LPPMs) have been proposed in the last decade. However, the effectiveness of these mechanisms usually rely on the tuning of a set of configuration parameters, often with a large range of possible values. In addition, this tuning depends on the expected privacy guarantee and utility of the resulting protected data, which form two conflicting assessment dimensions. Consequently, this parametrization task can be difficult for a non-expert system designer.

Few works have been done to automatically assess LPPMs in terms of privacy and utility as well as to assist a system designer in their configuration. Most of the existing LPPMs are not automated and focus on privacy guarantees [START_REF] Agir | User-side adaptive protection of location privacy in participatory sensing[END_REF][START_REF] Chatzikokolakis | Constructing elastic distinguishability metrics for location privacy[END_REF]. To the best of our knowledge, ALP [START_REF] Primault | Adaptive location privacy with alp[END_REF] is the only solution enabeling to address this problem. Specifically, ALP uses a greedy solution to possibly make the configuration parameters converge to values which aim to maximize or minimize given privacy or utility metrics. As far as we know, no solution provides the formal relationship between configuration parameters and privacy and utility metrics as we propose.

In this paper, we present a framework which aims to help in the fine-tuning of LPPMs according to a set of expectations in term of privacy and utility defined by a system designer. To achieve that, our solution adopts an automated approach which performs an in-depth analysis of LPPMs, and provides the formal relationship between their configuration parameters and both privacy and utility metrics. This framework is modular: by using different metrics it is possible to adapt the provided model to specific privacy and utility guarantees.

Through a real example, we illustrate the capability of our framework to configure Geo-Indistinguishability [START_REF] Andrés | Geo-indistinguishability: Differential Privacy for Location-based Systems[END_REF], a well known LPPM based on differential privacy, according to privacy and utility metrics.

PROBLEM STATEMENT

To illustrate the problem related to the configuration of a LPPM and the capability of our framework to solve it, we consider the case of using Geo-indistinguishability [START_REF] Andrés | Geo-indistinguishability: Differential Privacy for Location-based Systems[END_REF] (GEO-ǫ (in meters -1 , log scale) I for short) to protect a whole dataset containing mobility traces of taxi drivers around San Franciso, with objectives in terms of privacy and utility. GEO-I [START_REF] Andrés | Geo-indistinguishability: Differential Privacy for Location-based Systems[END_REF] follows the differential privacy model and is specifically designed to protect mobility data. To protect raw data, this LPPM adds random noise drawn from a Laplace distribution to the actual user location. GEO-I is parametrized by an parameter (expressed in meters -1 ), quantifying the amount of noise to add (the lower the , the higher the noise).

As privacy objective, we consider the retrieval in the protected data of at most 10% of the Points of interest (POIs) of users (meaningful locations where a user made a significant stop). This objective is evaluated through a privacy metric which quantifies the proportion of actual POIs retrieved from the protected data for each user.

As utility objective, we consider maintaining a similar location precision at the scale of a city block. More precisely, the difference between the area coverage of users in the actual mobility traces and their protected counterpart is expected to remain about the size of a city block block and no less accurate. This objective is quantified through a utility metric.

The challenge here for a system designer is to correctly configure the of GEO-I to ensure that the expected objectives are respected. Figure 1 shows the level of privacy and utility according to different values of (the vertical lines shows zones where metrics are not saturated). Results show that the privacy metric rapidly changes from 0 to 0.4 with an from 0.007 to 0.08, respectively (Figure 1a). The utility metric, in turn, evolves more slowly on a larger range from 0.2 to 1 with an from 10 -4 to 1, respectively (Figure 1b). Consequently, to ensure the expected privacy and utility objectives, the system designer has to choose carefully its . In our case, a = 0.01 minimizes the privacy while maximizing the utility. In other words, with l = 0.01 we ensure to a user that no more than 10% of her POIs can be retrieved while ensuring that 80% of her requests will concern the city block where she is.

OVERVIEW OF THE FRAMEWORK

Our framework proposes an automated approach to correctly configure LPPMs' parameters for given objectives in terms of privacy and utility, using a mathematical model of the LPPM behavior. Our solution proceeds in three automated steps:

1. First, the system needs to be defined: (1) the objective metrics for privacy (P r) and utility (U t), (2) the LPPM configuration parameters pi and their range of values, and (3) the properties of the dataset di that are likely to influence privacy and utility metrics (i.e., reflecting impactfully characteristics of users such as the uniqueness). All these properties pi, di are soundly chosen using a principal component analysis. For instance, in our illustration based on GEO-I, P r measures the proportion of POIs that can be retrieved, U t quantifies the difference in the service area coverage, only the parameter configures the LPPM, and no dataset properties is considered. 2. Then comes the modeling phase: experiments are automatically run where parameters pi and di vary in turn while evaluation metrics are measured. Based on this data, a mathematical relationship between privacy and utility metrics, configuration parameters, and dataset properties is computed as an invertible function: P r U t = f (p1, . . . , pn, d1, . . . , dm)

3. Finally, the LPPM configuration (i.e. the value of pi) is computed by inverting the f function, using the specified privacy and utility objectives and the evaluation of dataset properties di.

For instance, in our illustration based on GEO-I, P r and U t only depend on parameter . If we focus only on the interval where impacts the privacy and utility metrics (i.e. between the vertical lines in Figure 1) we can approximate the experimental curves with the following linear equations:

log( ) = P r -a b = U t -α β ,
with a = 0.84, b = 0.17, α = 1.21 and β = 0.09.

(

Hence, to ensure that her privacy and utility objectives are met, a system designer can easily leverage this model to determine the appropriate value of . Specifically, to guarantee 10% privacy, configuring = 0.01 ensures 80% utility.

CONCLUSION

In this paper we proposed a framework that enables easy and automated configuration of LPPMs according to privacy and utility objectives. Our preliminary results on the LPPM GEO-I show promising results as for the control of configuration parameter. Our future work will focus in testing other LPPMs and datasets, we also plan to extend our framework with more metrics and parameters.

  Demos '16 December 12-16 2016, Trento, Italy c 2016 Copyright held by the owner/author(s).

Figure 1 :

 1 Figure 1: GEO-I configuration parameter according to privacy and utility metrics


Francky Saint-Ange

Bruno Savoye

Laurent Michon

Patrick Bachèlery

Christine Deplus

Béatrice De Voogd

Jérôme Dyment

Eliane Le Drezen

Michel Voisset

Anne Le Friant

Patrick Bachelery

Georges Boudon

A volcaniclastic deep-sea fan off La Réunion Island (Indian Ocean): Gradualism versus catastrophism

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

INTRODUCTION

La Reunion Island is an active basaltic shield volcano, 220 km in diameter, rising ~7000 m from the Indian Ocean fl oor. The island is constructed on Paleocene oceanic crust and is considered to be the surface expression of a hotspot [START_REF] Duncan | The volcanic record of the Réunion hotspot in Duncan[END_REF]. The submarine fl anks of La Réunion are mainly debris avalanche deposits [START_REF] Oehler | Growth and collapse of the Reunion Island volcanoes[END_REF]. Little previous attention has been paid to the volcaniclastic sedimentation on the slope apron and the adjacent abyssal plain.

Studies on volcaniclastic sedimentation during the 1980s and 1990s were focused on shortterm, high-amplitude events, such as debris avalanches and pyroclastic fl ows [START_REF] Fisher | Marginal basin geology; Volcanic and associated sedimentary and tectonic processes in modern and ancient marginal basins[END_REF][START_REF] Moore | Prodigious submarine landslides on the Hawaiian Ridge[END_REF][START_REF] Schmincke | Volcanic evolution of Gran Canaria reconstructed from apron sediments: Synthesis of VICAP project drilling[END_REF][START_REF] Masson | Slope failures on the fl anks of the western Canary Islands[END_REF]. Very high-resolution studies carried on the submarine fl anks of La Palma and El Hierro in the Canary Archipelago [START_REF] Wynn | Turbidity current sediment waves on the submarine slopes of the western Canary Islands[END_REF], as well as on Stromboli volcano in the Mediterranean Sea [START_REF] Kidd | Sedimentary processes in the Stromboli Canyon and Marsili Basin, SE Tyrrhenian Sea: Results from side-scan sonar surveys[END_REF][START_REF] Casalbore | Morpho-sedimentary characteristics of the volcaniclastic apron around Stromboli volcano[END_REF], have shown evidence of sedimentary processes that are not linked with giant volcanic landslides. Recent advances on source-to-sink processes in volcanic setting [START_REF] Manville | Source to sink: A review of three decades of progress in the understanding of volcaniclastic processes, deposits, and hazards[END_REF] offer a more dynamic focus on volcaniclastic sedimentation by putting together the nature and dynamics of the volcanic products and the environmental setting.

To the authors' knowledge, no previous studies have examined the case of a torrential river supplying a Quaternary volcaniclastic turbidite system. We present here the case of the Cilaos deep-sea fan, the sediment source of which is a torrential river draining a dormant volcanic center on La Réunion. Our objective is to show that long-term erosion processes and high sediment supply in a volcanic-island setting can also form deep-sea fans similar to siliciclastic ones.

We combined acoustic backscatter images and swath bathymetry data (Simrad EM12D and EM120), 3.5 kHz echosounder profi les, and piston cores collected during cruises FOR-EVER (N/O L'Atalante) and ERODER 1 (N/O Beautemps-Beaupré) in 2006. Laser-diffraction grain-size analysis and X-ray imaging were performed on the cores. Due to the low content of bioclasts in our core sampling, radiocarbon dates were not conclusive.

SOURCE OF THE VOLCANICLASTIC TURBIDITE SYSTEM

The Cilaos deep-sea fan is connected to a small drainage basin (360 km 2 ) up to 3000 m elevation via the Rivière Saint-Etienne (Figs. 1 and2). The river has built an alluvial fan-delta no older than 350 ka, with a volume of ~25 km 3 (Saint-Ange, 2009). The river bed is composed of poorly sorted, coarse-grained sediments, ranging from sands to boulders (Fig. 3). The present river outlet is directly connected to a submarine valley. There is no shelf in front of † Deceased 21 August 2008. The sediments are transported by annual fl ash fl oods linked to the passage of tropical hurricanes and storms that increase the river discharge from 4.6 m 3 s -1 up to 1500 m 3 s -1 , with fl ood velocity often higher than 6 m s -1 (SOGREAH, 1998). The river discharge doubles for a 10 yr return period, and centennial fl oods are estimated to reach 5000 m 3 s -1 [START_REF] Sogreah | Société Grenobloise d'Etudes et d'Applications Hydrauliques[END_REF]. This type of fl ow leads to a mean solid load to the ocean of ~0.5 × 10 -3 km 3 in a few days, and for a centennial fl ood the solid load is estimated around 1-2 × 10 -3 km 3 [START_REF] Sogreah | Société Grenobloise d'Etudes et d'Applications Hydrauliques[END_REF]. Severity of fl ooding is illustrated by the 2007 fl oods induced by Hurricane Gamède, during which the fl ow reached 2000 m 3 s -1 (CGPC, 2007) and destroyed a major highway bridge (Fig. 3).

ABSTRACT

MORPHOLOGY AND NATURE OF THE CILAOS TURBIDITE SYSTEM Feeder Valleys

The Cilaos turbidite system is more than 280 km long and 100 km wide (Fig. 1). It starts at the coast with a set of 2-4-km-wide valleys (mean bank elevation ~150 m), which merge into a single wider valley at about -3000 m of water depth (Figs. 1 and4). These valleys are 70 km long and cut into debris avalanche deposits older than 500 ka [START_REF] Oehler | Growth and collapse of the Reunion Island volcanoes[END_REF]. Close to the coast, the valley fl oor gradient is ~8°, but it decreases quickly to reach 1° at the base of the volcanic edifi ce (Fig. 4; Fig. DR1). The high backscatter signal from valley fl oors is associated with very coarse sands and gravels, as revealed by piston core 03 (see Fig. 1 for location).

Morphology of the Fan

The upper part of the fan is characterized by fi elds of sediment waves (≥10 m high and 1 km of wavelength), the crest orientations of which are perpendicular to the valley axis. The sediment waves developed on slopes ranging between 2° and 0.4° (Fig. 4). They fi rst start to grow on the western part of the lower valley (Figs. 1 and4), where the height of the fl ank of the valley decreases signifi cantly.

The main deep-sea fan is developed at ~4000-4500 m below sea level (mbsl) on a complex abyssal plain morphology shaped by northeast-southwest-trending volcanic ridges (Fig. 1). The volcanic ridges divide the turbidite system into three parts (western, central, and eastern), presenting contrasting backscatter acoustic facies.

The eastern and central parts of the fan show few distinct architectural elements, mostly short channels observed in their upper areas only. In contrast, elements such as sediment waves and channels are well developed in the western part of the fan. Here, the upper fan shows a braided channel system (Fig. 1; Fig. DR2), which ends quickly or merges into a single channel that extends down to the lower part of the turbidite system (Figs. 1 and5). The 3.5 kHz data reveal it is as big as the siliciclastic fan of the Var (16,000 km 2 ) in the Mediterranean Sea, although its drainage basin is seven times smaller. The size of the fan in relation to its drainage system raises the question of rate of sediment supply.

Most hotspot-related volcanic islands in tropical or subtropical areas present three main similarities: (1) large production of volcaniclastic detritus linked to recurrent giant landslides [START_REF] Moore | Prodigious submarine landslides on the Hawaiian Ridge[END_REF][START_REF] Masson | Slope failures on the fl anks of the western Canary Islands[END_REF][START_REF] Oehler | Growth and collapse of the Reunion Island volcanoes[END_REF], (2) high magmatic production rates that frequently bring new material to erosion [START_REF] Duncan | The volcanic record of the Réunion hotspot in Duncan[END_REF][START_REF] Schmincke | Volcanic evolution of Gran Canaria reconstructed from apron sediments: Synthesis of VICAP project drilling[END_REF], and (3) among the highest worldwide erosion rates [START_REF] Louvat | Present denudation rates on the island of Réunion determined by river geochemistry: Basalt weathering and mass budget between chemical and mechanical erosions[END_REF]. The maximum denudation rates of some volcanic islands (Gran Canaria: 1.8 m k.y. -1 [START_REF] Menéndez | Fluvial dissection, isostatic uplift, and geomorphological evolution of volcanic islands[END_REF]; Kauai: 4 m k.y. -1 [Gayer et al., 2008]; La Réunion: 3.4 m k.y. -1 [START_REF] Louvat | Present denudation rates on the island of Réunion determined by river geochemistry: Basalt weathering and mass budget between chemical and mechanical erosions[END_REF]) are close to those estimated for active orogens (Taiwan: 3-6 m k.y. -1 [START_REF] Dadson | Links between erosion, runoff variability and seismicity in the Taiwan orogen[END_REF] Himalaya: 2-5 m k.y. -1 [START_REF] Burbank | Decoupling of erosion and precipitation in the Himalayas[END_REF]). In the case of La Réunion, we propose that the continuous magmatic activity, which lasted at least 2 m.y. between 2.1 and 0.012 Ma in the Cilaos area [START_REF] Mcdougall | The geochronology and evolution of the young volcanic island of Réunion, Indian Ocean[END_REF][START_REF] Deniel | New 230 Th-238 U and 14 C age determinations from Piton des Neiges volcano, Reunion-A revised chronology for the Differentiated Series[END_REF], combined with the high erosion rate were able to continuously provide enough material, thus explaining the high sediment supply with respect to the small size of the drainage basin.

This explanation highlights the fact that longterm erosion processes account for a signifi cant part of the volcaniclastic sediment output, and in the case of La Réunion, for the most signifi cant part of the late Quaternary.

Sediment Transfer to the Basin

The drainage basin and the alluvial fan combined represent an area of 445 km 2 . This value implies that the sediment storage capacity onshore represents only ~3% of the entire system. In view of the low storage capacity (Figs. 1 and2), the high erosion rate, and the river discharge, the sediments cannot stay a long time onshore and are most likely transferred to the basin during fl ood events related to storm rainfall. Because the river is permanently connected with the submarine valley, there is no problem of accommodation due to the sea-level drops or rises, as they do not signifi cantly impact the sediment transfer to the deep sea. All these parameters, together with the island's steep submarine upper slopes (Fig. 4) and the well-developed valleys cut into the debris avalanche deposits, provide an effi cient transfer mechanism for volcaniclastic sediments from subaerial regions to the abyssal plain. The recurrent rhythmic sequences observed in cores and the obvious channel systems indicate that turbidity currents arrive frequently in the upper part of the fan. The resulting sediment deposits do not show the characteristic signatures of turbidites related to giant volcanic landslides as identifi ed for other volcanic islands such as Canary Islands and Hawaiian Islands. In those areas, sediment record shows stacked subunits within a single thick turbidite bed (>1 m), indicating multiple stages of failure [START_REF] Garcia | Spatial variability of erosion rates inferred from the frequency distribution of cosmogenic 3 He in olivines from Hawaiian river sediments[END_REF][START_REF] Wynn | Canary Islands landslides and tsunami generation: Can we use turbidite deposits to interpret landslide processes?[END_REF]. This difference implies another mechanism for sediment transport than the catastrophic giant slopes failures.

In the case of La Réunion, the magnitude of river fl oods suggests that hyperpycnal fl ows contributed signifi cantly to the development and the feeding of the Cilaos turbidite system. This hypothesis is also supported by the braided plain observed in the upper western part of the fan, the growth and development of which required a sustained fl ow, with high discharge and predominance of bed-load material [START_REF] Hesse | Sandy submarine braid plains: Potential deep-water reservoirs[END_REF]. The hypothesis of hyperpycnal fl ow feeding a turbidite basin has also been proposed in Iceland, where turbidites related to Jökulhlaup events are considered to have been a signifi cant source of sediment to the Iceland Basin over the past 3 m.y. [START_REF] Elliott | Infl uence of sediment drift accumulation on the passage of gravity-driven sediment fl ows in the Iceland Basin[END_REF]. In addition, considering the high bed-load river type of La Rivière Saint-Etienne, small failures can occur immediately seaward of the main river mouth and secondary channels, as described for Squamish delta in British Columbia and similar settings [START_REF] Piper | Processes that initiate turbidity currents and their infl uence on turbidites: A marine geology perspective[END_REF].

CONSEQUENCES FOR VOLCANICLASTIC SEDIMENTATION MODELS

The existence of hyperpycnal fl ows feeding a volcaniclastic deep-sea fan introduces the idea of a gradual evolution of volcaniclastic sedimentation similar to the siliciclastic model. This model contrasts with the idea of catastrophic sedimentation related to fl ank collapses and pyroclastic fl ows. Following classifi cations made for turbidite systems by [START_REF] Richards | Submarine-fan systems: I. Characterisation and stratigraphic prediction[END_REF] and [START_REF] Piper | Sandy fans-From Amazon to Hueneme and beyond[END_REF], Cilaos is classifi ed as a sand-rich system. Commonly, sand-rich systems are found in active tectonic settings, such as along the California margin or in the Gulf of Corinth [START_REF] Piper | Sandy fans-From Amazon to Hueneme and beyond[END_REF]. Space available for development of fans in these settings is reduced, so they tend to be thicker rather than wide [START_REF] Richards | Submarine-fan systems: I. Characterisation and stratigraphic prediction[END_REF][START_REF] Piper | Sandy fans-From Amazon to Hueneme and beyond[END_REF]. In these classifi cations, the Cilaos deep-sea fan represents an alternative model and shows a simple and predictable pattern, as (1) it formed in an open-ocean setting isolated from siliciclastic source, (2) it is a wide volcaniclastic sandy system similar to the siliciclastic ones, and (3) it shows lateral and vertical variations of facies and architecture.

This Quaternary model echoes the thick and well-developed volcaniclastic turbidite systems described in ancient geological series like in the Kalgoorlie Sequence in Australia [START_REF] Krapez | Late Archaean deepmarine volcaniclastic sedimentation in an arcrelated basin: The Kalgoorlie Sequence of the Eastern Goldfi elds superterrane, Yilgarn craton[END_REF], or in the Lower Mesozoic Mineral King caldera complex in Sierra Nevada that the channel morphology becomes deeper and more confi ned downstream (Fig. 5). In the distal part of the turbidite system, the channel reaches its maximum size, where it is ~5 km wide and ~30 m deep (Fig. 5, profi le 20).

Nature of Sediments

The cores taken on overbank areas consist mainly of volcaniclastic detritus ranging from sands to clays, with silts and sands representing between 78% and 88% of the sediment. Sandy layers consist of dark basaltic fragments, brown glass, olivine, feldspar, pyroxene, and in some cases zeolites. True pelagic layers are rare, and hemipelagic layers are poor in bioclasts. Pumices are rare, but some are observed inside hemipelagic layers. The cores show classical turbidite sequences with traction deposits, such as ripple marks, laminated sand in places alternating with laminated bioclasts, and thick sand deposits, respectively, Tc, Tb, and Ta of the Bouma sequence (Fig. 4; Fig. DR3). These turbidites are rhythmic, with closely spaced sequences at the head of the fan and more interbedded hemipelagic sediment down the fan (Fig. 4; Fig. DR3).

SEDIMENTARY PROCESSES INVOLVED IN THE GROWTH OF THE TURBIDITE SYSTEM

Sediment Supply

The Cilaos deep-sea fan covers an area of ~15,000 km 2 and has a maximum thickness of 120 m [START_REF] Saint-Ange | La Sédimentation Volcanoclastique en Contexte de Point Chaud (Île de La Réunion[END_REF]. By comparison, (California;[START_REF] Busby-Spera | A sand-rich submarine fan in the Lower Mesozoic Mineral King Caldera complex, Sierra Nevada, California[END_REF]. In those areas, sediments were essentially supplied by reworking of volcaniclastic deposits. Analogue structures to those described in this article are observed, such as channelized and unchannelized deposits, cross-bedded structures, or braided plain channels [START_REF] Busby-Spera | A sand-rich submarine fan in the Lower Mesozoic Mineral King Caldera complex, Sierra Nevada, California[END_REF][START_REF] Krapez | Late Archaean deepmarine volcaniclastic sedimentation in an arcrelated basin: The Kalgoorlie Sequence of the Eastern Goldfi elds superterrane, Yilgarn craton[END_REF].

Quaternary deep-sea volcaniclastic sediments are widespread in many basins with high economic interest [START_REF] Fisher | Marginal basin geology; Volcanic and associated sedimentary and tectonic processes in modern and ancient marginal basins[END_REF][START_REF] Alibés | Determination of sediment volumes, accumulation rates and turbidite emplacement frequencies on the Madeira Abyssal Plain (NE Atlantic): A correlation between seismic and borehole data[END_REF][START_REF] Elliott | Infl uence of sediment drift accumulation on the passage of gravity-driven sediment fl ows in the Iceland Basin[END_REF], and they represent a complexity for the mechanism of siliciclastic sedimentation. In those areas, volcaniclastic sedimentation is mainly seen as linked with episodic and ephemeral massive processes [START_REF] Manville | Source to sink: A review of three decades of progress in the understanding of volcaniclastic processes, deposits, and hazards[END_REF]. Nevertheless the existence of the Cilaos turbidite system underlines the lack of understanding of Quaternary volcaniclastic sedimentation linked to long-term surface processes.

Flash fl oods, especially those that destroy bridges, are considered to be catastrophic events. Nevertheless, at the Quaternary scale and considering the phenomena of volcanic fl ank collapses, they are consistent events because they are linked with climate and long-term physical processes. Because the Cilaos fan is the result of fl ood supply of sediment by rivers, it appears to be a good example of evolution from a catastrophic sedimentary system model characteristic of volcanic environments to a gradual sedimentary system model characteristic of siliciclastic environments. It also provides a good opportunity to study in detail the different facies and architecture related to subaerial fl ood processes.

  A new geophysical data set off La Réunion Island (western Indian Ocean) reveals a large volcaniclastic submarine fan developing in an open-ocean setting. The fan is connected to a torrential river that fl oods during tropical cyclones. Sediment storage at the coast is limited, suggesting that the sediments are carried directly to the basin. The fan morphology and turbidites in cores lead us to classify it as a sand-rich system mainly fed by hyperpycnal fl ows. In the ancient geological record, there are many examples of thick volcaniclastic successions, but studies of modern analogues have emphasized mechanisms such as debris avalanches or direct pyroclastic fl ow into the sea. Because the Cilaos deep-sea fan is isolated from any continental source, it provides information on architecture and noncatastrophic processes in a volcaniclastic deep-sea fan.

Figure 1 .

 1 Figure 1. Acoustic backscatter images of southwestern part of La Réunion gathered during recent cruises FOREVER and ERODER. Red dots indicate location of cores, and blue lines indicate location of 3.5 kHz echosounder profi les shown in Figure 5. A: Close-up of La Rivière Saint-Etienne valleys. B: Location of Figure 2. Abbreviations: dad-debris avalanche deposits; dsf-deep-sea fan; vr-volcanic ridges.

Figure 2 .

 2 Figure 2. Map of alluvial fan of Rivière Saint-Etienne and bathymetry of upper part of slopes (contour interval: 10 m; artifact in drawing of contour lines highlights limit between swath and conventional bathymetry data sets). Onshore shaded topography was produced using a digital terrain model provided by Institut Géographique National (IGN), France. Conventional bathymetry data were provided by Service Hydrographique et Océanographique de la Marine (SHOM), France.

Figure 3 .

 3 Figure 3. Photo of Rivière Saint-Etienne after Hurricane Gamède (2007). Note size of boulders (cars for scale).

Figure 4 .

 4 Figure 4. Topographic profi le along Cilaos turbidite system (vertical exaggeration: 7×) and main sedimentary sequences observed along the system: a-gravel beds; b-ripple marks and laminated sands; c-laminated sands; d-sandy beds.

Figure 5 .

 5 Figure 5. Part of four 3.5 kHz echosounder profi les spanning Cilaos fan (see location on Fig. 1). Numbers indicate width (km) and height (m) of channels. Abbreviations: c-channel; cls-channel-levee system; dsc-discontinuous subsurface refl ectors; es-erosion structures; pc-paleochannel; sw-sediment waves; lba-low-backscatter area.

ACKNOWLEDGMENTS

We thank the captains and crews of the N/O L'Atalante and N/O Beautemps-Beaupré for their excellent work during cruises FOREVER and ERODER1, respectively. The technical and fi nancial support of L'Institut Francais de Recherche pour l'Exploitation de la Mer, Service Hydrographique et Océanographique de la Marine, and Centre National de la Recherche Scientifi que-Institut National des Sciences de l'Univers is also gratefully acknowledged. Thanks are due to Hervé Douris for the photos of la Rivière Saint-Étienne. An early draft of the manuscript was signifi cantly improved based on comments from D.J.W. Piper, P.J. Talling, J.-L. Schneider, R.B. Wynn, and D.G. Masson. We are also grateful to the two anonymous reviewers for their constructive reviews. This is Institut de Physique du Globe de Paris (IPGP) contribution 3025.


Caterina Calgaro
email: caterina.calgaro@math.univ-lille1.fr

Meriem Ezzoug
email: meriemezzoug@yahoo.fr

Ezzeddine Zahrouni
email: ezzeddine.zahrouni@fsm.rnu.tn

On the global existence of weak solution for a multiphasic incompressible fluid model with Korteweg stress

Keywords: Mixture theory, Kazhikhov-Smagulov model, Korteweg model, global existence result

In this paper, we study a multiphasic incompressible fluid model, called the Kazhikhov-Smagulov model, with a particular viscous stress tensor, introduced by Bresch and co-authors, and a specific diffusive interface term introduced for the first time by Korteweg in 1901. We prove that this model is globally well posed in a 3D bounded domain.

Introduction

In this paper, we are interested in the study of a system of PDEs describing the evolution of mixture flows. Let Ω be a bounded open set in R 3 with boundary Γ that is regular enough and let n be the outward unit normal on the boundary Γ. We denote by [0, T ] the time interval, for T > 0. We consider the following system of PDEs derived from the compressible Navier-Stokes equations with a Korteweg stress tensor. The closure of the system is given by a Fick's law which relates the velocity to derivatives of the density. This system, that we call the Kazhikhov-Smagulov-Korteweg (KSK) model, reads

         ρ ∂ t u + (u • ∇)u -λ(∇ρ • ∇)u - λ 2 div ρ∇u -ρ∇u T + ∇P = ρg -κ∆ρ∇ρ, ∂ t ρ + div ρu = λ∆ρ, div u = 0. (1)
With the notations Q T = (0, T ) × Ω and Σ = (0, T ) × Γ, the unknowns for this model are ρ : Q T → R the density of the fluid, u : Q T → R 3 the mean volume velocity field and P : Q T → R the pressure of the fluid (a modified pressure). Moreover, g stands for the gravity acceleration (but it can include further external forces) and the parameters λ > 0 and κ > 0 represent mass diffusion coefficient and Korteweg's constant, respectively. The KSK model ( 1) is completed with the following boundary and initial conditions u(t, x) = 0, ∂ρ ∂n (t, x) = 0, (t, x) ∈ Σ,

u(0, x) = u 0 (x), ρ(0, x) = ρ 0 (x), x ∈ Ω,

with the compatibility condition div u 0 = 0, where ρ 0 : Ω → R and u 0 : Ω → R 3 are given functions. Throughout this work, we assume the hypothesis

0 < m ≤ ρ 0 (x) ≤ M < +∞, x ∈ Ω. (4)
Let us mention some known results about the Kazhikhov-Smagulov model without the Korteweg stress tensor. Let us also mention other results about the incompressible Navier-Stokes equations with a Korteweg stress tensor, coupled with a convection or a convection-diffusion equation for the density of the mixture. First, we take κ = 0. In [START_REF] Bresch | Effects of density dependent viscosities on multiphasic incompressible fluid models[END_REF], Bresch et al. consider the problem (1) (without Korteweg term), ( 2) and [START_REF] Da Veiga | Diffusion on viscous fluids. Existence and asymptotic properties of solutions[END_REF]. They prove the global existence of weak solution without assuming small data and without any assumption on the diffusivity λ. In [START_REF] Antontsev | Boundary value problems in mechanics of nonhomogeneous fluids[END_REF][START_REF] Da Veiga | Diffusion on viscous fluids. Existence and asymptotic properties of solutions[END_REF], the equation (1) 1 is replaced by

ρ ∂ t u+(u•∇)u -λ(∇ρ•∇)u-λ(u•∇)∇ρ-µ ∆u+∇P + λ 2 ρ ∆ρ∇ρ+ ∇ρ•∇ ∇ρ- |∇ρ| 2 ρ ∇ρ = ρ g. (5)
In [START_REF] Antontsev | Boundary value problems in mechanics of nonhomogeneous fluids[END_REF], Antontsev et al. study the problem (5) -(1) 2 -(1) 3 , ( 2) and (3) without O(λ 2 ) terms in [START_REF] Bresch | On some compressible fluid models: Korteweg, lubrication and shallow water systems[END_REF].

Under assumption (4) and if constants λ, µ, m, M are such that λ < 2µ(M -m) -1 , they prove that there exists at least one weak solution of the problem, given u 0 ∈ H, ρ 0 ∈ H 1 (Ω) and g ∈ L 2 0, T ; L 2 (Ω) . If additionally u 0 ∈ V, then the problem has a unique strong solution locally in time. This solution is global in the two-dimensional case. In [START_REF] Da Veiga | Diffusion on viscous fluids. Existence and asymptotic properties of solutions[END_REF], Beirão da Veiga considers the same problem with all the terms in ( 5) and he proves the existence of a unique local solution for arbitrary initial data and external forces and the existence of a unique global strong solution for small initial data and external forces. In particular, if g = 0, the solution decay exponentially in time to the equilibrium solution with zero velocity field. Second, we consider κ > 0. In [START_REF] Kostin | Modelling of miscible liquids with the Korteweg stress[END_REF], Kostin et The aim of this work is to prove the existence of a global in time weak solution (see Definition 3.1) of the Kazhikhov-Smagulov-Korteweg model (1) for arbitrary initial data and external force field. Our main result reads: 4), T > 0 and g ∈ L 2 0, T ; L 2 (Ω) . Then there exists a weak solution (u, ρ) of (1) global in time such that

Theorem 1.1 Let u 0 ∈ H, ρ 0 ∈ H 1 (Ω) satisfy (
u ∈ L ∞ 0, T ; H ∩ L 2 0, T ; V , ρ ∈ L ∞ 0, T ; H 1 (Ω) ∩ L 2 0, T ; H 2
N , with finite and uniformly bounded energy such that ∀t ≤ T,

ρ(t)u(t) 2 L 2 (Ω) +κ ∇ρ(t) 2 L 2 (Ω) + t 0 λm 4 ∇u(s) 2 L 2 (Ω) +2κλ ∆ρ(s) 2 L 2 (Ω) ds ≤ √ ρ 0 u 0 2 L 2 (Ω) +κ ∇ρ 0 2 L 2 (Ω) + CM 2 λm T 0 g(s) 2 L 2 (Ω)
ds.

The outline of the paper is as follows. In section 2 we explain the derivation of the KSK model [START_REF] Antontsev | Boundary value problems in mechanics of nonhomogeneous fluids[END_REF]. In section 3 we present the functional setting. In section 4 we recall some preliminary results and in section 5 we prove the theorem 1.1.

Modeling of the KSK model

In this section we explain the derivation of the Kazhikhov-Smagulov-Korteweg model [START_REF] Antontsev | Boundary value problems in mechanics of nonhomogeneous fluids[END_REF]. The mixture of two fluids is described by the density ρ(t, x) ≥ 0, the mass velocity field v(t, x) and the pressure p(t, x), depending on the time and space variables (t, x) ∈ [0, T ] × Ω. Joseph [START_REF] Joseph | Fluid dynamics of two miscible liquids with diffusion and gradient stresses[END_REF], Galdi et al. [START_REF] Galdi | Mathematical problems for miscible, incompressible fluids with Korteweg stresses[END_REF] and Joseph et al. [START_REF] Joseph | Non-solenoidal velocity effects and Korteweg stresses in simple mixture of incompressible liquids[END_REF] have developed a theory of non-solenoidal velocity effects and Korteweg stresses in a simple mixture of incompressible liquids. According to these authors, we consider the Korteweg equations for generalized incompressible fluids whose density and volume change with the concentration of one component in the binary mixture, but not with the pressure. This theory may be framed in terms of a mass velocity field v which is not solenoidal (div v = 0), and a volume velocity field u which is solenoidal (div u = 0). Assuming that each fluid is incompressible, the mass density is conserved in the absence of diffusion. The theory of Korteweg, introduced in [START_REF] Korteweg | Sur la forme que prennent les équations du mouvement des fluides si l'on tient compte des forces capillaires causées par des variations de densité considérables mais continues et sur la théorie de la capillarité dans l'hypothèse d'une variation continue de la densité[END_REF] following earlier work by Van der Waals, considers the possibility that motions can be driven by additional stresses associated with gradients of density and gradients of the concentration of species. Dynamical effects which mimic surface tension can arise in thin mixing layers where the gradients of composition are large. It is worth pointing out that similar theories based on thermodynamic arguments are introduced by Brenner [START_REF] Brenner | Navier-Stokes revised[END_REF] (see also the references therein). The author emphasizes the assumption that density gradients give rise to stresses, even when these stresses are induced in single component liquids by temperature gradients. On one hand, in order to model the fluid capillarity effects, Korteweg introduced in the usual compressible fluid model a specific stress tensor which depends on density derivatives. Following the rigorous formulation introduced by Dunn and Serrin [START_REF] Dunn | On the thermomechanics of interstitial working[END_REF] (see also Bresch et al. [START_REF] Bresch | On some compressible fluid models: Korteweg, lubrication and shallow water systems[END_REF]) and neglecting thermal fluctuations, the model reads

∂ t ρ + div (ρv) = 0, ∂ t (ρv) + div (ρv ⊗ v) = ρg + div (S + K), (6)
with the viscous stress tensor S and the Korteweg stress tensor K given by:

S = (ν div v -p)I + 2µD(v), K = (α∆ρ + β|∇ρ| 2 )I + δ(∇ρ ⊗ ∇ρ) + γD 2 x ρ, (7)
where D(v) = (∇v + ∇v T )/2 is the strain tensor and D 2 x ρ is the hessian matrix of the density ρ. Here, the pressure p and the coefficients α, β, γ, δ, µ and ν are functions of ρ. In the viscous stress tensor S, we define the coefficients ν and µ as in [START_REF] Bresch | Effects of density dependent viscosities on multiphasic incompressible fluid models[END_REF] ν = -(µ 0 + λ), µ = µ 0 ρ, where µ 0 and λ are viscosities of the fluid with µ 0 ≥ 0 and 3 λ + 2µ 0 ≥ 0 and λ + 2 3 µ 0 is called the bulk viscosity. Therefore, the viscous stress tensor reads

S = -(µ 0 + λ) div v + p I + 2µ 0 ρD(v), (8)
and it follows that

div S = -(µ 0 + λ) ∇(div v) -∇p + 2µ 0 div ρD(v) . (9)
For the Korteweg stress tensor K, we consider the following relations:

α = κρ, β = κ 2 , δ = -κ, γ = 0,
for some constant κ > 0. This choice corresponds precisely to the Korteweg's original assumptions connected with the variational theory of Van der Waals. In this case, the Korteweg stress tensor reads

K = κ 2 (∆ρ 2 -|∇ρ| 2 )I -κ(∇ρ ⊗ ∇ρ), (10)
and we obtain

div K = κρ∇(∆ρ). (11)
On the other hand, the non-solenoidal velocity effects and the Korteweg stresses arise in process of slow diffusion on miscible incompressible fluids, for example water and glycerin. In the physical experiment given by Joseph [START_REF] Joseph | Fluid dynamics of two miscible liquids with diffusion and gradient stresses[END_REF], the presence of sharp interfaces resembles familiar shapes which can be seen in immiscible liquids. The two fluids are characterized by their reference mass density: we denote by ρ1 the density of the first specie and ρ2 the density of the second one. We also need the velocity field of each constituent: v 1 (t, x) and v 2 (t, x), respectively. We define the volume fraction of the first specie:

φ(t, x) = lim r→0
Volume occupied at time t by the first specie in B(x, r)

|B(x, r)| , with 0 ≤ φ(t, x) ≤ 1.
Therefore, assuming that each fluid is incompressible and keeps a constant mass density, the density of the mixture is defined by

ρ(t, x) = ρ2 1 -φ(t, x) :=ρ 2 (t,x) + ρ1 φ(t, x) :=ρ 1 (t,x) = ρ2 + (ρ 1 -ρ2 )φ(t, x).
Writing the mass conservation for the two species, we obtain

∂ t ρ + div (ρv) = 0, where ρv(t, x) = (ρ 2 v 2 + ρ 1 v 1 )(t, x)
defines the mean mass velocity (or barycentric velocity) v(t, x), which is not divergence free. By contrast, we define the mean volume velocity

u(t, x) = 1 -φ(t, x) v 2 (t, x) + φ(t, x)v 1 (t, x).
Using the definitions, we easily verify that the velocity field u is solenoidal.

According to Kazhikhov and Smagulov [START_REF] Kazhikhov | The correctness of boundary value problems in a diffusion model in an inhomogeneous fluid[END_REF] (see also [START_REF] Antontsev | Boundary value problems in mechanics of nonhomogeneous fluids[END_REF][START_REF] Brenner | Navier-Stokes revised[END_REF]) we consider the following non-standard constraint associated to the pressure p:

div v = -div λ∇ ln(ρ) , (12)
where λ > 0 is interpreted as a diffusion coefficient. This Fick's law describes the diffusive fluxes of one fluid into the other. Clearly, when we set

v = u -λ∇ ln(ρ), (13)
the relation yields [START_REF] Korteweg | Sur la forme que prennent les équations du mouvement des fluides si l'on tient compte des forces capillaires causées par des variations de densité considérables mais continues et sur la théorie de la capillarité dans l'hypothèse d'une variation continue de la densité[END_REF].

Let us prove that if we use the relation (13), then system (6) completed by ( 9) and ( 11) can be rewritten as the KSK model (1). More precisely, let us prove the following:

Lemma 2.1 Using the relation [START_REF] Kostin | Modelling of miscible liquids with the Korteweg stress[END_REF] with div u = 0 and λ = 2µ 0 , then system (6) may be rewritten as system (1), with the modified pressure

P = p - λ 2 (λ + 2 λ) ∆ ln(ρ) -κ ρ∆ρ.
Proof. We use [START_REF] Kostin | Modelling of miscible liquids with the Korteweg stress[END_REF] into equations involving ρ and v, in order to eliminate v. Clearly, the mixture density ρ satisfies the mass conservation and we obtain

∂ t ρ + div (ρu) = div (λ∇ρ). (14)
For the momentum equation, we develop each term of (6) 2 taking into account ( 9) and [START_REF] Joseph | Non-solenoidal velocity effects and Korteweg stresses in simple mixture of incompressible liquids[END_REF]. We start by using (6) 1 in order to write

∂ t (ρv) + div (ρv ⊗ v) = ρ ∂ t v + (v • ∇)v .
Then, using [START_REF] Kostin | Modelling of miscible liquids with the Korteweg stress[END_REF] we have

ρ ∂ t v + (v • ∇)v = ρ ∂ t u + (u • ∇)u -λρ∇( ∂ t ρ ρ ) -λρ(u • ∇)∇ ln ρ -λ(∇ρ • ∇)u + λ 2 (∇ρ • ∇)∇ ln ρ, with ρ∇( ∂ t ρ ρ ) = λ∇∆ρ -∇(u • ∇ρ) -λ ∆ρ ρ ∇ρ + (u • ∇ρ) ∇ρ ρ , ρ(u • ∇)∇ ln ρ = (u • ∇)∇ρ -(u • ∇ρ) ∇ρ ρ , (∇ρ • ∇)∇ ln ρ = ∇ρ • ∇ ∇ρ ρ - |∇ρ| 2 ρ 2 ∇ρ.
Also, taking into account that

(u • ∇)∇ρ -∇(u • ∇ρ) = -div ρ∇u T ,
we obtain

∂ t (ρv) + div (ρv ⊗ v) = ρ ∂ t u + (u • ∇)u -λ(∇ρ • ∇)u + λ div ρ∇u T + λ 2 ∆ρ ρ ∇ρ + ∇ρ • ∇ ∇ρ ρ -∇∆ρ - |∇ρ| 2 ρ 2 ∇ρ .
Furthermore, using ( 12) and ( 13), each term of div S reads

(µ 0 + λ) ∇(div v) = -λ(µ 0 + λ)∇∆ ln(ρ), µ 0 div ρ∇v = µ 0 div ρ∇u + µ 0 λ ∆ρ ρ ∇ρ + ∇ρ • ∇ ∇ρ ρ -∇∆ρ - |∇ρ| 2 ρ 2 ∇ρ , µ 0 div ρ∇v T = µ 0 div ρ∇u T + µ 0 λ ∆ρ ρ ∇ρ + ∇ρ • ∇ ∇ρ ρ -∇∆ρ - |∇ρ| 2 ρ 2 ∇ρ .
Finally, by gathering all these relations and (11) together, the momentum equation ( 6) 2 is rewritten as:

ρ ∂ t u + (u • ∇)u -λ(∇ρ • ∇)u + λ div ρ∇u T -2µ 0 div ρD(u) -λ(µ 0 + λ)∇∆ ln(ρ) + ∇p + λ 2 -2µ 0 λ ∆ρ ρ ∇ρ + ∇ρ • ∇ ∇ρ ρ -∇∆ρ - |∇ρ| 2 ρ 2 ∇ρ ( * ) = ρ g + κ ρ∇(∆ρ). ( 15
)
However, as in [START_REF] Bresch | Effects of density dependent viscosities on multiphasic incompressible fluid models[END_REF], equation ( 15) can be simplified if we make the special choice λ = 2µ 0 . For this choice, the O(λ 2 ) terms marked with (*) in [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] disappear. Moreover, we observe that κρ∇(∆ρ) = κ∇(ρ∆ρ) -κ∇ρ∆ρ.

In conclusion, equation (1) 1 is obtained by including all the gradient terms of (15) in the modified pressure P .

Remark 2.2

The model (1), without the Korteweg stress tensor (i.e. for κ = 0), is called the pollutant model (see [START_REF] Bresch | Effects of density dependent viscosities on multiphasic incompressible fluid models[END_REF]). In particular, (1) is obtained without any assumption on the diffusivity λ.

Remark 2.3

In [START_REF] Galdi | Mathematical problems for miscible, incompressible fluids with Korteweg stresses[END_REF][START_REF] Joseph | Fluid dynamics of two miscible liquids with diffusion and gradient stresses[END_REF][START_REF] Joseph | Non-solenoidal velocity effects and Korteweg stresses in simple mixture of incompressible liquids[END_REF] the solenoidal velocity u is related to v and the gradient of volume fraction φ by setting v = uλ(φ)∇φ, with a diffusion coefficient λ(φ) > 0. This relation is justified by considering that the total mass of solute plus solvent is conserved for each volume Ω, but solute and solvent diffuse in and out of Ω. Consequently, the volume fraction φ of the solute satisfies

∂ t φ + div (φv) = div λ(φ)∇φ . ( 16
)
The equivalence between ( 14) and ( 16) can be seen in [START_REF] Calgaro | Modeling and simulation of mixture flows: Application to powdersnow avalanches[END_REF].

Functional setup

Let us introduce the following functional spaces (see for instance [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF], [START_REF] Temam | Navier-Stokes equations, theory and numerical analysis[END_REF] for their properties):

V = u ∈ D(Ω) 3 : div u = 0 in Ω , V = u ∈ H 1 0 (Ω) : div u = 0 in Ω , H = u ∈ L 2 (Ω) : div u = 0 in Ω, u • n = 0 on Γ , H s N = ρ ∈ H s (Ω) : ∂ρ ∂n = 0 on Γ, Ω ρ(x) dx = Ω ρ 0 (x) dx , s ≥ 2.
V and H are the closures of V in H 1 0 (Ω) and L 2 (Ω) respectively.

Let us recall the definition of a weak solution for the KSK model [START_REF] Antontsev | Boundary value problems in mechanics of nonhomogeneous fluids[END_REF]. Such class of solutions can be found in [START_REF] Bresch | Effects of density dependent viscosities on multiphasic incompressible fluid models[END_REF] for the Kazhikhov-Smagulov type models and in [START_REF] Temam | Navier-Stokes equations, theory and numerical analysis[END_REF] for the incompressible Navier-Stokes equations.

Another equivalent weak formulation can be found in [START_REF] Antontsev | Boundary value problems in mechanics of nonhomogeneous fluids[END_REF].

Definition 3.1 A pair of functions (u, ρ) is called a weak solution of problem (1),( 2),(3) on Ω if and only if the following assumptions are satisfied:

1. u ∈ L ∞ 0, T ; H ∩ L 2 0, T ; V , ρ ∈ L ∞ 0, T ; H 1 (Ω) ∩ L 2 0, T ; H 2 N , 0 < m ≤ ρ(t, x) ≤ M < +∞, a.e. (t, x) ∈ Q T .
2. For all φ ∈ C 1 ([0, T ]; V) such that φ(T, .) = 0, one has:

T 0 -u, ρ∂ t φ + (ρu -λ∇ρ) • ∇ φ + λ 2 ρ(∇u -∇u T ), ∇φ dt = T 0 ρg -κ∆ρ∇ρ, φ dt + ρ 0 u 0 , φ(0) . ( 17
)
3. For all ϕ ∈ C 1 ([0, T ]; H 1 (Ω)) such that ϕ(T, .) = 0, one has:

T 0 u • ∇ρ, ϕ + λ ∇ρ, ∇ϕ -ρ, ∂ t ϕ dt = ρ 0 , ϕ(0) . ( 18
)
Remark 3.2 The pressure P associated with the weak solution (u, ρ) can be obtained using [START_REF] Sy | Local strong solution for the incompressible Korteweg model[END_REF] and the Rham's lemma [START_REF] Temam | Navier-Stokes equations, theory and numerical analysis[END_REF].

Preliminary results

In this section we recall some useful results concerning the convection-diffusion equation satisfied by the density function. This kind of results was used in the framework of compressible Navier-Stokes equations by Feireisl et al. [START_REF] Feireisl | On the existence of globally defined weak solutions to the Navier-Stokes equations[END_REF] and also by Bresch et al. [START_REF] Bresch | Effects of density dependent viscosities on multiphasic incompressible fluid models[END_REF] for the Kazhikhov-Smagulov model. Given the initial density ρ 0 and the velocity field u, we consider the density ρ as solution of the following Neumann problem:

         ∂ t ρ + u • ∇ρ = λ ∆ρ in Q T , ρ(0, x) = ρ 0 (x) in Ω, ∂ρ ∂n = 0 on Σ. (19)
The density ρ satisfies the maximum principle. This result is classical and may be found in [START_REF] Bresch | Effects of density dependent viscosities on multiphasic incompressible fluid models[END_REF].

Proposition 4.1 If (u, ρ) is a weak solution of (1), then

0 < m ≤ ρ(t, x) ≤ M < +∞ a.e. (t, x) ∈ Q T . ( 20
) Proposition 4.2 Let ρ 0 ∈ H 1 (Ω) satisfy (4) and u ∈ C [0, T ]; V ∩ H 2 (Ω) .
Then there exists a unique solution of (19

) such that ρ ∈ L ∞ 0, T ; H 1 (Ω) ∩ L 2 0, T ; H 2 N . Moreover, we have sup 0≤t≤T ρ(t) 2 L 2 (Ω) ≤ ρ 0 2 H 1 (Ω) , (21)
T 0 ∇ρ(t) 2 L 2 (Ω) dt ≤ 1 2λ ρ 0 2 H 1 (Ω) , (22)
sup 0≤t≤T ∇ρ(t) 2 L 2 (Ω) ≤ C λ ρ 0 2 H 1 (Ω) 1 + sup 0≤t≤T u(t) 2 L ∞ (Ω) , ( 23
) T 0 ∆ρ(t) 2 L 2 (Ω) dt ≤ C λ λ ρ 0 2 H 1 (Ω) 1 + sup 0≤t≤T u(t) 2 L ∞ (Ω) , (24)
where C λ is a positive constant depending only on λ.

For a given ρ 0 ∈ H 1 (Ω) which satisfies (4) and u ∈ C [0, T ]; V ∩ H 2 (Ω) , let ρ the solution obtained by Proposition 4.2. Hence, it is clear that the following map

S : C [0, T ]; V ∩ H 2 (Ω) -→ L ∞ 0, T ; H 1 (Ω) ∩ L 2 0, T ; H 2 N ,
such that ρ = Su, is well defined.

Proposition 4.3 Let ρ 0 ∈ H 1 (Ω) satisfy (4) and u 1 , u 2 ∈ C [0, T ]; V ∩ H 2 (Ω) . Set ρ = ρ 1 -ρ 2 = Su 1 -Su 2 and u = u 1 -u 2 ,
we have the following estimates:

sup 0≤t≤T ρ(t) 2 L 2 (Ω) + λ T 0 ∇ρ(t) 2 L 2 (Ω) dt ≤ M 2 λ T sup 0≤t≤T u(t) 2 L 2 (Ω) , (25)
sup 0≤t≤T ∇ρ(t) 2 L 2 (Ω) + λ T 0 ∆ρ(t) 2 L 2 (Ω) dt ≤ 2 λ T sup 0≤t≤T ∇ρ 1 (t) 2 L 2 (Ω) sup 0≤t≤T u(t) 2 L ∞ (Ω) + 2M 2 λ 3 T sup 0≤t≤T u 2 (t) 2 L ∞ (Ω) sup 0≤t≤T u(t) 2 L 2 (Ω) . (26)
Also, we recall that there exists an orthonormal basis of L 2 (Ω) defined by

ω k ∈ V ∩ H 2 (Ω) -P∆ω k = λ k ω k on Ω,
where P is the orthogonal projection operator of L 2 (Ω) onto H. For any n ∈ N * , we denote by X n the finite dimensional subspace of H such that

X n = Vect{ω k , k = 1, . . . , n},
and we consider the orthogonal projection

P n : L 2 (Ω) → X n defined by ∀w ∈ H, P n w, v = w, v , ∀v ∈ X n . ( 27
)
The operator P n satisfies the following property:

Proposition 4.4 P n is continuous on H and it can be extended as a continuous operator on V and V .

Following Feireisl et al. [START_REF] Feireisl | On the existence of globally defined weak solutions to the Navier-Stokes equations[END_REF], we introduce a family of operators M[ρ] : X n -→ X n defined by

M[ρ]v, ω = Ω ρ v • ω dx for all v, ω ∈ X n . ( 28
)
It is clear that if ρ ∈ L ∞ (Ω), then M[ρ] is well defined. Moreover, assumed m > 0, we set

D = ρ ∈ L ∞ (Ω); ρ(x) ≥ m > 0 .
Proposition 4.5 M[ρ] is one-to-one and its inverse satisfies

M[ρ] -1 L(Xn,Xn) ≤ inf x∈Ω ρ(x) -1 ∀ρ ∈ D, (29)
M[ρ 1 ] -1 -M[ρ 2 ] -1 L(Xn,Xn) ≤ C n m 2 ρ 1 -ρ 2 L 2 (Ω) ∀ρ 1 , ρ 2 ∈ D, ( 30
)
where C n is a constant depending on the dimension of X n .

5 Proof of Theorem 1.1

The Faedo-Galerkin method

The approximate solutions

(u n , ρ n ) ∈ C [0, T ]; X n × C [0, T ]; H 1 (Ω) ∩ H 2 N
we shall look for are required to satisfy

                         Ω ∂ t ρ n u n • vdx + Ω (u n • ∇ρ n ) u n • vdx -λ Ω ∆ρ n u n • vdx + Ω (ρ n u n -λ∇ρ n ) • ∇ u n • vdx - λ 2 Ω div ρ n ∇u n -ρ n ∇u T n • vdx = Ω ρ n g • vdx -κ Ω ∆ρ n ∇ρ n • vdx, ∀v ∈ X n , Ω ∂ t ρ n η dx + Ω u n • ∇ρ n η dx = λ Ω ∆ρ n η dx, ∀η ∈ H 1 (Ω), u n (0) = u 0n = P n u 0 , ρ n (0) = ρ 0 . (31) We set N [u n , ρ n ] = -(ρ n u n -λ∇ρ n )•∇ u n + λ 2 div ρ n ∇u n -ρ n ∇u T n -(u n •∇ρ n )u n +λ∆ρ n u n -κ∆ρ n ∇ρ n +ρ n g.
(32) Taking (31) 1 with v = ω k , for k = 1, . . . , n, and integrating in time between 0 and t ≤ T , the solution u n satisfies the following integral equations for k = 1, . . . , n:

Ω ρ n (t)u n (t) • ω k dx = Ω q 0 • ω k dx + t 0 Ω N [u n , ρ n ] • ω k dx ds, ( 33
)
where ρ n = Su n and q 0 = ρ 0 u 0n . Using ( 27) and (28), we rewrite (33) as follows:

M[ρ n (t)]u n (t), ω k = P n q 0 , ω k + P n t 0 N [u n (s), ρ n (s)] ds, ω k ,
for k = 1, . . . , n. Therefore, for all t ∈ [0, T ], we have

M[ρ n (t)]u n (t) = P n q 0 + P n t 0 N [u n (s), ρ n (s)] ds .
Since M[ρ n ] is invertible, then the resulting equation reads

u n ∈ C [0, T ]; X n , u n (t) = M[ρ n (t)] -1 P n q 0 + t 0 N [u n (s), ρ n (s)] ds . (34)
At this stage, u n appears as a fixed point of a suitable functional

Ψ Ψ : C [0, T ]; X n -→ C [0, T ]; X n u n -→ Ψ u n defined by Ψ u n (t) = M[ρ n (t)] -1 P n q 0 + t 0 N [u n (s), ρ n (s)] ds , for all t ∈ [0, T ].
Let X T be the Banach space C [0, T ]; X n endowed with the norm

u n X T = sup 0≤t≤T u n (t) L 2 (Ω) .
In order to apply the Banach fixed point theorem, we shall establish some uniform estimates for Ψ. In the sequel, C is constant and it changes from one line to another and may depend on different parameters that will be precise later on.

,

together with Proposition 4.1, we get

P n ∆ρ n ∇ρ n L 2 (Ω) ≤ C n √ M ∆ρ n 3/2 L 2 (Ω)
.

Employing the same arguments for the other terms of (32), we obtain

P n N [u n , ρ n ] L 2 (Ω) ≤ C u n 2 L 2 (Ω) + u n L 2 (Ω) ∇ρ n L 2 (Ω) + u n L 2 (Ω) + u n L 2 (Ω) ∆ρ n L 2 (Ω) + ∆ρ n 3/2 L 2 (Ω) + u n 2 L 2 (Ω) ∇ρ n L 2 (Ω) + g L 2 (Ω)
.

Using (29) in (37), we obtain for t ≤ T

Ψ u n (t) L 2 (Ω) ≤ M m u 0 L 2 (Ω) + C m T u n 2 X T + u n X T t 0 ∇ρ n L 2 (Ω) ds + u n 2 X T t 0 ∇ρ n L 2 (Ω) ds + T u n X T + u n X T t 0 ∆ρ n L 2 (Ω) ds + t 0 ∆ρ n 3/2 L 2 (Ω) ds + t 0 g L 2 (Ω) ds ,
and taking into account ( 22) and (24), we get

Ψ u n X T ≤ M m u 0 L 2 (Ω) + C m T u n 2 X T + T 1 2 √ 2λ ρ 0 H 1 (Ω) u n X T + C λ T 1 2 √ λ ρ 0 H 1 (Ω) u n X T 1+ u n 2 X T 1/2 + T u n X T + T 1 2 √ 2λ ρ 0 H 1 (Ω) u n 2 X T +T 1 4 C λ λ ρ 0 2 H 1 (Ω) (1+ u n 2 X T ) 3/4 + T 1 2 g L 2 (0,T ;L 2 (Ω))
.

Then, several calculations similar to those employed in (41) give us

T 0 P n N [u 1 n (t), ρ 1 n (t)]-N [u 2 n (t), ρ 2 n (t)] L 2 (Ω) dt ≤ C max(T, T 1 4 ) 1+ u 1 n 2 X T + u 2 n 2 X T u 1 n -u 2 n X T .
Finally, gathering all the preceding estimates together, we find (36).

From now on, we set

R = 2 M m u 0 L 2 (Ω) and B T R = u ∈ X T , u X T ≤ R .
Proposition 5.2 There exists T n ∈]0, 1[ small enough and u n ∈ B Tn R such that

u n = Ψ(u n ). Proof. Let 0 < T n < 1 such that max CT 1 4 n R + 1 R , CT 1 4 n 1+ u 0 L 2 (Ω) +2R 2 ≤ 1 2 .
Using Proposition 5.1, we prove that Ψ is a contraction mapping on B Tn R and we deduce the existence of a unique fixed point of Ψ.

It is clear that u n , the fixed point of Ψ obtained in Proposition 5.2, yields that (u n , ρ n = Su n ) is a local solution of the approximate problem (31). Now, we will prove that this local solution is, in fact, a global one. For this purpose we establish uniform estimates for (u n , ρ n ) with respect to time.

Proposition 5.3 There exists a constant C > 0 depending on ρ 0 , u 0 , g, M, m, λ, κ, such that for all t ∈ [0, T n )

m u n (t) 2 L 2 (Ω) + λm 4 t 0 ∇u n (s) 2 L 2 (Ω) ds ≤ C, (42)
κ ∇ρ n (t) 2 L 2 (Ω) + 2κλ t 0 ∆ρ n (s) 2 L 2 (Ω) ds ≤ C. (43)
Proof. Firstly, taking v = u n (t) in (31) 1 , we obtain

1 2 d dt Ω ρ n |u n | 2 dx + λ 2 Ω ρ n ∇u n -∇u T n : ∇u n dx = Ω ρ n g • u n dx -κ Ω ∆ρ n ∇ρ n • u n dx. (44)
Next, by (31) 2 we have

∂ t ρ n -λ∆ρ n = -u n • ∇ρ n , in the distribution sens on Q T . (45) Since (u n , ρ n ) ∈ C [0, T n ]; X n × L ∞ 0, T n ; H 1 (Ω) ∩ L 2 0, T n ; H 2 N , we conclude that ∂ t ρ n ∈ L 2 0, T n ; L 2 (Ω) .
Now, we multiply (45) by -κ ∆ρ n (t) and after integration by parts, we deduce that

κ 2 d dt ∇ρ n 2 L 2 (Ω) + κ λ ∆ρ n 2 L 2 (Ω) = κ Ω ∆ρ n u n • ∇ρ n dx. (46)
By adding (44) and ( 46), we get

d dt √ ρ n u n 2 L 2 (Ω) +κ ∇ρ n 2 L 2 (Ω) +2κλ ∆ρ n 2 L 2 (Ω) +λ Ω ρ n ∇u n -∇u T n : ∇u n dx = 2 Ω ρ n g•u n dx.
Using the so-called vorticity tensor W = ∇u n -∇u T n (see [START_REF] Bresch | Effects of density dependent viscosities on multiphasic incompressible fluid models[END_REF]), and observing that

Ω ρ n ∇u n -∇u T n : ∇u n dx = Ω ρ n W : ∇u n dx = 1 2 Ω ρ n |rot u n | 2 dx,
we find

d dt √ ρ n u n 2 L 2 (Ω) +κ ∇ρ n 2 L 2 (Ω) + 2κλ ∆ρ n 2 L 2 (Ω) + λ 2 Ω ρ n |rot u n | 2 dx = 2 Ω ρ n g • u n dx.
Since ∆ = ∇div -rot rot , we have

Ω |∇u n | 2 dx = Ω |rot u n | 2 dx.
Then, applying the lower bound of ρ n to

Ω ρ n |rot u n | 2 dx ≥ m Ω |rot u n | 2 dx
, and the above identity, we obtain

d dt √ ρ n u n 2 L 2 (Ω) +κ ∇ρ n 2 L 2 (Ω) + 2κλ ∆ρ n 2 L 2 (Ω) + λm 2 ∇u n 2 L 2 (Ω) ≤ 2 Ω ρ n g • u n dx.
In order to estimate the last term, we use Proposition 4.1, the Poincaré and Young inequalities, and we find

2 Ω ρ n g • u n dx ≤ λm 4 ∇u n 2 L 2 (Ω) + CM 2 λm g 2 L 2 (Ω)
.

Finally, we obtain for each t ∈ [0, T n )

d dt √ ρ n u n 2 L 2 (Ω) +κ ∇ρ n 2 L 2 (Ω) + 2κλ ∆ρ n 2 L 2 (Ω) + λm 4 ∇u n 2 L 2 (Ω) ≤ CM 2 λm g 2 L 2 (Ω)
.

Obviously, thanks to Proposition 5.3, we have the following result:

Corollary 5.4 (u n , ρ n ) is a global solution of (31) and for all T > 0,

(u n ) n is bounded in L ∞ 0, T ; H ∩ L 2 0, T ; V , (47)
(ρ n ) n is bounded in L ∞ 0, T ; H 1 (Ω) ∩ L 2 0, T ; H 2 N . (48)

Uniform estimates for time derivatives

Now, we want to establish uniform estimates for time derivatives ∂ t ρ n and ∂ t u n .

Proposition 5.5 Let T > 0, then the sequence

(∂ t ρ n ) n is bounded in L 4/3 0, T ; L 2 (Ω) .
Proof. Taking the L 2 -norm of ∂ t ρ n , using the Hölder and Gagliardo-Nirenberg inequalities and (38), we get

∂ t ρ n L 2 (Ω) ≤ λ ∆ρ n L 2 (Ω) + C 0 C GN u n 1/4 L 2 (Ω) ∇u n 3/4 L 2 (Ω) ρ n 1/2 L ∞ (Ω) ∆ρ n 1/2 L 2 (Ω)
.

Using the uniform estimate (42) and (20), we obtain

∂ t ρ n L 2 (Ω) ≤ λ ∆ρ n L 2 (Ω) + C ∇u n 3/4 L 2 (Ω) ∆ρ n 1/2 L 2 (Ω) . (49)
Using the Young inequality ab ≤ 1 2 (a 2 + b 2 ) in (49), we get

∂ t ρ n L 2 (Ω) ≤ λ ∆ρ n L 2 (Ω) + C ∇u n 3/2 L 2 (Ω) (50)
It is clear that the second hand side in L 4 3 (0, T ) since we have the uniform time estimates (42) and (43) for ∇u n L 2 (Ω) and ∆ρ n L 2 (Ω) respectively.

Following [START_REF] Bresch | Effects of density dependent viscosities on multiphasic incompressible fluid models[END_REF], we establish an estimation of the fractional time derivative of u n .

Proposition 5.6 Let 0 < δ < T such that T -δ 0 u n (t + δ) -u n (t) 2 L 2 (Ω) dt ≤ C δ 1 4 , (51)
with C independent of n and δ.

Proof. For all functions φ ∈ X T , the approximate solution (u n , ρ n ) satisfies the following equation:

d dτ Ω ρ n u n • φdx - Ω ρ n u n • ∂φ ∂τ dx - Ω ρ n (u n • ∇)φ • u n dx + λ Ω (∇ρ n • ∇)φ • u n dx + λ 2 Ω ρ n ∇u n : ∇φdx - λ 2 Ω ρ n ∇u T n : ∇φdx = Ω ρ n g • φdx -κ Ω ∆ρ n ∇ρ n • φdx. (52)
Integrating (52) with respect to τ between t and t + δ, and taking φ = u n (t + δ) -u n (t), we obtain

Ω ρ n (t + δ) u n (t + δ) -ρ n (t) u n (t) u n (t + δ) -u n (t) dx = t+δ t Ω ρ n (τ ) u n (τ ) • ∇ u n (t + δ) -u n (t) • u n (τ )
+ ∇u n (t) 2 L 2 (Ω)
Analogously, one can obtain the desired estimates of I j (t) terms, for j = 4, . . . , 7.

Finally, if we choose 0 < δ < 1 and thanks to Propositions 5.3 and 5.5, then by gathering all the above estimates together, we rewrite (54) as follows:

ρ n (t + δ) u n (t + δ) -u n (t) 2 L 2 (Ω)
≤ Cδ ∇u n (t + δ) 2

L 2 (Ω)

+ ∇u n (t) 2 L 2 (Ω)
.

Using the lower bound of ρ n and Proposition 5.3, we conclude the proof.

  al. study the mathematical model consisting of a convectiondiffusion equation for the concentration and homogeneous Navier-Stokes equations with the Korteweg stress. They prove the global existence and uniqueness of the solution for the initial-boundary value problem in a two-dimensional bounded domain and its asymptotic for large time. In [17], Sy et al. take the homogeneous Navier-Stokes equations with the Korteweg stress coupled with the mass conservation equation. Their result, given in a three-dimensional bounded domain, concerns the existence of a unique local solution for any initial data and the existence of a unique global solution for small initial data.

  (t + δ) u n (t + δ) -ρ n (t) u n (t) = ρ n (t + δ) u n (t + δ) -u n (t) + ρ n (t + δ) -ρ n (t) u n (t), the relation (53) becomes ρ n (t + δ) u n (t + δ) -u n (t) 2 (t + δ) -ρ n (t) u n (t + δ) -u n (t) • u n (t)dx n (τ ) u n (τ ) • ∇ u n (t + δ) -u n (t) • u n (τ ) dx dτ ∇ρ n (τ ) • ∇ u n (t + δ) -u n (t) • u n (τ ) dx dτ n (τ )∇u n (τ ) : ∇ u n (t + δ) -u n (t) dx dτ ∆ρ n (τ ) ∇ρ n (τ ) • u n (t + δ) -u n (t) dx dτ = I 1 (t) + I 2 (t) + I 3 (t) + I 4 (t) + I 5 (t) + I 6 (t) + I 7 (t).and using the Hölder and Young inequalities and the embedding H 1 (Ω) ⊂ L 4 (Ω), we obtain|I 1 (t)| ≤ Cδ

	dx dτ ∇ρ (53) Ω Ω Ω Ω t+δ t+δ t+δ t t λ Using the following identity -λ = -+ t -λ -2 t+δ t + λ 2 t+δ t Ω ρ n (τ )∇u T n (τ ) : ∇ u n (t + δ) -u n (t) dx dτ + t+δ t Ω ρ t+δ t Ω (54) t+δ t ∂ρ n ∂τ dτ, 1 4 t+δ t ∂ρ n ∂τ 4 3 L 2 (Ω) dτ 3 4 ∇u n (t + δ) 2 L 2 (Ω) + ∇u n (t) 2 L 2 (Ω) 1 4 t+δ t ∇u n (τ ) 2 L 2 (Ω) dτ 3 4 ∇u n (t + δ) 2 L 2 (Ω) + ∇u n (t) 2 L 2 (Ω) , ρ n L 2 (Ω) |I 3 (t)| ≤ Cδ 3 8 t+δ t ∇u n (τ ) 2 L 2 (Ω) dτ 3 8 1 t+δ t ∆ρ n (τ ) 2 L 2 (Ω) dτ 4 ∇u n (t + δ) 2 L 2 (Ω)

n (τ ) • ∇ u n (t + δ) -u n (t) • u n (τ ) dx dτ λ 2 t+δ t Ω ρ n (τ )∇u n (τ ) : ∇ u n (t + δ) -u n (t) dx dτ + λ 2 t+δ t Ω ρ n (τ )∇u T n (τ ) : ∇ u n (t + δ) -u n (t) dx dτ + t+δ t Ω ρ n (τ ) g(τ ) • u n (t + δ) -u n (t) dx dτ -κ t+δ t Ω ∆ρ n (τ ) ∇ρ n (τ ) • u n (t + δ) -u n (t) dx dτ. Ω ρ n ρ ρ n (τ ) g(τ ) • u n (t + δ) -u n (t) dx dτ -κ

First, let us estimate I 1 (t). Using Hölder inequality, we get

|I 1 (t)| ≤ ρ n (t + δ) -ρ n (t) L 2 (Ω) u n (t + δ) -u n (t) L 4 (Ω) u n (t) L 4 (Ω) .

Writing

ρ n (t + δ) -ρ n (t) = .

In the same manner, we obtain the following estimations:

|I 2 (t)| ≤ Cδ

Proposition 5.1 There exists a constant C > 0 depending on n, λ, κ, M , m, ρ 0 H 1 (Ω) , g L 2 (0,T ;L 2 (Ω)) , such that for all u n ∈ X T ,

and for all u 1 n , u 2 n ∈ X T ,

Proof. We have for all u n ∈ X T and for t ≤ T

We will use several times that on a finite dimensional subspace all the norms are equivalent. For instance we have

Moreover, using the continuity of P n in V , the embedding L 4/3 (Ω) ⊂ H -1 (Ω), the Hölder inequality and the following inequality:

Hence, we have

where C depend on n, λ, κ, M , m, ρ 0 H 1 (Ω) and g L 2 (0,T ;L 2 (Ω)) . Using the Young inequality in (39), we get easily

.

In order to prove (36), let us consider u 1 n , u 2 n ∈ X T and we set ρ

n the solutions obtained by Proposition 4.2. Then, we find

Thanks to Propositions 4.3 and 4.5, we obtain

We rewrite the last two terms of (40) as follows:

For the term K 1 , Proposition 4.5 gives

For the term K 2 , we have

The existence of weak solution (u, ρ)

The final task is to employ the preceding uniform estimates in order to pass to the limit in the approximate problem (31). When n → +∞, we have u 0n -→ u 0 in H strongly.

Thanks to (47) and (48), it is possible to choose the subsequences (u n ) n and (ρ n ) n such that

weakly-star, and

Using these convergence results, we are able to pass to the limit in the linear terms of (31). Finally, to be able to pass to the limit in the nonlinear terms of (31), it is necessary to use the following strong convergence: Proposition 5.7 There exists a subsequence (u n , ρ n ) n which converges strongly to (u, ρ) in L 2 0, T ; L 2 (Ω) × L 2 0, T ; H 1 (Ω) . Moreover (u, ρ) is a weak solution of (1).

Proof. Applying some compactness theorems [18, Chap.3, Theorem 2.1] for ρ n and [16, Theorem 5] for u n and using Propositions 5.5 and 5.6, we obtain the desired result.


Marc Bonnet

Rémi Cornaggia
email: remi.cornaggia@ensta-paristech.fr

Higher-order expansion of misfit functional for defect identification in elastic solids

Keywords: Topological derivative, identification, elastodynamics, asymptotic analysis

Higher-order expansion of misfit functional for defect identification in elastic solids

Marc Bonnet, Rémi Cornaggia

Problem statement. We consider a reference (i.e defect-free) 3D elastic solid Ω characterized by Hooke's tensor C and mass density ρ. The time-harmonic background displacement u then solves u, w C Ωω 2 (u, w) ρ Ω = F(w) ∀w ∈ W, (1) where ., . C D and (., .) ρ D denote the stiffness and mass bilinear forms associated to a given domain D characterized by (C, ρ), W ⊂ H 1 (Ω) is the function space incorporating the relevant essential boundary conditions (if any), the linear form F ∈ W defines the applied time-harmonic loading and ω is the angular frequency.

Assuming the presence of a defect inside Ω, and that we can measure the resulting displacement u ex on a surface Γ, we define the leastsquares cost functional J(w), with the elastodynamic displacement w associated to the given excitation and a known trial defect, by:

J(w) = 1 2 Γ |w(x) -u ex (x)| 2 dS x (2)
We now consider a specific trial defect B a = z + aB (Fig 1), centered at z ∈ Ω, of small size a and reference shape B. It is a perfectly bonded inclusion filled with a material characterized by its Hooke tensor C = C + ∆C and mass density ρ = ρ + ∆ρ. We denote u a the displacement in the perturbed domain, and v a = u au the displacement perturbation. J(u a ) admits the Ω : (C, ρ) exact expansion about u:

(C ⋆ , ρ ⋆ ) z B a (z) B a (z) = z + aB a Γ
J(u a ) = J(u) + J (u; v a ) + J (u; v a , v a ) = J(u) + Γ (u -u ex )•v a + 1 2 Γ |v a | 2 (3)
The goal is now to expand J(u a ) in powers of a.

Similar expansions have been studied in e.g. [START_REF] Bonnet | Inverse acoustic scattering by small-obstacle expansion of a misfit function[END_REF] for rigid obstacles in 3D acoustic media and [START_REF] Silva | Higher order topological derivatives in elasticity[END_REF] for holes in 2D elastic bodies. Define the adjoint field p as the solution of

p, w C Ω -ω 2 (p, w) ρ Ω = J (u; w) ∀w ∈ W. (4) We can then compute J (u; v a ) as J (u; v a ) = -p, u a ∆C Ba + ω 2 (p, u a ) ∆ρ Ba (5)
Expanding v a in powers of a is now needed. As we will see, J (u; v a , v a ) requires only the leading contribution of v a | Γ whereas a higher-order expansion of v a | Ba is needed for evaluating [START_REF] Silva | Higher order topological derivatives in elasticity[END_REF].

Expansion of the solution perturbation. Following e.g. [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF], v a solves the integro-differential Lippmann-Schwinger equation:

L a [v a ](x) = -u, G ∆C Ba + ω 2 (u, G) ∆ρ Ba (6) with L a [v](x) := v(x) + v, G ∆C Ba -ω 2 (v, G) ∆ρ Ba and G = G(•, x)
is the elastodynamic Green's tensor for a unit point force applied at x and satisfying homogeneous boundary conditions consistent with problem (1) on ∂Ω. Substituting the ansatz

v a (x) = aV 1 (x) + a 2 V 2 (x) + 1 2 a 3 V 3 (x) + 1 6 a 4 V 4 (x) + δ a (x),
x ∈ B a (7) (with x := (x-z)/a ∈ B) into (6) and expanding the resulting equation in powers of a (in particular using that G(ξ, x) = a -1 G ∞ ( ξ-x)+O( 1), with G ∞ denoting the static full-space Kelvin fundamental solution) yields a sequence of integral equations for the V j . These equations correspond to elastostatic problems for the normalized inclusion B embedded in an unbounded reference medium, and are solved with the help of Eshelby's equivalent inclusion method [START_REF] Mura | Micromechanics of defects in solids[END_REF].

The remainder δ a in (7) solves an integrodifferential equation of the form L a [δ a ] = γ a . The operator L a : H 1 (B a ) → H 1 (B a ) is shown to be invertible with bounded inverse, while γ a can be estimated as γ a H 1 (Ba) = O(a 11/2 ). Consequently, there exists a constant C > 0 independent of a such that

δ a H 1 (Ba) ≤ Ca 11/2 . ( 8
)
For x / ∈ B a , plugging (7) in the form v a (x) ≈ aV 1 (x) into (6) yields the outer expansion

v a (x) = -a 3 ∇u(z) : A : ∇G(z, x) -ω 2 ∆ρ|B|u(z)•G(z, x) + o(a 3 ), (9)
A being the elastic moment tensor associated B, C and ∆C [START_REF] Ammari | Polarization and moment tensors[END_REF][START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF].

Cost functional expansion. Substituting (7) into ( 5) and ( 9) into (3), J(u a ) is finally found to have an expansion of the form:

J(u a ) = J 6 (a, z) + o(a 6 ) (10) with J 6 (a, z) = J(u) + a 3 T 3 (z) + a 4 T 4 (z) + a 5 T 5 (z) + a 6 T 6 (z),
the o(a 6 ) estimate resulting from (8) and ( 9). The T j (z) are found to be given in terms of (i) the background field u and its derivatives at z, (ii) the adjoint field p and its derivatives at z, (iii) A and other elastic moment tensors that involve the material parameters, the shape B and the angular frequency ω, and (iv) the complementary part of

G, i.e. G -G ∞ = (G ∞,ω -G ∞ )+G C
, where G ∞,ω is the elastodynamic full-space fundamental solution and G C accounts for the boundedness of Ω. In particular, T 3 (z) is the well-known topological derivative:

T 3 (z) = -∇u : A : ∇p -ω 2 ∆ρ|B|u•p (z).
Moreover, T 4 (z) = 0 for any centrally-symmetric shape B. The complementary part G ∞,ω -G ∞ (known analytically) is involved in T 5 (z) and T 6 (z), while G C appears in T 6 (z) only. Since the exact computation of G C would require solving an elastodynamic problem on Ω for each trial location z, we plan to use an approximation method to save computational time.

Closed-form formulae for the T j can be obtained when B is spherical (for which case we provide explicit expressions) or ellipsoidal.

Identification. Following [START_REF] Bonnet | Inverse acoustic scattering by small-obstacle expansion of a misfit function[END_REF], estimates of the location z est and size a est of the real defect can then be sought as minimizers of J 6 (a, z), with z spanning a predefined sampling grid. This entails computing the T j (z) over the sampling grid and minimizing a → J 6 (a, z) for each z, the latter step being very fast and straightforward A preliminary example is set in free space (so that G C = 0) for a spherical scatterer of radius 0.1λ S illuminated by a plane P-wave travelling along the positive x-direction, with a discrete array of displacement sensors lying behind the scatterer. The above procedure yields the size estimate a est ≈ 0.105λ S ; moreover, the estimated location z est is found to be very close to the true center of the scatterer. The contour plot of J 6 (a est , z) (Fig 2 ) shows improved localisation (relative to the topological derivative T 3 (z)) for this partial-aperture configuration.

Figure 1 :

 1 Figure 1: Computational domain and inclusion.

Figure 2 :

 2 Figure 2: T 3 (z) and J 6 (a est , z) plotted in the (XZ) plane around the obstacle (in white).


Member, IEEE Driss Yousfi
email: dr_yousfi@yahoo.com

Abdallah Darkawi

Comparison of two Position and Speed Estimation Techniques used in PMSM Sensorless Vector Control

Keywords: PMSM vector control, sensorless, flux estimator, reduced observer, comparative study

The paper describes the comparison between two different high-performance techniques used for the sensorless estimation of the motor shaft position in Permanent Magnet Synchronous Motor Sensorless Drives. Rotor position and speed are estimated from measured terminal voltages and currents, and are used as feedback in a sensorless vector control scheme, achieving almost the same high-performance of a sensored drive. The paper point out the differences, by using experimental implementation, between an open-loop flux estimator based on the electrical model of the machine, and a reduced flux observer.

Introduction

Vector controlled Permancnt Magnet Synchronous Motor (PMSM) need an encoder or a resolver to correctly align the stator current vector. Such an electromechanical position transducer is not present in DC motor drives or in constant V/Hz induction motor drives, thus disadvantaging PMSM drives in medium and low power applications. To overcome this weakness, several sensorless strategies for PMSM vector control, based on different position and speed estimation techniques have been proposed in last years [START_REF] Johnson | Review of sensorless methods for brushless DC[END_REF]. Generally, they consist of processing stator voltage and current measurements to perform the sensorless estimation of the mechanical quantities Among the proposals in this field, two kinds of approaches seem to be most interesting since they do not need any additional hardware i.e. the electrical model based techniques [START_REF] Yousfi | Indirect Position and Speed Sensing for PMSM Sensorless Control[END_REF] and the state observer based techniques [START_REF] Pillay | Application characteristics of permanent magnet synchronous and brushless DC motors for servo drives[END_REF][START_REF] Safi | Analysis and simulation of the high-speed torque performance of brushless DC motor drives[END_REF][START_REF] Chen | Current-mode control for sensorless BDCM drive with intelligent commutation tuning[END_REF][START_REF] Cheng | Design of a sensorless commutation IC for BLDC motors[END_REF]. From the two mentioned families, the paper consider precisely:

-An open-loop electrical model based flux estimator that was developed in our laboratory [START_REF] Ribeiro | Comparision of carrier signal voltage and current injection for the estimation of flux angle and rotor position[END_REF]. -A reduced flux observer which is derived from an existing full order Luenberger observer [START_REF] Batzel | Electric propulsion with the sensorless permanent magnet synchronous motor: model and approach[END_REF][START_REF] Batzel | Electric Propulsion With Sensorless Permanent Magnet Synchronous Motor: Implementation and Performance[END_REF], and improved to achieve good performance.

The paper, then, wants to compare the two methods in terms of:

• Complexity of the method.

• Position and speed estimation performance.

• Sensorless speed range.

• Immunity against motor parameters uncertainty and electrical measurement errors. • Initial position error cancellation.

Modeling of the PMSM

In this section, a brief description of the PMSM mathematical models are described since both the investigated estimation methods need to manipulate the equations of the machine.

The models of the PMSM in the stationary frame (α-β) and the rotating frame (d-q) are respectively:

dt d i R v dt d i R v β β β α α α Φ Φ + = + = (1) ) sin( i L ) cos( i L m m θ Φ Φ θ Φ Φ β β β α α α + = + = (2)
,and

d q q q q d d d w dt d i R v w dt d i R v Φ Φ Φ Φ + + = - + = (3) q q q m d d d i L i L = + = Φ Φ Φ (4)
,where Φ m is the maximum phase flux linkage of the permanent magnet, L αβ and L dq the inductances in αβ and dq frames, R the winding resistance and θ the actual rotor angle.

These models form the basis for the presented rotor position and speed estimator. As the considered PMSM is a nonsalient machine, all the inductances in those equations are equal .i.e. L α =L β =L d =L q =L.

Open-loop Flux Estimator Based Technique

Position sensing via flux-linkage variation has been known for many years, but its successful implementation has become possible only in the last decade with the emergence of devices with sufficient real-time processing power.

The first investigated technique belongs to this category of estimators. It is based on an open-loop flux linkage estimation and PI-controller. Accurate flux linkages are derived from the back-emf by means of a modified integrator. The estimation algorithm is schematically illustrated in Fig. 1 and involves the following steps [START_REF] Yousfi | Indirect Position and Speed Sensing for PMSM Sensorless Control[END_REF]:

Step 1:

Estimation of the flux linkage in αβ-frame by measuring the phase voltages and currents:

∫ - = dt ) i R v ( ˆαβ αβ αβ Φ ( 5
)
Step 2: Estimation of the motor currents αβ i ˆ:

)] ŝin( [ L 1 i ˆ)] ĉos( [ L 1 i ˆm m θ Φ Φ θ Φ Φ β β α α - = - = (6)
Fig. 1: Flowchart of the position and speed estimator.

Step 3:

Calculation of the difference between the actual currents and the estimated ones in α-β and d-q frames using the Park transformation:

αβ αβ αβ ∆ i î i - = ( 7
) ) ĉos( i ) ŝin( i i q θ ∆ θ ∆ ∆ β α + - = (8)
This last current error is passed through a Low Pass filter to reduce the effect of the noise due to the differentiation.

Step 4:

Finally, a PI-Controller is used to cancel the obtained error by accelerating or decelerating the q d ˆ-frame. Thus, the PI-Controller output is the estimated speed w ˆ.

The position estimate θ ˆ can be deduced through an integration and corrected using Look Up Table procedure that gives a correction term ∆θ ; so:

θ ∆ θ θ + = ĉ (9)

Reduced Flux Observer Based Technique

State observers require the use of a relatively accurate motor model, the measurement of the motor currents (system output) and the knowledge of the feeding voltages (system input). The basic idea is to use the difference between the state variables and the estimated state variables to calculate the rotor position and speed, directly or through related variables. Several sort of approaches are reported in literature i.e. deterministic (Luenberger [START_REF] Pillay | Application characteristics of permanent magnet synchronous and brushless DC motors for servo drives[END_REF][START_REF] Safi | Analysis and simulation of the high-speed torque performance of brushless DC motor drives[END_REF]), non-linear (sliding mode [START_REF] Chen | Current-mode control for sensorless BDCM drive with intelligent commutation tuning[END_REF]) and stochastic (extended Kalman filters [START_REF] Cheng | Design of a sensorless commutation IC for BLDC motors[END_REF]) observers. They exhibit different peculiarities in term of algorithm complexity and sensitivity to parameter variation and noise. The used observer is based on the flux linkage model of the PMSM in a two phase stationary frame as shown in the following equation:

Cx y Bu Ax x = + = & (10) where             - - - = 0 w 0 0 w 0 0 0 0 0 0 0 A e e γ γ γ γ ,             = 0 0 0 0 1 0 0 1 B and         - - = L 1 0 L 1 0 0 L 1 0 L 1 C (11)
The input vector and the full state vector are respectively:

T u v v α β   =   and T m m x α β α β   = Φ Φ Φ Φ   .
The aimed states are: cos and sin

m m e m m e α β θ θ Φ = Φ Φ = Φ .
The output vector is the measured currents: In the state space model w e and θ e are the electrical rotor speed and position. γ = R / L is the inverse of the time constant. To consider the reduced version of the above observer a reduced vector and matrix are introduced as follow:

T y i i α β   =   .
A 11 = J w , A 12 = [0], A 21 = γ I, A 22 = -γ I, B 1 = [0], B 1 =I (12) T m m v x y α β α β     = Φ Φ Φ Φ =      

M

,where

m m v α β Φ   =   Φ   and y α β Φ   =   Φ         - =       = 0 1 1 0 J and 1 0 0 1 I If ẑ v Gy = -
, z the reduced observer is defined by:

[ ] y G z v ˆu ) B G B ( y A G A G ) A G A ( z ) A G A ( z r 2 r
+ = - + - + - + - = & (13)
v is the estimate state and G r is the reduced observer gain and is defined by :

1 2 1 2 2 1 r g g G g I g J g g -   = = +     (14)
Supposing the mechanical quantities slowly changing versus the electrical ones, the model can be then considered as a linear time varying system. Now with appropriate calculation of the matrixes G r , B 0 and K 0 the state error continually vanishes. The rotor position estimation is calculated by the arctg operation:

        Φ Φ = θ α β m m tan Arc ˆ ( 15
)
Fig. 2 shows the resulting observer scheme.

It is clear from the state equations [START_REF] Corley | Rotor position and velocity estimation for a salient-pole permanent magnet synchronous machine at standstill and high speeds[END_REF][START_REF] Ribeiro | Comparision of carrier signal voltage and current injection for the estimation of flux angle and rotor position[END_REF] that rotor speed is required first for the implementation of the flux observer. A simple manner to estimate its magnitude is given by [START_REF] Batzel | Electric propulsion with the sensorless permanent magnet synchronous motor: model and approach[END_REF]: Here, δ is an adjustment coefficient.

Experimental Results

The general block schematic of the laboratory setup used for experimental verification is shown in Fig. 3. It consists of a 500W three-pole PMSM which is fed by an IGBT Voltage Source Inverter. For estimation and control tasks, a dSpace DS1104 board has been used. In order to compare the estimated position and speed signals with the actual output, an incremental encoder is also mounted on the drive shaft. The positional error is measured to be within +/-1.5° in all conditions for the first method and +/-1° for the second. A minor prior adjustment is necessary to allow initial position error cancellation in open-loop flux estimator based technique.

Results Analysis and Comparison

When the mechanical sensor is replaced by an estimation technique that provides the rotor position and speed, all the control methods commonly used in sensored control of the motor are conserved in the sensorless scheme. That is the major advantage shared by the presented sensorless techniques. On the other hand, the major difference between the two techniques results from the flux estimation which is carried out by an open-loop calculation in the first method and by a closed-loop one in the second method. the flux observer based technique is a self-adjusting flux estimation method.

In the following section the effects of different factors on the rotor position and speed estimation accuracy are investigated.

Starting up

Starting difficulty is a challenging problem in sensorless operation of PMSM drives. For consistent starting, a separate technique was needed, often using the inductance-variation approach [START_REF] Takeshita | Sensorless salientpole PM synchronous motor drive in all speed ranges[END_REF] or the high frequency signal injection [START_REF] Corley | Rotor position and velocity estimation for a salient-pole permanent magnet synchronous machine at standstill and high speeds[END_REF][START_REF] Ribeiro | Comparision of carrier signal voltage and current injection for the estimation of flux angle and rotor position[END_REF]. In fact, a closed loop smooth and reliable starting can only be achieved when information regarding the absolute rotor position is available. For a simple starting, the test motor is started up in open-loop operation as a stepper motor with sinus-wave excitation until the rotor speed reaches convenient level for angle estimation. The lowest sensorless speed is 25 rpm for the first technique and 33 rpm for the second.

Speed Estimation

As described in section 3 the rotor speed, in the open-loop flux estimation based method, is obtained by closed-loop estimation using PI-Controller. That is the reason of the good estimation results of this quantity (Fig. 4). However, in the observer estimation technique, the speed is estimated separately and it is considered to be a system parameter [START_REF] Ribeiro | Comparision of carrier signal voltage and current injection for the estimation of flux angle and rotor position[END_REF]. Therefore the position estimation accuracy is dependent on the accuracy of the speed estimate, which depends itself on the PM flux linkage Φ m (16). Not including an adaptive correction scheme, PM flux linkage uncertainty would generate velocity estimation error, particularly in low speed rang. As, in such condition the back-emf signals derived from the measured current and voltage are feeble. Additional procedure is then necessary to remove the effects of this uncertainty.

Effect of Resistance Uncertainty:

In the PMSM, the most temperature sensitive parameter is the stator winding resistance, the fundamental parameter necessary for all the presented techniques. The effect of an uncertainty at the level of this parameter introduces a steady rotor angle estimation error. If the stator current is in phase with the motional back-emf, uncertainty in the stator resistance generates no error in the estimated rotor angle [START_REF] Batzel | Electric Propulsion With Sensorless Permanent Magnet Synchronous Motor: Implementation and Performance[END_REF].

Effect of Stator Voltage and Current Measurement Errors

Both the indirect position sensing methods require the measurement of the stator current and voltage waveforms. The use of these measurements in the rotor position estimator will result in angle estimation inaccuracy owing to the inevitable scaling or offset errors contained in the sensed quantities.

In practice, the open-loop integration in ( 5) is prone to errors caused by drift. Small offset signals in the measurements are summed over time, causing the integrator output to saturate. Furthermore, an initial value problem associated with the pure integrator occurs from the integral value at the initial instant.

Integrator drift can be reduced by improving the pure integrator with a low-pass filter with very low cutoff frequency. Another way to solve this problem consists in replacing the integrator by an alternative integrator structure [START_REF] Yousfi | Indirect Position and Speed Sensing for PMSM Sensorless Control[END_REF][START_REF] Hu | New integration algorithms for estimating motor flux over a wide speed range[END_REF]. But these modifications have two drawbacks: the inhibition of the flux estimator in low-speed range and the phase lags induced by the filter. Consequently an additional position error arises. Based on the above considerations, a Look Up Table procedure is implemented to compensate the position estimation error. The estimation algorithm becomes then more complicated.

In the second technique, selecting the magnitude of the closed-loop eigenvalue reals to be much greater than the angular velocity helps to reduce the effect of voltage measurement offsets on the rotor position estimation accuracy.

Conclusion

When the presented estimation methods are compared in their complete version, the reduced flux observer based technique seems more complex in term of formulation and practical implementation. But in actual fact, the closed-loop flux estimation achieved by this technique makes it more attractive; especially if the problem of the speed estimation in low speed rang is resolved by using a position derivation instead of the back-emf based calculation. Furthermore, with appropriate parameters adjusting the observer based technique presents a good immunity against motor parameters uncertainty and electrical measurement defaults.

Basically, the main limitations of the investigated solutions refer to standstill operation and safe starting due to the unobservability of rotor angle under those operating conditions. In many industry applications, initial positioning is not necessary and a simple open-loop starting up can be used. In the opposite case, starting method should be involved in the whole sensorless algorithm.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Block diagram of the reduced position observer.

Fig. 4 .

 4 Fig. 4. Experimental results of the open-loop flux estimator based technique.

Fig. 5 .

 5 Fig. 5. Experimental results of the reduced flux observer based technique.