Gianluca Frasca
email: gianluca.frasca@univ-rennes1.fr

Frédéric Gueydan

Marc Poujol

Fleurice Parat

Patrick Monie

Alexandre Pichat

Sophie Mazier

Jean-Pierre Brun

Patrick Monié

Fast switch from extensional exhumation to

INTRODUCTION

The Mediterranean is characterized by the coexistence of backarc basins and arcuate subduction zones caused by slab rollback (e.g., Calabria, Hellenic trench; Wortel and Spakman, 2000). The Gibraltar arc (Betics in Spain and Rif in Morocco) is marked by the exhumation of large bodies of sub-continental mantle, the Ronda Peridotites, whose role in the geodynamics of the western Mediterranean is still debated (Garrido et al., 2011;[START_REF] Van Hinsbergen | Origin and consequences of Western Mediterranean subduction, rollback, and slab segmentation[END_REF]. The polycyclic history of the Ronda lithospheric section makes it difficult to unequivocally link slab rollback evolution and sub-continental mantle exhumation, as proposed for the central Mediterranean backarc system [START_REF] Acosta-Vigil | The composition of nanogranitoids in migmatites overlying the Ronda peridotites (Betic Cordillera, S Spain): the anatectic history of a polymetamorphic basement[END_REF][START_REF] Sánchez-Rodríguez | Mesozoic formation of pyroxenites and gabbros in the Ronda area (southern Spain), followed by Early Miocene subduction metamorphism and emplacement into the middle crust: U-Pb sensitive high-resolution ion microprobe dating of zircon[END_REF][START_REF] Bache | Evolution of rifted continental margins: The case of the Gulf of Lions (Western Mediterranean Basin)[END_REF].

The Ronda Peridotites consist of a tectonic sliver within the Alboran crust, bounded by two major tectonic contacts. A Crust-Mantle Extensional Shear Zone at the top of the This article is protected by copyright. All rights reserved.

Ronda Peridotites accommodated extreme continental lithospheric thinning (CMESZ, red in Fig. 1; [START_REF] Afiri | Oligo-miocene exhumation of the Beni-Bousera peridotite through a lithosphere-scale extensional shear zone[END_REF][START_REF] Argles | Attenuation and excision of a crustal section during extensional exhumation: the Carratraca massif, Betic Cordillera, southern Spain[END_REF]Frasca et al., 2016;[START_REF] Précigout | Deformation and exhumation of the Ronda Peridotite (Spain)[END_REF]. The crustal rocks overlying the peridotites comprise a strongly attenuated crustal section [START_REF] Balanyá | Alternating contractional and extensional events in the Alpujarride nappes of the Alboran domain (Betics, Gibraltar arc)[END_REF]. The metamorphic grade decreases continuously from granulite facies to unmetamorphosed rocks in the 4 km thick crustal sequence above the Ronda Peridotites [START_REF] Argles | Attenuation and excision of a crustal section during extensional exhumation: the Carratraca massif, Betic Cordillera, southern Spain[END_REF][START_REF] Balanyá | Alternating contractional and extensional events in the Alpujarride nappes of the Alboran domain (Betics, Gibraltar arc)[END_REF][START_REF] Barich | Microstructures and petrology of melt inclusions in the anatectic sequence of Jubrique (Betic Cordillera, S Spain): Implications for crustal anatexis[END_REF]. Furthermore, shearing in the granulites was coeval and consistent with garnet-spinel mylonitization at the top of the mantle [START_REF] Précigout | Deformation and exhumation of the Ronda Peridotite (Spain)[END_REF]Frasca et al., 2016). However, the age of this continental lithospheric thinning is still debated. In the granulites, garnet rims and part of the matrix crystallized during the Tertiary [START_REF] Massonne | Wealth of P-T-t information in medium-high grade metapelites: Example from the Jubrique Unit of the Betic Cordillera[END_REF]Platt et al., 2003b;[START_REF] Whitehouse | Dating high-grade metamorphism -constraints from rare-earth elements in zircon and garnet[END_REF]. A much older partial melting event has been proposed from ca. 1.4-1.2 GPa to 0.6-0.5 GPa at around 800ºC for part of the Ronda crustal section in Hercynian times [START_REF] Barich | Microstructures and petrology of melt inclusions in the anatectic sequence of Jubrique (Betic Cordillera, S Spain): Implications for crustal anatexis[END_REF][START_REF] Ruiz Cruz | Garnet variety and zircon ages in UHP meta-sedimentary rocks from the Jubrique Zone (Alpujárride complex, Betic Cordillera, Spain): evidence for a pre-alpine emplacement of the Ronda Peridotites[END_REF]. Furthermore, Jurassic ages have been obtained in pyroxenite layers within the peridotites (Sanchez-Rodriguez and Gebauer, 1999).

The Ronda Peridotites Thrust (RPT, blue in Fig. 1; [START_REF] Frasca | Structural record of Lower Miocene westward Alboran Domain motion in the Western Betics (southern Spain)[END_REF][START_REF] Vitale | Deformation characterization of a regional thrust zone in the northern Rif (Chefchaouen, Morocco)[END_REF] marks the second tectonic contact of the peridotites. The rocks below the Ronda Peridotites are highly diverse, with variable degrees of metamorphism from low-grade limestones to strongly mylonitized marbles (see summary in [START_REF] Frasca | Structural record of Lower Miocene westward Alboran Domain motion in the Western Betics (southern Spain)[END_REF]). Migmatites in the footwall and granitic dikes intruding the mantle and the crustal hanging wall yield a cluster of crystallization ages between 23 and 19 Ma [START_REF] Monié | Cooling and exhumation of the Western Betic Cordillera, 40 Ar/ 39 Ar thermochronological constraints on a collapsed terrane[END_REF][START_REF] Sosson | Late exhumation stages of the Alpujarride complex (western Betic Cordilleras, Spain): new thermochronological and structural data on Los Reales and Ojen nappes[END_REF][START_REF] Sánchez-Rodríguez | Mesozoic formation of pyroxenites and gabbros in the Ronda area (southern Spain), followed by Early Miocene subduction metamorphism and emplacement into the middle crust: U-Pb sensitive high-resolution ion microprobe dating of zircon[END_REF]Gebauer, 2000, Esteban et al., 2011;[START_REF] Rossetti | Timing and modes of granite magmatism in the core of the Alboran Domain, Rif chain, northern Morocco: Implications for the alpine evolution of Western Mediterranean[END_REF]. A thermal overprint on Hercynian migmatites in the thrust footwall during the Tertiary has been also reported [START_REF] Barich | Microstructures and petrology of melt inclusions in the anatectic sequence of Jubrique (Betic Cordillera, S Spain): Implications for crustal anatexis[END_REF]. A major shortening event during the Miocene therefore triggered the final crustal emplacement of the Ronda Peridotites [START_REF] Tubía | Sequential development of the metamorphic aureole beneath the Ronda Peridotites and its bearing on the tectonic evolution of the Betic Cordillera[END_REF][START_REF] Précigout | Deformation and exhumation of the Ronda Peridotite (Spain)[END_REF].

Accepted Article

This article is protected by copyright. All rights reserved. These first order geological features have yielded a diversity of tectonic models for the exhumation of the Ronda Peridotites from mantle depths to the surface. The three main models are: 1/ exhumation during the Hercynian orogeny (collision/collapse) with only passive reworking during the Alpine orogeny [START_REF] Massonne | Wealth of P-T-t information in medium-high grade metapelites: Example from the Jubrique Unit of the Betic Cordillera[END_REF][START_REF] Ruiz Cruz | Garnet variety and zircon ages in UHP meta-sedimentary rocks from the Jubrique Zone (Alpujárride complex, Betic Cordillera, Spain): evidence for a pre-alpine emplacement of the Ronda Peridotites[END_REF]; 2/ exhumation during the Tethysian rifting, followed by the subduction of the Tethyan passive margin (and its lithospheric mantle) during the complex Alpine phase [START_REF] Tubía | Sequential development of the metamorphic aureole beneath the Ronda Peridotites and its bearing on the tectonic evolution of the Betic Cordillera[END_REF]Vissers et al., 1993); and 3/ Oligo-Miocene exhumation during westward Alboran slab roll-back, followed by rift-inversion during the lower Miocene (Garrido et al., 2011;[START_REF] Frasca | Structural record of Lower Miocene westward Alboran Domain motion in the Western Betics (southern Spain)[END_REF][START_REF] Précigout | Deformation and exhumation of the Ronda Peridotite (Spain)[END_REF].

In order to discriminate between these models, we provide new ages (U-Pb LA-ICP-MS and Ar-Ar) obtained on syndeformational structures related to both thinning and thrusting. The ages clearly separate two successive crustal melting events structurally associated with the Ronda Peridotites. We compare the new ages with the available datasets in the region and discuss their geodynamic implications.

TWO SUCCESSIVE CRUSTAL MELTING EVENTS

Structural setting

Two types of intrusion were sampled from the Ronda Peridotites and their crustal envelope (for positions of the samples, see the regional map and the N-S cross-section in Fig. 1A,C and a more precise structural location on the schema in Fig. 2A): 1-a syndeformational melt vein in the deepest part of the crust overlying the peridotites, in the Crust-Mantle Extensional Shear Zone (GF13-7; Figs. 2B,C,D); and 2-four leucocratic dikes crosscutting the main foliations either of the Ronda Peridotites or of the overlying deep crust (GF13-6, GF13-9, GF13-10 and GF13-5a) (Figs. 2G,H,I).

Accepted Article

This article is protected by copyright. All rights reserved.

The syn-deformational intrusions (e.g. GF 13-7) are structurally associated with a topto-the-W shearing along the Crust-Mantle Extensional Shear Zone (Frasca et al., 2016;[START_REF] Précigout | Deformation and exhumation of the Ronda Peridotite (Spain)[END_REF]. Syn-deformational intrusions are located either in veins that rotated during shearing (Fig. 2B) or along C′-type shear bands in the high-temperature granulitic foliation (Fig. 2C). C′-type shear bands are slightly oblique to the regional high-temperature foliation, while veins cut the foliation at higher angles (see stereoplots in Figs. 2E,F).

In contrast, the undeformed dikes intruding the Ronda Peridotites and the crustal envelope (GF 13-5a, GF13-6, GF 13-9 and GF 13-10) crosscut the regional foliation (Figs. 2G,H) and show a random distribution in the stereographic projection (Fig. 2L). These dikes are therefore unrelated to the high-temperature foliations (stereoplot in Fig. 2L) and their source is structurally identified below the Ronda Peridotites (Fig. 2A). This second type of intrusion was thus emplaced during the motion of the Ronda Peridotite Thrust.

Petrography and geochemistry

A magmatic origin for the syn-deformational and undeformed sampled intrusions is clearly indicated by interlocking, graphic and granophyric textures (Figs. 2D,I; see also [START_REF] Cuevas | Tectonic implications of the granite dike swarm in the Ronda Peridotites (Betic Cordilleras, southern Spain)[END_REF] and [START_REF] Barich | Microstructures and petrology of melt inclusions in the anatectic sequence of Jubrique (Betic Cordillera, S Spain): Implications for crustal anatexis[END_REF]). U-Pb LA-ICP-MS dating and Nd and Sr isotope analyses were performed at Géosciences Rennes following the procedures described in [START_REF] Ballouard | Tectonic record, magmatic history and hydrothermal alteration in the Hercynian Guérande leucogranite[END_REF] (see Tables DR1 1 andDR2 1 for details and operating conditions of the LA-ICP-MS measurements 1 ).

The isotopic signature of the dikes is similar to those of other Alboran Domain granites [START_REF] Rossetti | Early Miocene strike-slip tectonics and granite emplacement in the Alboran Domain (Rif chain, Morocco): significance for the geodynamic evolution of Western Mediterranean[END_REF] and references therein). The high 87 Sr/ 86 Sr (0.715037, GF13-9, and 0.718605, GF13-7) and the low 143 Nd/ 144 Nd (0.512020, GF13-9, and 0.512026, GF13-7) values indicate a predominantly crustal origin for both types of dikes, compatible with dehydration melting of a metasedimentary source (Table 1). However, the different structural

Accepted Article

This article is protected by copyright. All rights reserved. contexts (Fig. 2A) and the large differences in their rare earth element geochemistry (REE, Fig. 3A) suggest different timings and different sources for the melt, as supported by their respective TDM (i.e. Time Depleted Mantle model age) values of 1270 and 2600 Ma (Tables DR3 1 andDR5 1 ).

In situ melting in the lower crust above the peridotites at ca. 22.5 Ma

The syn-deformational intrusion in a shear zone located inside the granulites (GF13-7) is peraluminous, mainly composed of K-feldspar, quartz and plagioclase, with accessory tourmaline, and presents a concave REE pattern (Fig. 3A). This enrichment in heavy REE suggests that garnet was consumed in the melting reaction (Fig. 3A and Table 1).

Cathodoluminescence (CL) images of the zircon grains (Fig. 3B insets) show distinct cores and rims for some of them, although, when euhedral, the grains are fairly homogeneous.

The cores yield dates ranging from Neoarchaean to Cretaceous (Table DR3 1 ) in good agreement with the TDM value for this sample (Table DR3 1 ). The rims and the homogeneous grains give consistent Miocene dates. In a Tera-Wasserburg diagram, these analyses plot in a concordant to sub-concordant position (Fig. 3B). The concordant data yield a concordia age of 22.52±0.16 Ma (MSWD=1.4, n=16), consistent with a mean 206 Pb/ 238 U age of 22.52±0.18 (MSWD=1.6, n=20).

Undeformed dikes: Melt extraction from the metamorphic sole below the Ronda

Peridotites at ca. 20 Ma

Undeformed dikes intruding the peridotites and the overlying crustal rocks (GF13-6, GF13-5A) have a chondrite-normalized REE pattern enriched in light REE and depleted in heavy REE compared with the extensional-related sample (GF 13-7), suggesting a lower

Accepted Article

This article is protected by copyright. All rights reserved. degree of partial melting and different sources (Fig. 3A). Our samples GF13-6 and GF13-5a present REE patterns similar to those of [START_REF] Rossetti | Early Miocene strike-slip tectonics and granite emplacement in the Alboran Domain (Rif chain, Morocco): significance for the geodynamic evolution of Western Mediterranean[END_REF] from the Rif (Fig. 3A). These rocks also contain slightly more plagioclase and biotite than sample GF13-7, and often titanite.

The dike from which GF13-5a was sampled crosscuts the granulite foliation. Of the zircon grains studied from this sample, three inherited cores give a Palaeozoic date of 337.2±0.13 Ma (MSWD=0.21) (Table DR4 1 ), while the remaining 15 data points are concordant (Fig. 3C) and yield a concordia age of 20.66±0.13 Ma (MSWD=1.06; n=15) with a mean 206 Pb/ 238 U age of 20.65±0.13 Ma (MSWD=0.17; n =15) (Table DR4 1 ). The inherited components are far less abundant than in the extensional-related sample GF13-7 and indicate a different origin, suggesting two distinct partial melting events. Samples GF13-9 and GF13-10 consistently show inherited components similar to sample GF13-5a, suggesting a different origin from sample GF13-7. In sample GF13-9, analyses plot in a concordant to discordant position depending on the degree of Pb loss and/or the presence of common Pb (Fig. 3D). A group of 10 concordant analyses yield a concordia age of 19.98±0.27 Ma (MSWD=2.2; n=10) and a mean 206 Pb/ 238 U age of 20.14±0.25 Ma (MSWD=2.4, n=14) (Table DR5 1 ).

For sample GF13-10, zircon CL images show rare grains with distinct cores and rims.

Only 3 cores were dated and give Palaeozoic dates. The remaining data points give Miocene dates. The most concordant data (Fig. 3E) define a concordia age of 20.18±0.18 Ma (MSWD=1.9; n=18), in agreement with the mean 206 Pb/ 238 U age of 20.08±0.20 Ma (MSWD=2.6; n=23) (Table DR6 1 ).

Accepted Article

This article is protected by copyright. All rights reserved.

AGES OF EXTENSION AND THRUSTING IN SEMI-BRITTLE STRUCTURES

Lower-crustal melting during thinning is coeval with upper-crustal shearing accommodated by C′-type extensional shear planes (Figs. 4A,B and structural details in Frasca et al., (2016)). 40 Ar/ 39 Ar step-heating analysis of the white micas along these C′-type shear planes (Sample B-1505, Fig. 4C) yields a discordant age spectrum due to the mixing of two generations of mica, with a young component dated around 20 Ma and an older one of 24 Ma or older.

Therefore, the intermediate segment at 22.59±0.13 Ma is meaningless (Fig. 4C). The riftrelated structures are then crosscut by thrust-related structures (thrusts and dextral strike-slip faults, see details in [START_REF] Frasca | Structural record of Lower Miocene westward Alboran Domain motion in the Western Betics (southern Spain)[END_REF]). 40 Ar/ 39 Ar step-heating analysis of the white micas extracted from the matrix of the tectonic breccia (Sample B-1503, Fig. 4D) along one of the dextral strike-slip lateral ramps of the Ronda Peridotites Thrust (see stereoplot in Fig. 4B)

gives a plateau age of 20.16 ± 0.11 Ma at the end of degassing, after the release of an excess argon component for 60% of 39 Ar (see the age spectrum in Fig. 4D, complete isotopic results in Table DR7 1 and method details in Table DR1 1 ). Note that the two generations of micas in sample B1505 (rift-related structure, Fig. 4C) may reflect these two stages of deformation, with the old ages at ca. 24 Ma related to rifting and the young ages at ca. 20 Ma related to thrusting.

DISCUSSION AND CONCLUSION

The coupling of geochronological and geochemical data with the regional structures allows us to identify and date two distinct tectonic events during the Lower Miocene, related to mantle exhumation and its final thrust emplacement: 1-lithospheric thinning in a backarc setting leading to in-situ melting at ca. 22.5 Ma, coherent with C′-type extensional shear bands in the mid/upper crust (Ar-Ar ages); and 2-rift-inversion marked by a hot thrusting of

Accepted Article

This article is protected by copyright. All rights reserved.

the Ronda Peridotites and the overlying crustal envelope onto the continental crust, dated at ca. 20 Ma by undeformed dikes cross-cutting the former HT foliation (in both the peridotites and the overlying high-grade crustal rocks) and by brittle thrust-related faults. The differences in geochemistry, age and structural setting unequivocally separate these two different crustal melting events (Fig. 5).

The new age of ca. 22.5 Ma for the rift-related structures is compatible with previous geochronological data attesting to Oligo-Miocene rifting in the Alboran Domain: U-Pb ages of 22.0±0.3Ma and 22.7±0.3Ma, obtained on zircon rims from granulites (Platt et al., 2003a) in Ronda and Beni Bousera (Moroccan equivalent of the Ronda peridotites); LA-ICP-MS ages of 21.3±0.3Ma obtained on syn-foliation monazites within the granulites, and an 40 Ar/ 39 Ar biotite age of 22.5±0.5 Ma from granulites overlying the Beni Bousera Peridotites [START_REF] Frasca | Structural record of Lower Miocene westward Alboran Domain motion in the Western Betics (southern Spain)[END_REF]Michard et al., 2006).

The new ages of ca. 20 Ma for the dikes intruding the thinned lithosphere (GF13-5a, GF13-9 and GF13-10) are relatively younger than suggested by some previous U-Pb data (Fig. 5). The older ages of dikes intruding the peridotites (21.5±3.8 and 22.6±1.8 Ma; [START_REF] Esteban | A revised Aquitanian age for the emplacement of the Ronda peridotites (Betic Cordilleras, southern Spain)[END_REF] suggest a possible onset of thrusting at 22-21 Ma, although the large errors affecting the concordia ages and the limited number of analyzed zircon grains do not permit precise dating. Zircon ages of 21.79±0.6 (younger intersection on the concordia diagram) and 22.3±0.2 Ma (weighted 206 Pb/ 238 U mean age), reported from dikes in the thinned continental crust above the peridotites, have been interpreted as a metamorphic imprint on Hercynian partial melting products or as thrust-related intrusion [START_REF] Sánchez-Navas | Pre-Alpine discordant granitic dikes in the metamorphic core of the Betic Cordillera: tectonic implications[END_REF][START_REF] Rossetti | Timing and modes of granite magmatism in the core of the Alboran Domain, Rif chain, northern Morocco: Implications for the alpine evolution of Western Mediterranean[END_REF]. Ar-Ar data on biotites from undeformed granites [START_REF] Monié | Cooling and exhumation of the Western Betic Cordillera, 40 Ar/ 39 Ar thermochronological constraints on a collapsed terrane[END_REF][START_REF] Sosson | Late exhumation stages of the Alpujarride complex (western Betic Cordilleras, Spain): new thermochronological and structural data on Los Reales and Ojen nappes[END_REF] and the Ar-Ar plateau age on white micas from thrust-related structures presented in this study support thrusting at 21-20 Ma. Finally, a major change in sedimentation in both the external and the internal domains at around 20 Ma is interpreted as

Accepted Article

This article is protected by copyright. All rights reserved. a marker of the onset of thrusting in the Western Betics (see review in [START_REF] Frasca | Structural record of Lower Miocene westward Alboran Domain motion in the Western Betics (southern Spain)[END_REF]).

The compilation of available geochronological data for the Western Alboran domain (western Betics and Rif) shows a bimodal distribution of Ar-Ar ages on micas, zircon fission track (ZFT) and apatite fission track (AFT) ages, suggesting two distinct cooling events related to the two crustal-melting events reported here (Fig. 5, Table DR8 1 ). A first cooling event occurred between 25 and 22-20 Ma and may mark the late stage of continental rifting (Fig. 5, see discussion in Frasca et al. (2016)). However, the majority of ages using mediumto low-temperature thermochronometers are younger than 20 Ma (Fig. 5). The absence of older ages in the thrust footwall (see ages in dark grey in Fig. 5; [START_REF] Monié | Cooling and exhumation of the Western Betic Cordillera, 40 Ar/ 39 Ar thermochronological constraints on a collapsed terrane[END_REF][START_REF] Esteban | Insights into extensional events in the Betic Cordilleras, southern Spain: New fission-track and U-Pb SHRIMP analyses[END_REF] supports a reset of these low-temperature thermochronometers during heating at ca. 21 Ma, as a consequence of the hot emplacement of the peridotites nappe on a colder crustal unit (Fig. 5), followed by a general fast cooling. The characteristic time for cooling (t c ) is controlled by the thickness of the thrust sheet (the Ronda Peridotites and their crustal envelope) and its thermal diffusivity, t c =d 2 /Κ, where K=10 -6 m 2 /s [START_REF] Turcotte | Geodynamics: Applications of continuum physics to geological problems[END_REF]. Taking 10 km as an average thrust-sheet thickness yields a characteristic time of 3.17 Ma for the conductive cooling of the thrust sheet, explaining the regional fast cooling starting at ca. 21 Ma. At the regional scale, the tectonic history of the Ronda Peridotites is the result, during slab rollback, of a fast switch at ca. 21 Ma from exhumation driven by back-arc extension (i.e. rifting) to a thrust emplacement of the overriding plate onto the Iberian margin. This sequence of events, from rifting between 30 and 21 Ma to thrusting at ca. 21 Ma, is consistent with the geodynamic model of [START_REF] Duggen | Magmatic evolution of the Alboran region: The role of subduction in forming the Western Mediterranean and causing the Messinian salinity crisis[END_REF], based on geochemical tracers: subduction processes (marked by the Malaga tholeiite dikes) followed by nappe emplacement (marked by crustal melting). A remaining key issue is to identify the exact location of the rift inversion, either to the north of the Algerian margin or close to the present-day position of the

Accepted Article

This article is protected by copyright. All rights reserved. (G) Undeformed granite intrusion in the peridotites. (H) Dike cross-cutting the regional foliation in the granulite of the lower crust. (I) Granophyric textures of different sizes testifying to the different phases of crystallization of K-feldspar-quartz intergrowths in a melt (GF13-5a); crossed-polarized light. (L) Stereoplot of the dikes and sills cutting through the peridotites and dikes cutting the regional foliation in the lower crust (Carratraca area). GPS locations of the sampled outcrops can be found in Table DR9 1 .

Accepted Article

This article is protected by copyright. All rights reserved. rift-related La Robla normal fault and El Chenil Low Angle Normal Fault (LANF). Blue: riftinversion-related compressional and strike-slip faults (modified after Frasca et al., 2016). For location and details of colours see Fig. 1. GPS locations of the sampled outcrops can be found in Table DR9 1 andDR7 in "supplementary material" 1 ).

Figure 5.

Frequency distribution of the ages available for the western Betics, Rif, ODP Leg 161 and ages obtained in our study (symbols as in Fig. 1). The black line is the best approximation of the trend of the data (polynomial degree 7), excluding the extremity values.

From top to bottom: i) U-Pb on zircons [START_REF] Platt | Early Miocene high-temperature metamorphism and rapid exhumation in the Betic Cordillera (Spain): evidence from U-Pb zircon ages[END_REF][START_REF] Sánchez-Rodríguez | Mesozoic formation of pyroxenites and gabbros in the Ronda area (southern Spain), followed by Early Miocene subduction metamorphism and emplacement into the middle crust: U-Pb sensitive high-resolution ion microprobe dating of zircon[END_REF]Platt et al., 2003b;[START_REF] Esteban | Timing and origin of zircon-bearing chlorite schists in the Ronda peridotites (Betic Cordilleras, Southern Spain)[END_REF][START_REF] Esteban | A revised Aquitanian age for the emplacement of the Ronda peridotites (Betic Cordilleras, southern Spain)[END_REF]; ii) 40 Ar/ 39 Ar on micas [START_REF] Chalouan | The Rif belt[END_REF][START_REF] Monié | Cooling and exhumation of the Western Betic Cordillera, 40 Ar/ 39 Ar thermochronological constraints on a collapsed terrane[END_REF]Kelley and Platt, 1999;[START_REF] Sosson | Late exhumation stages of the Alpujarride complex (western Betic Cordilleras, Spain): new thermochronological and structural data on Los Reales and Ojen nappes[END_REF]Platt et al., 2003a;[START_REF] Rossetti | Timing and modes of granite magmatism in the core of the Alboran Domain, Rif chain, northern Morocco: Implications for the alpine evolution of Western Mediterranean[END_REF]Frasca et al., 2016); and iii) fission tracks on zircon (Platt et al., 2003a;Esteban et al., 2004[START_REF] Esteban | Metamorfismo, exhumación y termocronología de la Unidad de Yunquera (Alpujárrides occidentales, Cordilleras Béticas)[END_REF][START_REF] Esteban | Timing and origin of zircon-bearing chlorite schists in the Ronda peridotites (Betic Cordilleras, Southern Spain)[END_REF]Azdimousa et al., 2013) and apatite (Hurford et al., 1999;[START_REF] Sosson | Late exhumation stages of the Alpujarride complex (western Betic Cordilleras, Spain): new thermochronological and structural data on Los Reales and Ojen nappes[END_REF]Platt et al., 2003a;Esteban et al., 2004[START_REF] Esteban | Metamorfismo, exhumación y termocronología de la Unidad de Yunquera (Alpujárrides occidentales, Cordilleras Béticas)[END_REF]Azdimousa et al., 2013). Pale and dark grey colours in the histograms correspond to samples above and below the RPT, respectively. The bimodal distributions of Ar-Ar, ZFT and AFT ages are well correlated with the two Lower Miocene tectonic events (syn-rift in red and rift inversion in blue; see sketches on top).

Accepted Article

This article is protected by copyright. All rights reserved.

Figure 1 .

 1 Figure 1. (A) Tectonic map of the Ronda Peridotites massifs (from Frasca et al. 2015),

Figure 2 .

 2 Figure 2. (A) Schematic relationships (not to scale) between peridotites, overlying rocks and

Figure 3 .

 3 Figure 3. (A) Different chondrite-normalized patterns for the collected samples and those of

Figure 4 .

 4 Figure 4. (A) 3D visualization of the southeastern termination of the Carratraca massif, with

  . (B) Stereoplots of data collected in the area of (A) (lower hemisphere; Schmidt net): foliations in the upper crustal rocks, stretching lineations and associated C′type shear bands with the arrow indicating the sense of shear (OSXStereonet: version 9.2.0. http://www.ux.uis.no/~nestor/work/programs.html); and rift-inversion-related fault surfaces (great circles) with slip vectors (blue arrows) (FSA software by Célérier (2013): version 35.2). (C,D) 40 Ar/ 39 Ar age spectra for white micas extracted from samples B-1503 (rift-related structure) and B-1505 (thrust-related structure; see method and analytical data in Tables DR1

Table 1 .

 1 Major, trace element and isotopic compositions of the melt products sampled in the Ronda Peridotites and their crustal envelope. ASI = Al 2 O 3 /(Na 2 O-1.67P 2 O 5 +K 2 O+CaO), molar; A/NK = Al 2 O 3 /(Na 2 O+K 2 O), molar; A/2CNK = Al 2 O 3 /(CaO+Na 2 O+K 2 O), molar.

	Sample	GF13-7	GF13-6	GF13-5A
	SiO 2 (wt.%)	73.66	75.99	83.05
	TiO 2	0.14	0.06	0.09
	Al 2 O 3	14.00	13.24	8.19
	FeOt	1.81	0.68	0.72
	MnO	0.02	0.01	0.01
	MgO	0.36	0.48	0.18
	CaO	1.33	1.95	0.61
	Na 2 O	3.61	4.80	1.51
	K 2 O	2.64	0.74	4.31
	P 2 O 5	0.14	0.05	0.03
	LOI	1.57	0.87	0.93
	Total	99.28	98.87	99.63
	ASI	1.27	1.09	1.00
	A/NK	1.59	1.52	1.15
	A/2CNK	1.03	0.84	0.87
	Be (ppm)	< 1	5	2
	Sc	7	6	< 1
	V	16	5	8
	Cr	120	190	150
	Co	4	3	2
	Ni	< 20	< 20	< 20
	Rb	96	20	144
	Sr	405	204	82
	Y	24	13	3
	Zr	66	42	41
	Nb	11	7	4
	Sn	8	5	2
	Cs	1.4	< 0.5	2
	Ba	407	878	271
	La	11.1	13.6	20.2
	Ce	22.8	27.1	42.6
	Pr	2.62	3.03	4.97
	Nd	9.7	10.2	18.2
	Sm	2.8	2.9	4.1
	Eu	1.03	0.85	0.78
	Gd	2.8	2.3	2
	Tb	0.6	0.4	0.2
	Dy	4.2	2.6	0.8
	Ho	0.9	0.5	0.1
	Er	2.6	1.3	0.3
	Tm	0.48	0.21	0.05
	Yb	3.9	1.3	0.2
	Lu	0.62	0.17	0.04
	Hf	1.9	1.8	1.5
	Ta	2.1	1.8	0.3
	W	13	< 1	< 1
	Tl	0.3	< 0.1	0.5
	Pb	23	12	32
	Bi	< 0.4	< 0.4	< 0.4
	Th	3.4	5.8	12.2
	U	4.4	6.6	13.3
	87 Sr/ 86 Sr (standard)	0.718605	0.720831	0.719869
	±2σ	0.000013	0.000012	0.000012
	143 Nd/ 144 Nd (standard)	0.512026	0.512068	0.512140
	±2σ 147 Sm/ 144 Nd	0.000005 0.15528	0.000004 0.14569	0.000006 0.12805
	εNd(t)	-11.26	-10.99

Acknowledgements

This work was funded by the European Union FP7 Marie Curie ITN "TOPOMOD", contract 264517. Thanks to A. Delplanque and P. Frasca for help in drawing Fig. 5 and to Bernard Célérier for useful discussions. Very constructive comments made by David Chew, an anonymous reviewer, the Associate Editor and Editor helped to improve the manuscript.

1 Supplementary material

Accepted Article

This article is protected by copyright. All rights reserved.

Frasca, G., Gueydan, F., Brun, J.-P., Monié, P., 2016, Deformation mechanisms in a continental rift up to mantle exhumation. Field evidence from the western Betics, Spain. Marine and Petroleum Geology, v. 76, p. 1-19. Garrido, C. J., Gueydan, F., Booth-Rea, G., Précigout, J., Hidas, K., Padrón-Navarta, J. A., and Marchesi, C., 2011, Garnet lherzolite and garnet-spinel mylonite in the Ronda peridotite:

Vestiges of Oligocene backarc mantle lithospheric extension in the Western Mediterranean. Geology, v. 39, no. 10, Gueydan, F., Pitra P., Afiri, A., Poujol, M., Essaifi, A., and Paquette, J.-L, 2015, Oligo-Miocene thinning of the Beni Bousera peridotites and their Variscan crustal host rocks, Internal Rif, Morocco. Tectonics, 34, doi: 10.1002/2014TC003769. Hidas, K., Booth-Rea, G., Garrido, C. J., Martínez-Martínez, J. M., Padrón-Navarta, J. A., Konc, Z., Giaconia, F., Frets, E., and Marchesi, C., 2013, Backarc basin inversion and subcontinental mantle emplacement in the crust: kilometre-scale folding and shearing at the base of the proto-Alborán lithospheric mantle (Betic Cordillera, southern Spain). Journal of the Geological Society, v. 170, no. 1, p. 47-55.

Hurford, A. J., Platt, J. P., and Carter, A., 1999, Fission-track analysis of samples from the Alboran Sea basement. In Zahn, R., Comas, M. C. and Klaus, A., editors. Proceedings of the Ocean Drilling Program, Scientific Results, v. 161, no. 21, Kelley, S.P., and Platt, J.P., 1999, Ar-ar dating of biotite and muscovite from Alboran basement samples. In Zahn, R., Comas, M. C., and Klaus, A., editors. Proceedings of the Ocean Drilling Program, Scientific Results, v. 161, no. 22, p. 301-305. Lundeen, M., 1978, Emplacement of the Ronda Peridotite, Sierra Bermeja, Spain. Geological Society of America Bulletin, v. 89, p. 172-180. Marchesi, C., Garrido, C. J., Bosch, D., Bodinier, J.-L., Hidas, K., Padrón-Navarta, J.A., and Gervilla, F., 2012, A Late Oligocene Suprasubduction Setting in the Westernmost

Accepted Article

This article is protected by copyright. All rights reserved.

Mediterranean Revealed by Intrusive Pyroxenite Dikes in the Ronda Peridotite (southern Spain). The Journal of Geology, v. 120, no. 2,

Accepted Article

This article is protected by copyright. All rights reserved.

Gibraltar arc. These two hypotheses would imply either a large (e.g., Hidas et al., 2013) or a limited (e.g., [START_REF] Frasca | Structural record of Lower Miocene westward Alboran Domain motion in the Western Betics (southern Spain)[END_REF] amount of rollback since 20 Ma, respectively.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.


Kevin Schmidmayer
email: kevin.schmidmayer@univ-amu.fr

Fabien Petitpas
email: fabien.petitpas@univ-amu.fr

Eric Daniel
email: eric.daniel@univ-amu.fr

Nicolas Favrie
email: nicolas.favrie@univ-amu.fr

Sergey Gavrilyuk
email: sergey.gavrilyuk@univ-amu.fr

A

Keywords: diffuse interface, Godunov type methods, hyperbolic systems, multiphase flows, shock waves, surface tension

Introduction

The breakup of liquid droplets induced by high speed flows has a wide range of engineering and scientific applications and has given rise to a large number of publications. In some cases, this phenomenon causes damages as for example when droplets are impacting aircrafts in supersonic flight causing erosion of its surface (Engel [8], Joseph et al. [START_REF] Joseph | Breakup of a liquid drop suddenly exposed to a high-speed airstream[END_REF], Igra and Takayama [START_REF] Igra | Investigation of aerodynamic breakup of a cylindrical water droplet[END_REF][START_REF] Igra | Numerical simulation of shock wave interaction with a water column[END_REF]). Studying of droplets behavior in a high speed flow may also be encountered when security issues are considered as, for example, for shock wave attenuation (Chauvin et al. [START_REF] Chauvin | Shock waves in sprays: numerical study of secondary atomization and experimental comparison[END_REF][START_REF] Chauvin | Experimental investigation of the propagation of a planar shock wave through a two-phase gas-liquid medium[END_REF]). Other applications can be found in explosive science or in combustion systems where a liquid jet is atomized (Welch and Boyle [START_REF] Welch | New turbines to enable efficient geothermal power plants[END_REF], Meng and Colonius [START_REF] Meng | Numerical simulations of the early stages of high-speed droplet breakup[END_REF], Devassy et al. [START_REF] Devassy | Atomization modelling of liquid jets using a two-surface-density approach[END_REF]). Detailed reviews on droplet breakup can be found in Pilch and Erdman [START_REF] Pilch | Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop[END_REF], Wierzba and Takayama [START_REF] Wierzba | Experimental investigation on liquid droplet breakup in a gas stream[END_REF], Hsiang and Faeth [START_REF] Hsiang | Near-limit drop deformation and secondary breakup[END_REF].

Concerning numerical simulations, the breakup study is usually focused on the first stages of droplet deformation when Richtmyer-Meshkov and/or Rayleigh-Taylor instabilities appear (Yang et al. [START_REF] Yang | A model for characterization of a vortex pair formed by shock passage over a light-gas inhomogeneity[END_REF], Quirk and Karni [START_REF] Quirk | On the dynamics of a shock-bubble interaction[END_REF], Layes and Le Metayer [START_REF] Layes | Quantitative numerical and experimental studies of the shock accelerated heterogeneous bubbles motion[END_REF], Meng and Colonius [START_REF] Meng | Numerical simulations of the early stages of high-speed droplet breakup[END_REF]), but not on the further stages when capillary and/or viscous effects become significant.

In the last decades, several theoretical studies have been performed to treat capillary effects in multiphase flows. The seminal work of Brackbill et al. [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF] succeeded in transforming a surface force into a volume force, quite easy to treat as a source term in a multiphase flow model. The surface tension volume force is expressed thanks to a color function c (x). This approach has been used in Chen and Doolen [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF], Sussman et al. [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF], Gueyffier et al. [START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF], Osher and Fedkiw [START_REF] Osher | Level set methods: An overview and some recent results[END_REF][START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF], Tryggvason et al. [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF], Périgaud and Saurel [START_REF] Périgaud | A compressible flow model with capillary effects[END_REF], Le Martelot et al. [START_REF] Le Martelot | Towards the direct numerical simulation of nucleate boiling flows[END_REF] where capillary effects are added into the momentum and the energy equations.

The aim of this work is to develop a mathematical model for fluid flows with capillary effects that is hyperbolic, verifies conservation principles and entropy inequality together with a suitable numerical method capable to treat the effect of the flow on the droplet from the short time scale when the shock wave interacts with the droplet to the long time scale when capillary effects become significant. We focus in this study on multiphase compressible fluid flows only. Viscous and heat conduction are not taken into account and will be a part of future works. Some ideas on the treatment of heat conduction in multiphase compressible flows can already be found in [START_REF] Petitpas | A discrete method to treat heat conduction in compressible two-phase flows[END_REF].

Section 2 presents the Brackbill et al. [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF] method to treat the surface tension and a review of existing models with a conservative form of the capillary terms. In Section 3, the new model with capillary effects is presented. The model is in agreement with the conservation principles and with the second law of thermodynamics. It is shown that the model is weakly hyperbolic. It has two sound characteristics associated with the classical compression waves and two new sound characteristics associated with the capillary effects. However, for multiple contact characteristics one eigenvector is always missing. Section 4 is devoted to the building of a numerical method able to solve capillary terms in a conservative manner. The method is based on split models that are separately hyperbolic or weakly hyperbolic. These submodels are solved thanks to adequate numerical schemes. Section 5 presents the validation of the method on 2D test cases. It shows that the model and the numerical method are able to treat accurately both capillary effects and shock wave propagation. Quantitative comparisons are done with other methods based on source terms integration to show the importance of the conservative formulation. To illustrate the capabilities of the model, the aerodynamic breakup of a water column induced by a shock wave is numerically solved and is compared with experiments. In Appendix, the model derivation is given.

Compressible two-phase capillary flows: state of the art

Surface tension force and color function

The main difficulty in modeling the capillary effects is about considering a surface force in numerical models that solve volume average quantities. The seminal work of Brackbill et al. [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF], called CSF (Continuum Surface Force) method, succeeded to do it by using a color function, c (x). Thanks to this function, the surface tension volume force is then expressed:

F v (x) = σκ (x) ∇c (x) [c] ,
where σ is the surface tension coefficient and κ (x) the local curvature of the interface defined by: κ (x) = -∇ • n (x),

where n (x) is the normal vector to the interface between the both phases:

n (x) = ∇c (x) ∇c (x) .
The color function c (x) allows locations of the different fluids and the interface. c (x) is defined as:

c (x) =    c 1 in fluid 1, c 2
in fluid 2, c 1 c (x) c 2 in the transition region.

(1)

In the transition region c (x) is given by interpolation, meaning that the interface has a non zero thickness. [c] = c 2 -c 1 is the jump of the color function.

It is assumed that the color function obeys a transport equation [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF]:

∂c (x) ∂t + u I • ∇c (x) = 0,
where u I is the interface velocity. Numerical results using this force can be found in [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF][START_REF] Le Martelot | Towards the direct numerical simulation of nucleate boiling flows[END_REF][START_REF] Osher | Level set methods: An overview and some recent results[END_REF][START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF][START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF][START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF]. In these references, the surface tension force is treated as source terms in the momentum and the energy equations. Nevertheless, this treatment of capillary effects violates conservation principles.

Review of existing compressible models with capillary effects

Two family of methods are available to treat interface problems.

• The first family of methods considers interfaces as sharp. Sharp interfaces can be obtained using interface-tracking methods, where usually a level set function tracks the interface (Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF]). However, such formulations often involve slight modifications of the governing equations. For example a pressure evolution equation can replace the energy equation around interfaces (Karni [25,[START_REF] Karni | Hybrid multifluid algorithms[END_REF]). In the Ghost Fluid Method (Fedkiw et al. [START_REF] Fedkiw | A non oscillatory Eulerian approach to interfaces in multimaterial flows (The Ghost Fluid Method)[END_REF]) and its simplified version (Koren et al. [START_REF] Koren | Riemann-problem and level-set approaches for homentropic two-fluid computations[END_REF]), thermodynamically similar variables are added across interfaces to complete stencils. These methods do not generate spurious oscillations at interfaces, however, they are basically not conservative (Liu et al. [START_REF] Liu | Ghost fluid method for strong shock impacting on material interface[END_REF]). Even if progress has been recently done concerning this aspect (Hu et al. [START_REF] Hu | A conservative interface method for compressible flows[END_REF], Luo et al. [START_REF] Luo | A conservative sharp interface method for incompressible multiphase flows[END_REF], Han et al. [START_REF] Han | Scale separation for multi-scale modeling of free-surface and two-phase flows with the conservative sharp interface method[END_REF], Schranner et al. [START_REF] Schranner | On the convergence of the weakly compressible sharp-interface method for two-phase flows[END_REF]), that makes them less desirable for problems where shock waves are involved.

• The second family of methods is described below and is called diffuse interface methods.

In this class of methods, interfaces are not explicitly tracked but allowed to diffuse numerically (Abgrall and Karni [START_REF] Abgrall | Computations of compressible mutlifluids[END_REF], Saurel and Abgrall [START_REF] Saurel | A simple method for compressible multifluid flows[END_REF]). These methods are particularly interesting because they are able to deal with dynamic appearance and disappearance of interfaces. Moreover, this is also the only class of models where the thermodynamics of mixture cells is well defined, thanks to a specific equation of state for each phase (liquid or gas).

The study of capillary effects within the framework of the diffuse interface methods is based on the generalization of the Allaire et al. [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF] model. This model originally does not include capillary effects. The second law of thermodynamics is verified only if the condition of thermal equilibrium T = T k , with k = {1, 2}, is retained. Périgaud and Saurel [START_REF] Périgaud | A compressible flow model with capillary effects[END_REF] extended this model by including the capillary effects. As a result, the surface tension volume force appears as a flux term in the momentum equation as well as the work of this force F v (x) • u in the total energy equation. In this reference, as Gueyffier et al. [START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF] did within the incompressible flows framework, a conservative formulation was obtained:

           ∂α 1 ∂t + u • ∇α 1 = 0, ∂α k ρ k ∂t + ∇ • (α k ρ k u) = 0, ∂ρu ∂t + ∇ • ρu ⊗ u + P I -σ ∇α 1 I -∇α 1 ⊗∇α 1 ∇α 1 = 0, ∂ρE+εσ ∂t + ∇ • u (ρE + ε σ + P ) -σ ∇α 1 I -∇α 1 ⊗∇α 1 ∇α 1 • u = 0,
where α k and ρ k are the volume fraction and the density of phase k. ρ, u, P , E = e + 1 2 u 2 and e are respectively the mixture variables for density, velocity, pressure, total energy and internal energy. The specific mixture internal energy is defined as e = k Y k e k (ρ k , P ) and each fluid is governed by its own equation of state (EOS) e k = e k (ρ k , P ). Introducing the capillary effects in fluxes (conservative formulation) leads to a new term in the mixture total energy equation, ε σ = σ ∇α 1 which is defined as a capillary potential energy. Note that α 1 obeys a transport equation and plays the role of the color function (1) because u I = u when the velocity equilibrium is considered. In [START_REF] Périgaud | A compressible flow model with capillary effects[END_REF], the authors used a pressure equilibrium closure, allowing the resolution of interface problems dedicated to high speed flows. Nevertheless, this closure does not allow to respect the second law of thermodynamics for the mixture.

Le Martelot et al. [START_REF] Le Martelot | Towards the direct numerical simulation of nucleate boiling flows[END_REF] corrected this drawback by including the capillary effects in a velocity, pressure and temperature equilibrium system of equations by considering the mass fraction as the color function: the second law of thermodynamics is then respected. However, the characteristic velocities were estimated under assumption that the local curvature is fixed. Even if such an approach is physically reasonable, it does not allow us to rigorously estimate the propagation speeds. Also, the thermal equilibrium closure is reasonable for the description of the boiling phenomenon, but it is too drastic for droplet breakup study under high speed flows.

It becomes then obvious that the introduction of a new model having clearly defined sound speeds with the capillary effects in conservative form that can still be compatible with the mixture entropy increase without the temperature equilibrium assumption is necessary.

Mechanical equilibrium model with capillary effects

The model

The new model we propose is reminiscent that of the Kapila et al. model [START_REF] Kapila | Two-phase modeling of DDT in granular materials: Reduced equations[END_REF]. This last model has been shown to be suitable to treat interface problems between compressible fluids. Multiple extensions of this model have been developped to solve several concrete problems dealing with phase transition (Massoni et al. [START_REF] Massoni | Proposition de methodes et modeles Euleriens pour les problemes a interfaces entre fluides compressibles en presence de transfert de chaleur[END_REF], Saurel et al. [START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF], cavitation (Petitpas et al. [START_REF] Petitpas | Diffuse interface models for high speed cavitating underwater systems[END_REF]), detonation in high energetic materials (Petitpas et al. [START_REF] Petitpas | Modelling detonation waves in condensed energetic materials: Multiphase CJ conditions and multidimensional computations[END_REF]), solid-fluid interaction and compaction of granular media (Favrie and Gavrilyuk [START_REF] Favrie | Dynamic compaction of granular materials[END_REF][START_REF] Favrie | Diffuse interface model for compressible fluidcompressible elastic-plastic solid interaction[END_REF]) and low Mach number flows (Murrone and Guillard [START_REF] Murrone | Behavior of upwind scheme in the low mach number limit: Iii. preconditioned dissipation for a five equation two phase model[END_REF]). It considers compressible two-phase flows in mechanical equilibrium (pressure and velocity equilibrium). In this model, the volume fraction obeys the following equation coming from the pressure equilibrium condition:

∂α 1 ∂t + u • ∇α 1 = K∇ • u, (2)
where the term K∇ • u accounts for the differences in the acoustic behaviour of both phases. K is given by:

K = ρ 2 a 2 2 -ρ 1 a 2 1 ρ 2 a 2 2 α 2 + ρ 1 a 2 1 α 1 , a k being the speed of sound of phase k.
The capillary effects are then added in conservative form. Because the color function is a purely geometric variable, a supplementary equation for the color function is added. It will be shown in the following that the introduction of this supplementary equation is an important ingredient to perform the hyperbolicity study. Applying the Hamilton principle for the model derivation (see Appendix for details), the model becomes:

               ∂α 1 ∂t + u • ∇α 1 -K∇ • u = 0, ∂α k ρ k ∂t + ∇ • (α k ρ k u) = 0, ∂ρu ∂t + ∇ • ρu ⊗ u + P I + Ω = 0, ∂ρE+εσ ∂t + ∇ • (ρE + ε σ + P ) u + Ω • u = 0, ∂c ∂t + u • ∇c = 0, (3)
with Ω being the capillary tensor given by:

Ω = -σ ∇c I - ∇c ⊗ ∇c ∇c . (4)
The capillary energy is equal to ε σ = σ ∇c and the color function is normalized by its jump:

c = c/[c].
The surface tension terms do not affect the pressure and the entropy equations. The pressure equation evolution then reads:

dP dt + ρa 2 ∇ • (u) = 0.
where a is the Wood mixture speed of sound [START_REF] Wood | A textbook of sound[END_REF]:

a 2 = ρ k α k ρ k a 2 k -1 . (5)
As in Kapila et al. [START_REF] Kapila | Two-phase modeling of DDT in granular materials: Reduced equations[END_REF], the entropy equations remain unchanged in continuous motion:

ds k dt = 0, with the material derivative operator d(•)/dt = ∂(•)/∂t + u • ∇(•).
The mixture entropy is equation obviously assured:

ds dt = k dY k s k dt = 0.

Hyperbolicity

In this section, the study of the hyperbolicity of system (3) is done. Thanks to the rotational invariance of the equations, the study can be reduced to the study of only 1D equation.

Primitive form

The model needs to be transformed into a vector form:

∂W ∂t + A (W) ∂W ∂x = 0. ( 6
)
The color function is first rewritten by taking its gradient to obtain a conservative equation for w = ∇c:

∂w ∂t + ∇ (u • w) = 0.
Or, in an equivalent form:

∂w ∂t + ∂w ∂x T • u + ∂u ∂x T • w = 0. ( 7
)
Since w is a gradient, we have:

curl (w) = 0, (8)
i.e. ∂w ∂x

T = ∂w ∂x .
Thus, equation [START_REF] Devassy | Atomization modelling of liquid jets using a two-surface-density approach[END_REF] with the constraint (8) becomes:

∂w ∂t + ∂w ∂x • u + ∂u ∂x T • w = 0. (9)
Constraint (8) also appears in solid mechanics where the formulation of the extended system under a "curl" constraint is important in the hyperbolicity study (see Ndanou et al. [START_REF] Ndanou | Criterion of hyperbolicity in hyperelasticity in the case of the stored energy in separable form[END_REF]). Equation ( 9) may be developed in three dimensions, where we denote u = (u, v, w) T , w = (w 1 , w 2 , w 3 ) T and x = (x, y, z) T :

     ∂w 1 ∂t + ∂w 1 ∂x u + ∂w 1 ∂y v
In the case where all the variables depend only on (t, x), (10) reduces to:

   ∂w 1 ∂t + ∂w 1 ∂x u + ∂u ∂x w 1 + ∂v ∂x w 2 + ∂w ∂x w 3 = 0, ∂w 2 ∂t + ∂w 2 ∂x u = 0, ∂w 3 ∂t + ∂w 3 ∂x u = 0.
The vector w is introduced to rewrite components for the capillary tensor Ω (4) in the x-direction:

Ω 11 = σ -w 2 2 -w 2 3 w 2 1 + w 2 2 + w 2 3
,

Ω 12 = σ w 1 w 2 w 2 1 + w 2 2 + w 2 3
,

Ω 13 = σ w 1 w 3 w 2 1 + w 2 2 + w 2 3 .
Finally, the system rewritten in vector form reads for two phases:

                                             ∂α 1 ∂t + u ∂α 1 ∂x -K ∂u ∂x = 0,
∂w 3 ∂x = 0, ∂P ∂t + u ∂P ∂x + ρa 2 ∂u ∂x = 0, ∂w 1 ∂t + u ∂w 1 ∂x + w 1 ∂u ∂x + w 2 ∂v ∂x + w 3 ∂w ∂x = 0, ∂w 2 ∂t + u ∂w 2 ∂x = 0, ∂w 3 ∂t + u ∂w 3 ∂x = 0, ∂s 1 ∂t + u ∂s 1 ∂x = 0, ∂s 2 ∂t + u ∂s 2 ∂x = 0, ∂Y 1 ∂t + u ∂Y 1 ∂x = 0, ∂c ∂t + u ∂c ∂x = 0.

Eigenvalues

The vector W and the matrix A in ( 6) are then defined by:

W = [α 1 , u, v, w, P, w 1 , w 2 , w 3 , s 1 , s 2 , Y 1 , c] T , A (W) =                      u -K 0 0 0 0 0 0 0 0 0 0 u 0 0 1 ρ 1 ρ ∂Ω 11 ∂w 1 1 ρ ∂Ω 11 ∂w 2 1 ρ ∂Ω 11 ∂w 3 0 0 0 0 0 u 0 0 1 ρ ∂Ω 12 ∂w 1 1 ρ ∂Ω 12 ∂w 2 1 ρ ∂Ω 12 ∂w 3 0 0 0 0 0 0 u 0 1 ρ ∂Ω 13 ∂w 1 1 ρ ∂Ω 13 ∂w 2 1 ρ ∂Ω 13 ∂w 3 0 0 0 0 ρa 2 0 0 u 0 0 0 0 0 0 0 w 1 w 2 w 3 0 u 0 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 0 u                     
The determinant of matrix A can be calculated:

det A (W) -λI = (u -λ) 8 (u -λ) 4 + L(u -λ) 2 + M ,
where:

L = -a 2 - 1 ρ w 1 ∂Ω 11 ∂w 1 + w 2 ∂Ω 12 ∂w 1 + w 3 ∂Ω 13 ∂w 1 , M = a 2 ρ w 2 ∂Ω 12 ∂w 1 + w 3 ∂Ω 13 ∂w 1 .
8 real eigenvalues are straightforward found:

λ 1,2,3,4,5,6,7,8 = u.
The 4 other eigenvalues are the roots of the quadratic equation:

X 2 + LX + M = 0, (11)
with X = (u -λ) 2 . If the discriminant ∆ is positive, the roots of equation ( 11) are real and complex numbers otherwise. Defining the components of the normal vector n k = w k / w , after some calculations the discriminant is:

∆ = a 2 + σ w ρ n 2 2 + n 2 3 2 -4a 2 σ w ρ n 2 2 + n 2 .
But: Hence:

n 2 2 + n 2
∆ a 2 + σ w ρ n 2 2 + n 2 3 2 -4a 2 σ w ρ n 2 2 + n 2 3 = a 2 - σ w ρ n 2 2 + n 2 3 2 0.
The roots of the quadratic equation ( 11) are then all real. Thus, the corresponding eigenvalues of model ( 3) are real, explicit and given by:

λ 9,10 = u ± a s , (12)
λ 11,12 = u ± a c , (13)
where:

a 2 s = a 2 + b + (a 2 + b) 2 -4a 2 b (n 2 2 + n 2 3 ) 2 , a 2 c = a 2 + b -(a 2 + b) 2 -4a 2 b (n 2 2 + n 2 3 ) 2 ,
a is the Wood mixture speed of sound previously defined ( 5) and:

b = σ w ρ n 2 2 + n 2 3 .
Finally, there are 8 multiple eigenvalues λ = u and 4 eigenvalues corresponding to sound waves [START_REF] Fedkiw | A non oscillatory Eulerian approach to interfaces in multimaterial flows (The Ghost Fluid Method)[END_REF] and capillary waves (13) (see Figure 1).

The system is hyperbolic if the multiple eigenvalues λ = u have exactly 8 linearly independent eigenvectors. One can prove that it is not the case, one eigenvector is missing (see Appendix), so the system is only weakly hyperbolic. The whole system of the eigenvectors will not be given because this equilibrium system will not be numerically solved in this form (see Section 4 for discussion of the numerical method).

Numerical resolution of model (3)

Basic ideas

The numerical resolution of model ( 3) represents a challenge regarding the two following points :

• The first difficulty is due to the K∇ • u term in the volume fraction evolution equation. This is the precious ingredient leading to the respect of the mixture entropy equation when a pressure equilibrium assumption is retained. Nevertheless, the presence of this non-conservative term considerably complicates the numerical method which crucially depends on the choice of appropriate Rankine-Hugoniot relations (see [START_REF] Gavrilyuk | Rankine-Hugoniot relations for shocks in heterogeneous mixtures[END_REF] and [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] for details). This is a reason why a non-equilibrium pressure model is preferred with a pressure relaxation term instead of the non-conservative term K∇ • u in the volume fraction equation. This model is presented below.

• The second difficulty is in the simultaneous treatment of 5 waves which are present in the model. This difficulty is circumvented by the use of split models that are proven to be hyperbolic and weakly hyperbolic respectively.

Pressure relaxation model with capillary effects

Because the non-conservative volume fraction equation presents a major problem regarding numerical resolution, the following weakly hyperbolic non-equilibrium pressure model is proposed:

                     ∂α 1 ∂t + u • ∇α 1 = µ (P 1 -P 2 ) , ∂α 1 ρ 1 ∂t + ∇ • (α 1 ρ 1 u) = 0, ∂α 2 ρ 2 ∂t + ∇ • (α 2 ρ 2 u) = 0, ∂ρu ∂t + ∇ • ρu ⊗ u + P I + Ω = 0, ∂α 1 ρ 1 e 1 ∂t + ∇ • (α 1 ρ 1 e 1 u) + α 1 P 1 ∇ • u = -µP I (P 1 -P 2 ) , ∂α 2 ρ 2 e 2 ∂t + ∇ • (α 2 ρ 2 e 2 u) + α 2 P 2 ∇ • u = µP I (P 1 -P 2 ) , ∂c ∂t + u • ∇c = 0, ( 14
)
where µ is the pressure relaxation coefficient,

P I = Z 2 P 1 +Z 1 P 2 Z 1 +Z 2
(see [START_REF] Saurel | A multiphase model with internal degrees of freedom: Application to shock-bubble interaction[END_REF] for details) and Z k = ρ k a k is the acoustic impedance of the phase k. The mixture pressure is given by:

P = α 1 P 1 + α 2 P 2 .
Due to the condition P 1 = P 2 in this model, the total energy equation of the mixture is replaced by the internal energy equation for each phase. Nevertheless, the mixture total energy equation of the system can be written in usual form:

∂ρE + ε σ ∂t + ∇ • (ρE + ε σ + P ) u + Ω • u = 0. ( 15
)
The equation ( 15) is redundant when both phasic internal energy equations are solved, but it will appear to be an important ingredient for numerical method to ensure the energy conservation and to preserve a correct treatment of shock waves.

Following Section 3.2, one can prove that model ( 14) is also weakly hyperbolic.

One can note that the surface tension effects are missing in the phasic energy equations since it is only a mixture characteristic.

The entropy equations for system [START_REF] Gavrilyuk | Rankine-Hugoniot relations for shocks in heterogeneous mixtures[END_REF] can also be written:

α 1 ρ 1 T 1 ds 1 dt = µ(P 1 -P 2 ) 2 Z 1 Z 1 + Z 2 , α 2 ρ 2 T 2 ds 2 dt = µ(P 1 -P 2 ) 2 Z 2 Z 1 + Z 2 ,
that insures the mixture entropy (s = Y 1 s 1 + Y 2 s 2 ) increases. Model (3) can be recovered as the asymptotic limit of Model ( 14) when the pressure of each phase tends to be equal. A special splitting procedure will be done for the numerical resolution of model ( 14).

Splitting procedure

Model ( 14) without the relaxation terms is split in two submodels. The first submodel does not take into account the surface tension terms. The hyperbolicity is then easily verified. The second submodel contains the only capillary terms and is proven weakly hyperbolic. Such an approach was first proposed in Favrie et al. [START_REF] Favrie | A thermodynamically compatible splitting procedure in hyperelasticity[END_REF] where the splitting procedure was used to separate the treatment of longitudinal and shear waves in hyperelasticity. The submodels are presented below only in the x-direction.

Hyperbolic submodel 1

The first submodel is similar to that presented in [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] with additional decoupled equations for the gradient of the color function:

                                                                   ∂α 1 ∂t + u ∂α 1 ∂x = 0, ∂α 1 ρ 1 ∂t + ∂α 1 ρ 1 u ∂x = 0, ∂α 2 ρ 2 ∂t + ∂α 2 ρ 2 u ∂x = 0, ∂ρu ∂t + ∂ρu 2 + α 1 P 1 + α 2 P 2 ∂x = 0, ∂ρv ∂t + ∂ρuv ∂x = 0, ∂ρw ∂t + ∂ρuw ∂x = 0, ∂α 1 ρ 1 e 1 ∂t + ∂α 1 ρ 1 e 1 u ∂x + α 1 P 1 ∂u ∂x = 0, ∂α 2 ρ 2 e 2 ∂t + ∂α 2 ρ 2 e 2 u ∂x + α 2 P 2 ∂u ∂x = 0, ∂w 1 ∂t + ∂w 1 u ∂x = 0, ∂w 2 ∂t + u ∂w 2 ∂x = 0, ∂w 3 ∂t + u ∂w 3 ∂x = 0. ( 16
)
This system describes only the transport and the compression waves. The equation for w 1 is taken in conservative form to let the possibility to consider weak solutions. The other terms in this equation will be treated in the second submodel.

Model ( 16) may be rewritten in a vector form [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF] with:

W = [α 1 , s 1 , s 2 , u, v, w, P 1 , P 2 , w 1 , w 2 , w 3 ] T A (W) =                   u 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 P 1 -P 2 ρ 0 0 u 0 0 α 1 ρ α 2 ρ 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 ρ 1 a 2 1 0 0 u 0 0 0 0 0 0 ρ 2 a 2 2 0 0 0 u 0 0 0 0 0 w 1 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 u                  
The eigenvalues of the system are:

λ 1,2,3,4,5,6,7,8,9 = u, λ 10 = u -a f , λ 11 = u + a f ,
where a f is the frozen mixture sound speed:

a 2 f = Y 1 a 2 1 + Y 2 a 2 2 .
The hyperbolicity of this first submodel is proven in [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF].

Weakly hyperbolic submodel 2

The second submodel is:

                                         ∂α 1 ∂t = 0, ∂α 1 ρ 1 ∂t = 0,
∂w 3 ∂x = 0, ∂α 1 ρ 1 e 1 ∂t = 0, ∂α 2 ρ 2 e 2 ∂t = 0, ∂w 1 ∂t + w 2 ∂v ∂x + w 3 ∂w ∂x = 0, ∂w 2 ∂t = 0, ∂w 3 ∂t = 0. ( 17
)
This second system describes the capillary effects. Also, the non-conservative product in the equation for w 1 is well defined because w 2 and w 3 are continuous through the shock. Model (17) may be rewritten in a vector form [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF] with:

W = [α 1 , s 1 , s 2 , u, v, w, P 1 , P 2 , w 1 , w 2 , w 3 ] T A (W) =                    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 ρ ∂Ω 11 ∂w 1 1 ρ ∂Ω 11 ∂w 2 1 ρ ∂Ω 11 ∂w 3 0 0 0 0 0 0 0 0 1 ρ ∂Ω 12 ∂w 1 1 ρ ∂Ω 12 ∂w 2 1 ρ ∂Ω 12 ∂w 3 0 0 0 0 0 0 0 0 1 ρ ∂Ω 13 ∂w 1 1 ρ ∂Ω 13 ∂w 2 1 ρ ∂Ω 13 ∂w 3
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 w 2 w 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

                  
The eigenvalues of the system are:

λ 1,2,3,4,5,6,7,8,9 = 0,

λ 10 = -n 2 2 + n 2 3 σ w ρ , λ 11 = n 2 2 + n 2 3 σ w ρ .
The eigenvalues of the system are all real but, as in the case of model ( 3), there are not as much linearly independent eigenvectors as eigenvalues: one eigenvector is always missing. So, the second submodel is weakly hyperbolic.

Numerical Method

Finally, the numerical method is presented as a 3-step method. Each step is successively performed in order to circumvent specific numerical problems:

• First, the hyperbolic non-equilibrium pressure model ( 16) is solved using a Godunovtype method.

• Second, model ( 17) is solved. A specific attention is paid to the choice for the flux terms in order to ensure the momentum and energy conservation.

• Third, a relaxation procedure leads to the pressure equilibrium.

The chain of these three steps is equivalent to solve model [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF]. Each step of the method is presented in details hereafter.

The full system of equations is first rewritten in the following vector form:

∂U ∂t + ∂ (F x h (U) + F x c (U)) ∂x + ∂ (F y h (U) + F y c (U)) ∂y + ∂ (F z h (U) + F z c (U)) ∂z + H nc ∇ • u = H relax .
The vector U contains the unknown quantities defined in the system:

U = [α 1 , α 1 ρ 1 , α 2 ρ 2 , ρu, ρv, ρw, α 1 ρ 1 e 1 , α 2 ρ 2 e 2 , c, ρE + ε σ ] T
The vectors F β h (U), F β c (U), H nc and H relax , with β = {x, y, z}, contain respectively the hydrodynamic fluxes, the capillary fluxes, the non-conservative terms and the relaxation terms:

F x h (U) =                 α 1 u α 1 ρ 1 u α 2 ρ 2 u ρu 2 + P ρuv ρuw α 1 ρ 1 e 1 u α 2 ρ 2 e 2 u cu (ρE + P ) u                 F y h (U) =                 α 1 v α 1 ρ 1 v α 2 ρ 2 v ρuv ρv 2 + P ρvw α 1 ρ 1 e 1 v α 2 ρ 2 e 2 v cv (ρE + P ) v                 F z h (U) =                 α 1 w α 1 ρ 1 w α 2 ρ 2 w ρuw ρvw ρw 2 + P α 1 ρ 1 e 1 w α 2 ρ 2 e 2 w cw (ρE + P ) w                 F x c (U) =                 0 0 0 Ω 11 Ω 12 Ω 13 0 0 0 ε σ u + Ω 11 u + Ω 12 v + Ω 13 w                 F y c (U) =                 0 0 0 Ω Ω Ω 0 0 0 ε σ u + Ω 21 u + Ω 22 v + Ω 23 w                 F z c (U) =                 0 0 0 Ω 31 Ω 32 Ω 33 0 0 0 ε σ u + Ω 31 u + Ω 32 v + Ω 33 w                 H nc =                 -α 1 0 0 0 0 0 α 1 P 1 α 2 P 2 -c 0                 H relax =                 µ (P 1 -P 2 ) 0 0 0 0 0 -µP I (P 1 -P 2 ) µP I (P 1 -P 2 ) 0 0                
It is important to note that the additional equation for the mixture total energy has been added to the system for numerical solution purposes. This equation is obviously in agreement with the complete system [START_REF] Gavrilyuk | Rankine-Hugoniot relations for shocks in heterogeneous mixtures[END_REF] and will be necessary to correct the energy conservation equation during the final relaxation step.

The unknown vector U n+1 is obtained from the initial condition U n by application of the three successive operators according to the sequence:

U n+1 = L relax L cap L hyper (U n ) .
Each step of the numerical method corresponds to the application of one of the three operators detailed below in a cartesian 2D framework.

Hyperbolic operator

The application of the first operator L hyper corresponds to the resolution of the hyperbolic submodel (16) using a Godunov-type method [START_REF] Godunov | A finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF] extended to an high-order scheme with a MUSCL-Hancock procedure. The solution for this step is given for the cell (i, j) by:

U hyper i,j = U n i,j -G∆t, G =                1 ∆x F * h;i+ 1 2 ,j U n+ 1 2 i,j,R , U n+ 1 2 i+1,j,L -F * h;i-1 2 ,j U n+ 1 2 i-1,j,R , U n+ 1 2 i,j,L + 1 ∆y F * h;i,j+ 1 2 U n+ 1 2 i,j,T , U n+ 1 2 i,j+1,B -F * h;i,j-1 2 U n+ 1 2 i,j-1,T , U n+ 1 2 i,j,B + H n nc;i,j ∆x u * i+ 1 2 ,j U n+ 1 2 i,j,R , U n+ 1 2 i+1,j,L -u * i-1 2 ,j U n+ 1 2 i-1,j,R , U n+ 1 2 i,j,L + H n nc;i,j ∆y v * i,j+ 1 2 U n+ 1 2 i,j,T , U n+ 1 2 i,j+1,B -v * i,j-1 2 U n+ 1 2 i,j-1,T , U n+ 1 2 i,j,B               
.

The superscript * represents the solution of the Riemann problem on the corresponding cell boundary using the extrapolated values to the cell boundary U n+ 1 2 in the case of high order method (the subscripts R, L, T , B represent respectively the right, left, top and bottom neighbouring cells). These fluxes can be computed by any Riemann solver. Here a Harten-Lax-van Leer Contact (HLLC) approximate Riemann solver is used. Details of this method can be found in Saurel et al. [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF]. Since the time marching scheme is explicit, the time step obeys a classical Courant-Friedrichs-Lewy (CFL) criterion.

Capillary operator

The second operator L cap introduces capillary terms and corresponds to the resolution of submodel [START_REF] Han | Scale separation for multi-scale modeling of free-surface and two-phase flows with the conservative sharp interface method[END_REF]:

U cap = L cap U hyper .
Application of the operator L cap to vector U hyper is done in the finite volume framework:

U cap i,j = U hyper i,j -∆t    1 ∆x F x c;i+ 1 2 ,j U hyper i+ 1 2 ,j -F x c;i-1 2 ,j U hyper i-1 2 ,j + 1 ∆y F y c;i,j+ 1 2 U hyper i,j+ 1 2 -F y c;i,j-1 2 U hyper i,j-1 2    .
The only equations affected by the capillary effects are the momentum and the total energy equations. These equations are developed in two dimensions:

(ρu) cap i,j = (ρu) hyper i,j -σ∆t        - w hyper i+ 1 2 ,j -w hyper i-1 2 ,j ∆x + w 2 1 w hyper i+ 1 2 ,j - w 2 1 w hyper i-1 2 ,j ∆x + ( w 1 w 2 w ) hyper i,j+ 1 2 -( w 1 w 2 w ) hyper i,j-1 2 ∆y        , (ρv) cap i,j = (ρv) hyper i,j -σ∆t        - w hyper i,j+ 1 2 -w hyper i,j-1 2 ∆y + ( w 1 w 2 w ) hyper i+ 1 2 ,j -( w 1 w 2 w ) hyper i-1 2 ,j ∆x + w 2 2 w hyper i,j+ 1 2 - w 2 2 w hyper i,j-1 2 ∆y        , (ρE + ε σ ) cap i,j = (ρE + ε σ ) hyper i,j -σ∆t          w 2 1 u w hyper i+ 1 2 ,j - w 2 1 u w hyper i-1 2 ,j ∆x + ( w 1 w 2 v w ) hyper i+ 1 2 ,j -( w 1 w 2 v w ) hyper i-1 2 ,j ∆x + ( w 1 w 2 u w ) hyper i,j+ 1 2 -( w 1 w 2 u w ) hyper i,j-1 2 ∆y + w 2 2 v w hyper i,j+ 1 2 - w 2 2 v w hyper i,j-1 2 ∆y          .
The different capillary terms at the cell boundary are obtained by the means of an arithmetic average of quantities of neighbouring cells. The vector w requires derivatives of the color function which are computed by using second-order finite difference approximations.

Pressure relaxation operator

To go through the operators chain, the solution at time n + 1 is obtained by a pressure relaxation algorithm and corrects the components of U cap :

U n+1 = L relax (U cap ) .
The details about the pressure relaxation algorithm as well as the correction procedure used to guarantee total energy conservation can be found in Saurel et al. [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF].

Numerical results and validations

In this section, 2D test cases are proposed to prove that the model and the numerical method are able to treat accurately both capillary effects and shock wave propagation.

In each presented cases, the equation of state (EOS) for the air obeys to the ideal gas law:

P air = (γ air -1) ρ air e air , with γ air = 1.4. The liquid obeys the stiffened gas EOS:

P liquid = (γ liquid -1) ρ liquid e liquid -γ liquid P ∞,liquid ,
where the stiffened gas EOS parameters are :

• γ liquid = 2.1 and P ∞,liquid = 1.10 6 P a for the capillary verification test cases of Section 5.1 and 5.2.

• γ liquid = 4.4 and P ∞,liquid = 6.10 8 P a for water in the compressible validation test case of Section 5.3.

Static capillary effects validation

2D tests are proposed to study the ability of the method to treat accurately capillary effects. The proposed method described in the present paper is compared with previous existing methods that use source terms integration. A particular attention is paid to the convergence, the numerical stability, the conservative property, the accuracy with respect to the Laplace pressure law and the dynamics of capillary flows. It is shown that the present method is able to simulate capillary effects in presence of pressure waves. In the following, "conservative method" will always mean the numerical method of Section 4 to distinguish it from the source terms integration method that do not conserve neither the momentum nor the energy.

The first test cases are static tests and they consist in finding the equilibrium state of a cylindrical column of liquid (2D liquid droplet) placed in air. The initial density in air is equal to 1 kg.m -3 and 1000 kg.m -3 in the liquid. The radius of the liquid column R is set to 0.15m and the surface tension coefficient is equal to 800N.m -1 . These unrealistic values are chosen to magnify the model properties. The pressure is initially uniform in the whole domain (75cm x 75cm) and is set to 1 bar (the initial conditions are presented in Figure 2).

Accuracy regarding Laplace law and mesh convergence

The following test cases treating mesh convergence and accuracy regarding the Laplace law are performed using outgoing pressure wave boundary conditions. It consists in imposing the Neumann boundary conditions for the pressure and imposing in-outgoing conditions expressed in terms of the Riemann invariants. Because of the capillary effects, the pressure of the liquid column converges to satisfy the Laplace law. The expression of the theoretical pressure jump of Laplace in 2 dimensions is:

[P ] = σ R ,
where [P ] expresses the pressure jump between inside and outside the droplet, here P liquid -P air . The pressure convergence is analyzed using the relative residual:

= M ax |P n i,j -P n-1 i,j | P n i,j .
The convergence is considered to be reached when the criterion 1.10 -4 is verified. An example of evolution of this relative residual with a mesh of 120x120 cells is presented in Figure 3 for both conservative method and source terms integration method. It is clear that the pressure convergence is obtained faster using the conservative method than using the source terms integration one.

The pressure profiles through the liquid column is presented in Figure 4 for the mesh with 120x120 computational cells. Pressure profiles are presented for both methods after the pressure convergence time mentioned previously. The new conservative method has a better accuracy than the source terms integration one.

The mesh convergence study is performed on 4 different meshes containing 60x60 up to 120x120 computational cells. For each chosen mesh, the relative residual is used to match pressure convergence. Then the relative pressure error ξ is computed by:

ξ = |[P num ] -[P th ]| [P th ] ,
where [P num ] is the pressure jump between the averaged pressure inside the column (P num = α liquid P liquid / α liquid when α liquid 0.99) and the pressure of the air P air . The theoretical pressure jump [P th ] is calculated using the following relation:

[P th ] = σ π S ,
where S represents the 2D volume of the droplet. It is determined doing the sum of the volume fraction of liquid over the whole domain multiplied by the cell volume: S = α liquid S cell . The relative pressure error ξ is averaged in time and plotted versus grid resolution in Figure 5 for the conservative method proposed in the paper. Grid convergence is validated and the variation of ξ over time, when convergence is reached, decreases with the cells refinement.

Stability

Regarding the stability of the method, contours of the liquid column are presented in Figure 6. Wall boundary conditions are used in this example on a mesh containing 120x120 together with cartesian mesh is on the left hand side; Contours after convergence is reached: source terms integration method is in the middle (red contour) and conservative method is on the right hand side (blue contour) for physical times about 3.68s and 2.32s, respectively. With source terms integration method, the position of the column as already changed when convergence is reached.

cells. The left hand side of the figure shows the initial contour together with the mesh. Effect of the mesh on the initial circle is noticeable. The two other pictures show contour positions immediately after the convergence is reached, corresponding to a physical time of about 3.68s for the source term integration method and 2.32s for the conservative method. In this test case, if the source term integration method is used, one can observe a slight move of the gravity center of the column before a converged solution is reached. In fact the gravity center starts to move at a physical time about 2.57s. At the converged time (3.68s), the position of the gravity center is still incorrect (picture in the middle of Figure 6). This drift may be considered as a visible result of the non-conservation of the momentum. This is crippling, for example, when studying interaction between a droplet and a shock wave (see the last subsection): the droplet must be at rest before the interaction with the shock wave. With the conservative method of the present paper, it becomes possible to consider the treatment of this kind of problem.

Nevertheless, the gravity center starts to move after a long time in the two methods. In this example, with the conservative method, it is about 4.52s so well after the convergence is reached (2.32s). But these displacements cannot be compared with the drift observed for the source terms integration method. In other words, one can consider that the new method maintains the position of the gravity center at the converged solution. For longer time, with the accumulation of rounding errors, some spurious motions can be seen as traces of numerical instabilities. They do not have the same origin as the drift observed with the source terms integration method as soon as the calculation begins.

Energy conservation

The conservative property of the new model is verified by plotting the evolution of the relative error for the total mixture energy over the whole domain for the liquid column using wall boundary conditions with a mesh containing 120x120 cells (Figure 7). The energy is obviously conserved with the conservative method but not with the source terms integration 7: Evolution of the relative error for the total mixture energy over the whole domain for the liquid column with wall boundary conditions. Expression for the evolution of the relative error is: 100 (E (t) -E (t 0 )) /E (t 0 ). Source terms integration (red) and conservative (blue) methods.

one. Moreover, for the source terms integration method, the error does not saturate and then still increases when the convergence is reached. It is also noticeable that the difference in energies appears even initially and comes from the capillary energy ε σ which does not appear in the source term form expression of the total mixture energy equation. Note that its value is here: ε σ = 2.10 7 J.m -3 .

Dynamic capillary effects validation

The next two test cases for the dynamic verification of surface tension are slightly different from the previous ones.

Recovering a circle shape

In this test case we are looking for the shape recovering and the volume evolution of an initially square liquid column, with L square = 0.2m, placed in an atmospheric air with the same pressure P air . Because of the capillary terms the solution should evolve to a circle shape of the liquid interface and the pressure should verify the Laplace jump: P f inal = P air + σ/R. The successive dynamic stages leading to the converged circle shape are presented in Figure 8.

Moreover, the mass conservation implies a variation of the liquid column volume. This variation is following an isentropic transformation leading to the final volume:

V f inal = V initial P initial + P ∞ P f inal + P ∞ 1 γ
. Initially the area of the square column is V initial = 0.04m 2 and it converges to the theoretical value of V f inal = 0.039878m 2 with a precision equal to 0.011% at time t = 3s with a mesh made of 200x200 cells. In Figure 9, the relative variation of the liquid 2D volume ((V -V f inal ) /V f inal ) is plotted versus time for 3 different meshes. Grid convergence toward the theoretical value is observed.

Uniform velocity flow

In this section, we focus on the ability of the method to maintain a cylindrical column with a correct pressure jump in a uniform velocity flow. A sketch of the initial condition is presented in Figure 10. It consists in a liquid column with an initial radius of R = 0.15m moving in a uniform velocity field (u = 50m.s -1 ). At the initialization, the liquid column is in equilibrium with the ambient air. Simulations are performed on 4 different meshes containing 180x60 up to 360x120 computational cells for a physical domain of 2.25m x 0.75m.

In Figure 11, the relative pressure error ξ defined previously is averaged in time and plotted versus grid resolution. Mesh convergence is observed and results show that the Laplace pressure jump is well verified.

During the simulations, the position of the liquid column is in perfect agreement with the theoretical one (Figure 12). It is then important to note that surface tension effects have influence only on the shape of the column, since initially there is a grid impact (see Figure 6 for an example), but not on the position. This validates the method for this simple advection test case.

Shock wave interaction with a water column

We focus here on the early stages of droplet breakup in the high speed flow behind a shock wave.

The numerical simulations are performed in two dimensions to simulate the interaction of an air flow and a water column. Comparisons with Igra and Takayma experiments [START_REF] Igra | Numerical simulation of shock wave interaction with a water column[END_REF][START_REF] Igra | A study of shock wave loading on a cylindrical water column[END_REF] are carried out.

A cylindrical water column with initial diameter of D = 6.4mm is exposed to a shock wave at Mach number 1.3 in atmospheric air (see Figure 13 for initialization sketch). The initial densities are ρ air = 1.2kgm -3 and ρ water = 1000kg.m -3 . The corresponding initial Weber number (W e = ρ air u 2 D/σ) in these conditions is 3690. 2D computations are performed on a cartesian mesh containing 3200x1200 cells representing a physical domain of 220mm x 82.5mm. Shocked air is entering at the left. Outgoing pressure wave boundary condition is used at the right and the Newmann boundary conditions is used at the top and the bottom.

The results are presented at three different stages. Each stage is representative of a typical physical effect:

• The first stage is presented in Figure 14 where the flow dynamics is governed by the wave pattern. This is a totally compressible stage corresponding to the shock-column interaction for a large Weber number. Comparisons with experiments are presented in Figure 15.

• After the shock wave propagation, the column shape is changing due to inertial effects. Liquid-gas interface is stretched and the first liquid filaments are appearing. The interface location is compared with experiments in Figure 16.

• At the third stage, the breakup begins because of high inertia of the water column. This effect has already been observed in bubble breakup (Yang et al. [START_REF] Yang | A model for characterization of a vortex pair formed by shock passage over a light-gas inhomogeneity[END_REF], Quirk and Karni [START_REF] Quirk | On the dynamics of a shock-bubble interaction[END_REF], Layes and Le Metayer [START_REF] Layes | Quantitative numerical and experimental studies of the shock accelerated heterogeneous bubbles motion[END_REF]). The Weber number is lower (due to diameter change and also because the velocity gap is reducing) and the capillary effects tend to tear the filaments out the main column (last pictures of Figure 17).

In Figure 18, results obtained for a flow without capillary effects (left picture) and a flow with capillary effects (right picture) are compared at time t = 1.2ms. One can see that the capillary effects begin to have influence on the flow topology, particularly on the filaments shape where the surface tension effects tend to break the filaments to form little droplets.

Conclusion

A new multiphase model treating interface problems and capillary effects has been derived. Hyperbolicity study has been completely done and a splitting numerical method that guarantees conservation of the mass, the momemtum and the energy has been built. Comparison with "source terms integration method" and with experiments have shown the advantages of the new model and numerical method. The simulation of an aerodynamic breakup of a water column induced by a shock wave showed that the method is able to treat accurately both pressure waves interaction and capillary effects in the same formulation. Future works will include viscous effects, heat conduction effects and massive simulations including a large number of droplets in 3D. Here ρ = α 1 ρ 1 + α 2 ρ 2 is the mixture density, e = Y 1 e 1 + Y 2 e 2 is the mixture specific energy, Y i = α i ρ i /ρ are the mass fractions, σ is the surface tension coefficient, c is the color function.

Consider the Hamilton action:

a = t 2 t 1 D LdtdD, (18)
where t 1 and t 2 are fixed time instants, and D is a material fluid volume.

The governing equations are extremal curves of the Hamilton action under the following constraints:

∂ρ ∂t + ∇ • (ρu) = 0, ∂Y 1 ∂t + u • ∇ (Y 1 ) = 0, ∂s 1 ∂t + u • ∇ (s 1 ) = 0, ∂s 2 ∂t + u • ∇ (s 2 ) = 0, ∂c ∂t + u • ∇ (c) = 0. (19)
The Eulerian variations of the unknown variables in terms of virtual displacement δx are given by (see Gavrilyuk [START_REF] Gavrilyuk | Multiphase flow modeling via Hamilton's principle[END_REF] for details):

δρ = -∇ • (ρδx) , δY 1 = -∇Y 1 • δx, δs 1 = -∇s 1 • δx, δs 2 = -∇s 2 • δx, δc = -∇c • δx, δu = ∂δx ∂t + ∂δx ∂x • u - ∂u ∂x • δx.
Taking the variation of the Hamilton action (18) under the boundary conditions:

δx| ∂D = 0, δx| t=t 1 = 0, δx| t=t 2 = 0,
one can obtain the momentum equation:

ρ ∂u ∂t + ∂u ∂x • u + ∇P = -σ (∇ • n) ∇c,
with n = ∇c/ ∇c . The equation is exactly the same as in Brackbill et al. [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF] for the one component case. It can be written in conservative form:

∂ρu ∂t + ∇ • ρu ⊗ u + P I + Ω = 0,
where the capillary tensor Ω is given by:

Ω = -σ ∇c I -n ⊗ n .
The variation of the Hamilton action with respect to the volume fraction gives as the pressure equilibrium condition:

P 1 = P 2 = P = α 1 P 1 + α 2 P 2 .
resulting to the non-conservative equation (2) for α 1 (see Gavrilyuk [START_REF] Gavrilyuk | Multiphase flow modeling via Hamilton's principle[END_REF] for details).

Using the constraints on the entropies and mass fractions of each phase, the mixture entropy equation is given by: ∂s ∂t + u • ∇ (s) = 0.

The mixture entropy, mass and momentum equations admit the following mixture total energy equation:

∂ρE ∂t + ∇ • ((ρE + P ) u) = -∇ • Ω • u.
To transform the previous equation to a fully conservative form, the following relation is used:

∇ • Ω • u = ∇ • Ω • u -Ω : ∂u ∂x ,
where the right term still needs to be transformed. It can be developed as:

σ ∇c I - ∇c ⊗ ∇c ∇c : ∂u ∂x = σ ∇c ∇ • u -∇c • ∂u ∂x • ∇c ∇c . ( 20
)
Using the material derivative equation of ∇c :

d ∇c dt = d∇c dt • ∇c ∇c , (21)
with the material derivative operator d(•)/dt = ∂(•)/∂t + u • ∇(•) and the Schwarz theorem, we have:

d∇c dt = ∇ dc dt -∇c • ∂u ∂x . (22)
Hence with equations ( 19) and [START_REF] Igra | A study of shock wave loading on a cylindrical water column[END_REF], equation ( 21) is rewritten as:

d ∇c dt = -∇c • ∂u ∂x • ∇c ∇c .
Thus relation ( 20) is now developed as:

σ ∇c I - ∇c ⊗ ∇c ∇c : ∂u ∂x = σ ∇c ∇ • u + d ∇c dt = ∂σ ∇c ∂t + ∇ • (σ ∇c u) .
Finally, the energy conservation law is:

∂ρE + ε σ ∂t + ∇ • (ρE + ε σ + P ) u + Ω • u = 0,
with the capillary energy term ε σ = σ ∇c .

Appendix II: Eigenvectors of model [START_REF] Brackbill | A continuum method for modelling surface tension[END_REF] Remind that model (3) can be written under a vector form (6) with the vector W and the matrix A defined by: 0 0 0 0 ρa 2 0 0 u 0 0 0 0 0 0 0 w 1 w 2 w 3 0 u 0 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 0 u 0 0 0 0 0 0 0 0 0 0 0 u

W = [α 1 , u, v,
                    
The corresponding eigenvalues are real, explicit and given by: The system is hyperbolic if the multiple eigenvalues λ = u have exactly linearly independent eigenvectors R. 5 of them are straightforward:

R 1 =                     1 0 0 0 0 0 0 0 0 0 0 0                     R 2 =                     0 0 0 0 0 0 0 0 1 0 0 0                     R =                     0 0 0 0 0 0 0 0 0 1 0 0                     R 4 =                     0 0 0 0 0 0 0 0 0 0 1 0                     R 5 =                                        
And the 3 last eigenvectors for λ = u have to be found using the following reduced matrix: The corresponding reduced eigenvectors r are the solutions of a -uI r = 0, where r is written as: r = [r 1 , r 2 , r 3 , r 4 , r 5 , r 6 , r 7 ] T .

3 distinct systems of equations results from this equality:

ρa 2 r 1 = 0,

  

A solution of the first two equations is: r 1 = 0, r 2 = -ηw 3 and r = ηw 2 , where η is any real constant. If the matrix C is invertible, the components r 5 , r 6 and r 7 are uniquely defined for a given r 4 . First, the derivatives of the capillary tensor terms are:

∂Ω 11 ∂w 1 = σ ρ
3 ).

The determinant is non-zero as long as n 1 is non-vanishing, which is numerically always the case. Thus: 

r =           0 -ηw 3 ηw 2 r 4 -r 4 C -1   1 ρ 0 0            
Thus, only a two-parameter family of eigenvectors exists. Linearly independent eigenvectors can be given as:

R 6 =                     0 0 -w 3 w 2 0 0 0 0 0 0 0 0                     R 7 =                     0 0 0 0 1 -C -1   1 ρ 0 0   0 0 0 0                    
and thus 1 eigenvector is missing to have a hyperbolic model. Model (3) is then weakly hyperbolic.

Figure 1 :

 1 Figure 1: Model (3) admits two types of acoustic waves.

Figure 2 :

 2 Figure 2: Sketch of the initial conditions for the simulations of a cylindrical liquid column (2D liquid droplet) placed in air.

Figure 3 :

 3 Figure3: Relative residual function of the physical time for the 2D liquid column test case with a mesh of 120x120 cells. Source terms integration (red) and conservative (blue) methods. Convergence is reached after 2 seconds with the conservative method compared to more than 6 seconds with the source terms integration method.

Figure 4 :

 4 Figure 4: Pressure jump in a centered section of a steady 2D liquid column in air. σ = 800N.m -1 , R = 0.1496m and so [P ] = 5347P a. Theoretical (black lines) and converged simulations results with source terms integration (red crosses) and conservative (blue dots) methods are shown at times 6.68s and 2.07s, respectively.

Figure 5 :

 5 Figure 5: Relative pressure error ξ averaged in time when the convergence is reached versus the number of cells in one direction N for the 2D liquid column test case. Grid convergence toward the theoretical pressure jump of Laplace is shown for the conservative method proposed in the paper. The bars show the domain of variation of ξ in time after the convergence criteria is reached, associated to each mesh. The interval of variation always decreases with the cells refinement.

Figure 6 :

 6 Figure 6: Magnified view of column contours through color function (c = 0.5). Initial contour (black contour)together with cartesian mesh is on the left hand side; Contours after convergence is reached: source terms integration method is in the middle (red contour) and conservative method is on the right hand side (blue contour) for physical times about 3.68s and 2.32s, respectively. With source terms integration method, the position of the column as already changed when convergence is reached.

Figure

  Figure7: Evolution of the relative error for the total mixture energy over the whole domain for the liquid column with wall boundary conditions. Expression for the evolution of the relative error is: 100 (E (t) -E (t 0 )) /E (t 0 ). Source terms integration (red) and conservative (blue) methods.

Figure 8 :

 8 Figure 8: Magnified schlieren images of the mixture density of the 2D liquid column under a square shape in an air environment converging to a circle shape. Initially σ = 800N.m -1 and L square = 0.2m. Simulation at different times with a mesh of 150x150 cells: (a) initialization, (b) 5ms, (c) 15ms, (d) 20ms, (e) 30ms, (f) 35ms, (g) 45ms, (h) 50ms, (i) 105ms and (j) 3000ms.

Figure 9 :

 9 Figure 9: Evolution of the relative variation of the 2D volume of the liquid column ((V -V f inal ) /V f inal ) for 3 different meshes. The final isentropic 2D volume (S = 0.039878m 2 ) is numerically recovered with an error of 0.011% for the mesh with 200x200 cells.

Figure 10 :

 10 Figure 10: Sketch of the initial conditions for the simulation of a cylindrical liquid column (2D liquid droplet) in air, coupled with advection of the whole domain.

Figure 11 :

 11 Figure 11: Relative pressure error ξ averaged in time versus the number of cells in the y-direction N y for the 2D liquid column test case with surface tension and advection of the whole domain. Mesh convergence is observed. The results show that the Laplace pressure jump is well verified using the new method.

Figure 12 :

 12 Figure 12: Evolution of the position of the liquid column along the x axis minus its initial position for the 2D liquid column test case with surface tension and advection of the whole domain. Numerical results using the new method are in perfect agreements with the theoretical one.

Figure 13 :

 13 Figure 13: Sketch of the initialization for interaction of a shock wave with a water column.

Figure 14 :

 14 Figure 14: Magnified view of the early stages of the aerodynamic breakup of the water column in the highspeed airstream behind the shock wave of Mach number 1.3 in atmospheric air propagating from the left to the right. 2D simulation with a mesh of 3200x1200 cells and schlieren images of the mixture density are given at the following times: (a) 0µs, (b) 10µs, (c) 20µs, (d) 24µs, (e) 37µs and (f) 55µs.

Figure 15 :

 15 Figure 15: Magnified view of the early stages of the aerodynamic breakup of the water column in the highspeed airstream behind the shock wave of Mach number 1.3 in atmospheric air propagating from the left to the right. 2D simulation on the left with a mesh of 3200x1200 cells and a schlieren image of the mixture density. Experiment on the right from Igra and Takayama [22] with an interferogram on the upper half and an unreconstructed hologram on the lower half. Results are shown at times 15µs (top) and 55µs (bottom).

Figure 16 :

 16 Figure 16: Magnified view of the aerodynamic breakup of the water column in the high-speed airstream behind the shock wave of Mach number 1.3 in atmospheric air propagating from the left to the right. 2D simulation on the left with a mesh of 3200x1200 cells and a schlieren image of the mixture density. Experiment on the right from Igra and Takayama [22]. Results are shown at times 200µs (top) and 300µs (bottom).

Figure 17 :

 17 Figure 17: Magnified view of the aerodynamic breakup of the water column in the high-speed airstream behind the shock wave of Mach number 1.3 in atmospheric air propagating from the left to the right. 2D simulation with a mesh of 3200x1200 cells and schlieren images of the mixture density (upper half) coupling with volume fraction of water visualization (lower half, in white the air and in blue the water) are given at the following times: (a) 160µs, (b) 320µs, (c) 420µs, (d) 500µs, (e) 660µs, (f) 820µs, (g) 1000µs and (h) 1200µs.

Figure 18 :

 18 Figure 18: Magnified view of the aerodynamic breakup of the water column in the high-speed airstream behind the shock wave of Mach number 1.3 in atmospheric air propagating from the left to the right. 2D simulations respectively without (left) and with (right) capillary effects taken into account, with a mesh of 3200x1200 cells and a schlieren image of the mixture density. Results are shown at time 1.2ms.

= ρ u 2 2 -

 2 Appendix I: Variational principle for compressible mixtures with capillary effectsConsider a continuum characterized by the Lagrangian: L e -σ ∇c .

  w, P, w 1 , w 2 , w 3 , s 1 , s 2 , Y 1 , c]

λ 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 a 2 +a 2

 1234567822 = u, λ 9,10 = u ± a s , λ 11,12 = u ± a c , b + (a 2 + b) 2 -4a 2 b (n 2 2 + b -(a 2 + b) 2 -4a 2 b (n 2 2 + n 2 3 ) 2 ,a is the Wood mixture speed of sound previously defined (5) and:

w 1 r 1 + w 2 r 2 + w 3 r 3

 123 the matrix with the only capillary terms:



  Finally the resulting reduced eigenvectors are:

1.

Acknowledgments

Authors are particularly grateful to Dan Igra for providing experimental data and to Sarah Hank for fruitful discussions.