Stefano Catalano
email: catalano@unict.it

Salvatore Torrisi

Giuseppe Tortorici

Gino Romagnoli

ACTIVE FOLDING ALONG A RIFT-FLANK: THE CATANIA REGION CASE HISTORY (SE SICILY)

Keywords: active tectonics, contractional structures, morphological responses, eastern Sicily, riftflank deformation and seismotectonics

In this paper a new kinematic and dynamic model on the Recent and active contractional deformation of the Catania region, eastern Sicily, is discussed. The study area represents one of the most seismically active region of the Mediterranean, located at the intersection between the front of the Sicily collision belt and the seismogenic Siculo-Calabrian Rift-Zone. The analysed contractional tectonics form an active triangle zone that originated from the tectonic inversion of a Lower-Middle Pleistocene extensional basin, which was located at the northern edge of the African foreland. The triangle zone consists of two antithetic ENE-WSW oriented thrust-ramps that show evidence of motion during the Holocene and bound a folded belt that involves alluvial deposits as young as 40 ka. These contractional structures represent the final product of the positive tectonic inversion of extensional features located, in the Hyblean Plateau in SE Sicily, along the flank of the active riftzone. The Late Quaternary motions along the inverted structures was accommodated to the west by a major N-S oriented left-lateral fault zone, which separates the active contractional domains from the adjacent sectors of the African margin. As a whole, the Late Quaternary contractional tectonics of SE Sicily have been related to a NW-verging crustal stacking, related to a Mantle intrusion beneath the Hyblean Plateau that developed as effect of the rift flank deformation. The crustal lineaments, which compose the new kinematic model, represents potentially active seismogenic sources that might be considered in the frame of the seismotectonic picture of the Catania region.

Introduction

In the Catania region, in eastern Sicily (Fig. 1), several Authors [START_REF] Labaume | Early deformation at a submarine compressive front: The Quaternary Catania foredeep south of Mt. Etna, Sicily, Italy[END_REF][START_REF] Lanzafame | North-South compression in the Mt. Etna region (Sicily) spatial and temporal distribution[END_REF][START_REF] Monaco | Tettonica pleistocenica nell'area a sud dell'Etna (Sicilia orientale)[END_REF][START_REF] Monaco | Sudden change in the Late Quaternary tectonic regime in eastern Sicily: Evidences from geological and geomorphological features[END_REF] have documented the growth of Middle Pleistocene (> 200 ka) contractional structures, along the southern margin of Mt. Etna A c c e p t e d M a n u s c r i p t Volcano, at the front of the Maghrebian Collision Belt (Fig. 1). More recent studies [START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF]2007) described the occurrence of younger ENE-WSW oriented folds that, coherently with the active NNW-SSE oriented Nubia-Eurasia convergence, deform alluvial terraces as young as the 40 ka. This younger fold system developed along the western flank of the roughly N-S oriented Siculo-Calabrian Rift Zone (SCRZ) (Monaco & Tortorici, 2000;Catalano et al., 2008a), which consists of NNW-SSE trending major extensional faults located in the immediate Ionian off-shore area (Fig. 1). The seismological data [START_REF] Neri | Spatial variation of seismogenic stress orientations in Sicily, south Italy[END_REF] and stress-in-situ measurements along the western flank of the rift zone (Bousquet et al., 1988;[START_REF] Ragg | Patterns of tectonic stress in Sicily from borehole breakout observations and finite element modelling[END_REF] demonstrate the occurrence of a still active compressional regime, which is centred in the Mt. Etna region. This suggests that the growth of contractional structures can be still active and related to the seismicity of the Catania region. In this study, the geometry and the morphological evidence of the Late Quaternary contractional structures of the Catania region are analysed in order to define the displacement-rate and the mechanism of their evolution in the frame of the regional dynamics and seismotectonics of this densely populated area. The methodology includs morphometric analyses, applied on the Late Pleistocene alluvial terraces and the Holocene alluvial plain of the Simeto and Dittaino rivers (Catania Plain), and the study of channel morphology along the two main streams. The morphological data are augmented by surface geological information, derived from 1:10.000 scale field mapping, combined with the results of electrical resistivity tomography, performed across the Holocene alluvial plain.

Regional setting

Eastern Sicily is located along the Nubia-Eurasia convergent boundary [START_REF] Dewey | Kinematics of the western Mediterranean. Alpine Tectonics[END_REF][START_REF] Serpelloni | Kinematics of the Western Africa-Eurasia plate boundary from focal mechanisms and GPS data[END_REF], where the Neogene-Quaternary southeastern-verging Maghrebian Thrust Belt originated from the tectonic inversion of the former Africa continental palaeomargin (see inset in Fig. 1). The geometry of the collision belt consists of two imbricated continental blocks (Panormide and Pelagian blocks) and a thin-skinned Neogene-Quaternary accretionary wedge, A c c e p t e d M a n u s c r i p t which is sandwiched between them [START_REF] Roure | Late Cenozoic subduction complex of Sicily[END_REF])(Fig. 2). A main E-W oriented, N-dipping crustal oblique (dextral) thrust (Mt. Kumeta-Alcantara Line of [START_REF] Ghisetti | Different styles of deformation in the Calabrian Arc (Southern Italy); implication for a sismotectonics zoning[END_REF], which crops out in northern Sicily, represents the leading edge of the Panormide Block, also playing the role of back-stop of the Neogene-Quaternary accretionary wedge [START_REF] Roure | Late Cenozoic subduction complex of Sicily[END_REF]. This is emplaced on the Pelagian continental crust that, being buried beneath the entire eastern Sicily, culminates in the south-eastern sectors of the island, where it forms the Hyblean foreland domain [START_REF] Ben Avraham | Principali domini strutturali dalla collisione continentale neogenico-quaternaria nel Mediterraneo centrale[END_REF]. The present geometry of the Pelagian Block is the result of the Neogene-Quaternary flexure beneath the SE-verging orogenic belt. During Quaternary; this process was associated with the reactivaction of a main N-S oriented dextral transform fault zone, the Scicli Line [START_REF] Grasso | The western margin of the Hyblean Plateau: a neotectonic transform system on the SE Sicilian foreland[END_REF])(Fig. 1). To the west of this fault zone, the Hyblean Foreland sequences have also been flexured to form the Early Pleistocene NE-SW oriented Gela Foredeep [START_REF] Grasso | The western margin of the Hyblean Plateau: a neotectonic transform system on the SE Sicilian foreland[END_REF][START_REF] Grasso | Pliocene-Pleistocene tectonics on western margin of the Hyblean plateau and the Vittoria Plain (SE Sicily)[END_REF]. To the east of the Scicli Line, the Hyblean Foreland sequences form a wide plateau, bordered to the north by two Middle Pleistocene collapsed basins, the Scordia-Lentini Graben and the Simeto Graben (Fig. 1) [START_REF] Ghisetti | The structural features of the Hyblean Plateau and the Mount Judica area (South-Eastern Sicily): a microtectonic contribution to the deformational history of the Calabrian Arc[END_REF] that are superimposed on the previous foredeep basin developed in the Catania region. Recent studies [START_REF] Bousquet | Compression and Quaternary tectonic inversion on the Northern edge of the Hyblean Mountains, foreland of the Appenine-Maghrebian chain in Eastern Sicily (Italy): geodynamic implications for Mt[END_REF]Catalano et al., 2008b) have documented a Late Quaternary tectonic inversion along both the Scicli Line, from right-to left-lateral, and along the margins of the associated Scordia-Lentini Graben, from extensional to contractional. As a consequence, during the Late Quaternary, the eastern portions of the Hyblean Plateau have been displaced to the northwest, towards the Catania Plain, at a rate estimated to be about 1.4 mm/a (Catalano et al., 2008b).

The entire collision belt of eastern Sicily is confined, to the east, by the Siculo-Calabrian Rift Zone (SCRZ) (Monaco & Tortorici, 2000); this is the tectonic expression of an incipient divergent margin, which separates the Adria-Ionian microplate from the colliding Tyrrhenian and Pelagian lithospheric domains [START_REF] Serpelloni | Kinematics of the Western Africa-Eurasia plate boundary from focal mechanisms and GPS data[END_REF] (Fig. 1). The rift zone, in the Catania area, is composed of NNW trending oblique (dextral) extensional faults that are distributed from the eastern flank of Mt.Etna to the Ionian off-shore of the Hyblean Plateau, [START_REF] Monaco | Late Quaternary slip rates on the Acireale-Piedimonte normal faults and tectonic origin of Mt. Etna (Sicily)[END_REF][START_REF] Monaco | Tettonica pleistocenica nell'area a sud dell'Etna (Sicilia orientale)[END_REF] A c c e p t e d M a n u s c r i p t Tortorici, 2000;Catalano et al., 2008a).

Late Quaternary deformation and seismotectonics of Sicily

Sicily has been affected by intense Late Quaternary deformation, associated with high-level seismicity, which can be related to distinct seismogenic domains [START_REF] Neri | Spatial variation of seismogenic stress orientations in Sicily, south Italy[END_REF]. Low to medium magnitude events (4<M<6) affect the deformation belt, which is related to the Eurasia-Nubian convergent boundary (see inset a in Fig. 1). The seismicity is concentrated along the Tyrrhenian off-shore and the northwestern sectors of the island, where available focal mechanisms constrain a stress field characterised by a N130° oriented maximum compression. Similarly, focal mechanisms indicating motions along NE-SW oriented thrust fault have been reconstructed for the 1968 Belice seismic sequence (4.6≤M≤5.4; [START_REF] Anderson | Active tectonics in the Adriatic region[END_REF][START_REF] Monaco | Active thrust tectonics in western Sicily (Southern Italy): the 1968 Belice earthquake sequence[END_REF] that affected the external areas of the orogen, in the southwestern Sicily (see inset a in Fig. 1).

In eastern Sicily, the present-day stress field has a different orientation because the σ hmax follows the trend of the SCRZ, being NE-SW, in the Eolian Islands, and NNW-SSE, in the Mt.

Etna region [START_REF] Neri | Spatial variation of seismogenic stress orientations in Sicily, south Italy[END_REF]. A N145° oriented σ 1 has been obtained by the inversion of seismological data in the Hyblean region [START_REF] Musumeci | Stress Directions and Shear-Wave Anisotropy: Observations from Local Earthquakes in Southeastern Sicily, Italy[END_REF]. This evidence is in good agreement with the focal mechanism of the 1990 Augusta earthquake [START_REF] Amato | Main shock and aftershocks of the December 13, 1990, Eastern Sicily earthquake[END_REF], that indicates a NW-SE compression and a NE-SW tension (see inset a in Fig. 1).

Seismological data at a regional scale are consistent with the geodetic deformation that can be obtained by GPS-determined velocities with respect to stable Eurasia. Motions are constantly parallel to the maximum compression of the stress tensor. At sites LAMP and TRAP (see inset in Geodetic deformation also indicates that about 2.6 mm/a of shortening, representing the 63% of the entire NW-SE oriented Africa-Eurasia convergence, is accommodated by folding and thrusting within the Sicily collision belt (NOTO site vs. USTICA site; see inset a of Fig. 1), while the residual convergence (1.5 mm/a) measured at Ustica site must be accommodated further to the north.

In-situ stress measurements also fit well the seismological and geodetic data, as they indicate mean σ hmax orientation of about N131°, in western Sicily, and N148° in the Hyblean foreland areas [START_REF] Ragg | Patterns of tectonic stress in Sicily from borehole breakout observations and finite element modelling[END_REF].

Nevertheless, the main destructive (M> 6.5) historical earthquakes of Sicily (A.D. 1169, 1542, 1693) have been localised along the Ionian coast of the island [START_REF] Postpischl | Catalogo dei terremoti italiani dall[END_REF][START_REF] Boschi | Tectonic stresses and volcanism: in-situ stress measurements and neotectonic investigations in the Etna area (Italy)[END_REF][START_REF] Monaco | Tettonica pleistocenica nell'area a sud dell'Etna (Sicilia orientale)[END_REF][START_REF] Bianca | Quaternary normal faulting in southeastern Sicily (Italy): A seismic source for the 1693 large earthquake[END_REF]Monaco & Tortorici, 2000). They define a major seismogenic zone which is related to the occurrence of the active normal fault segments of the SCRZ (Fig. 1). Available crustal focal mechanisms in the Straits of Messina area [START_REF] Anderson | Active tectonics in the Adriatic region[END_REF]; see inset a in Fig. 1), coupled with kinematic analyses on fault planes along the eastern flank of Mt.Etna [START_REF] Monaco | Late Quaternary slip rates on the Acireale-Piedimonte normal faults and tectonic origin of Mt. Etna (Sicily)[END_REF]Monaco & Tortorici, 2000), indicate that motion along the normal fault belt occurred in response of a N100° oriented maximum extension. The analysis of seismic lines across the offshore active fault segments [START_REF] Bianca | Quaternary normal faulting in southeastern Sicily (Italy): A seismic source for the 1693 large earthquake[END_REF] and morpho-structural investigations carried out along the on-shore fault planes [START_REF] Monaco | Late Quaternary slip rates on the Acireale-Piedimonte normal faults and tectonic origin of Mt. Etna (Sicily)[END_REF] suggested Late Quaternary average extensional-rates ranging from 1.7 to 3.7 mm/a. These data are in good agreement with GPS measurements carried out in southern Italy (Hollestein et al., 2003). In more detail, a mean active extension of about 3.5 mm/a has been calculated in the Messina Straits Region (D'Agostino & [START_REF] Selvaggi | Crustal motion along the Eurasia-Nubia plate boundary in the Calabrian arc and Sicily and active extension in the Messina Straits from GPS measurements[END_REF].

Late Quaternary evolution of the Catania region

The town of Catania is located on the northern margin of the NE-SW oriented Middle A c c e p t e d M a n u s c r i p t Pleistocene Simeto Graben (Fig. 1) that originated at the front of the Maghrebian Thrust Belt and is filled by more than 1000 m of pelagic sediments (1.2 to 0.6 Ma-old;[START_REF] Di Stefano | Long-term uplift rate of the volcano basement (southern Italy) based on biochronological data from Pleistocene sediments[END_REF].

The seismic lines across the basin and bore-hole data [START_REF] Torelli | 5) Early-Middle Pleistocene deposits of the Gela Foredeep; (6) Undifferentiated Meso-Cenozoic Maghrebian allochthonous units; (7) Pliocene-Early Pleistocene volcanics of the Hyblean Plateau; (8) Meso-Cenozoic foreland sequences of the western Hyblean Plateau; (9) Meso-Cenozoic foreland sequences of the eastern Hyblean Plateau; (10) buried front of the Maghrebian allochthonous units; (11) Pliocene-Early Pleistocene fault; (12) Middle Pleistocene extensional fault (Simeto and Scordia-Lentini Graben); (13) Late Quaternary and active tectonics and volcano-tectonics: strike-slip fault (a); normal fault of the "Siculo-Calabrian Rift Zone" (b), thrust (c), anticline (d) dry and eruptive fissure (e); (14[END_REF]Fig. 3) emphasize that the Simeto Graben sequence seals the tectonic contact of the frontal nappe of the allochthonous edifice (Gela Nappe) upon the Early Pleistocene horizons of the Catania Foredeep (P-Q in Fig. 3) [START_REF] Lentini | Tentativo di elaborazione di uno schema strutturale della Sicilia Orientale[END_REF][START_REF] Lentini | Tentativo di elaborazione di uno schema strutturale della Sicilia Orientale[END_REF][START_REF] Bianchi | Sicilia orientale: Profilo geologico Nebrodi-Iblei[END_REF]Lentini et al., 1996 and references therein), which are located at the top of the flexured Hyblean carbonate sequences [START_REF] Roure | Late Cenozoic subduction complex of Sicily[END_REF][START_REF] Lentini | Main structural domains of the central mediterranean region and their Neogene tectonic evolution[END_REF](Fig. 3). The southern margin of the Simeto Graben is represented by the Primosole Horst (Fig. 3), where Middle Pleistocene extensional features are still preserved [START_REF] Torelli | 5) Early-Middle Pleistocene deposits of the Gela Foredeep; (6) Undifferentiated Meso-Cenozoic Maghrebian allochthonous units; (7) Pliocene-Early Pleistocene volcanics of the Hyblean Plateau; (8) Meso-Cenozoic foreland sequences of the western Hyblean Plateau; (9) Meso-Cenozoic foreland sequences of the eastern Hyblean Plateau; (10) buried front of the Maghrebian allochthonous units; (11) Pliocene-Early Pleistocene fault; (12) Middle Pleistocene extensional fault (Simeto and Scordia-Lentini Graben); (13) Late Quaternary and active tectonics and volcano-tectonics: strike-slip fault (a); normal fault of the "Siculo-Calabrian Rift Zone" (b), thrust (c), anticline (d) dry and eruptive fissure (e); (14[END_REF].

Conversely, the northern border of the basin has been completely obliterated by younger (200 ka) ramp-thrusts that bound the Dittaino-Simeto and the Terreforti highs [START_REF] Monaco | Tettonica pleistocenica nell'area a sud dell'Etna (Sicilia orientale)[END_REF][START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF]. Seismic lines across the deposits of the Simeto Graben image an incipient NW-ward gently dipping detachment that ramps from the base of the allochthonous units through the Middle Pleistocene horizons. This detachment has been recognised at a depth of about 2000 m, beneath the Catania Plain and it approaches the surface near the coast (Fig. 3; [START_REF] Torelli | 5) Early-Middle Pleistocene deposits of the Gela Foredeep; (6) Undifferentiated Meso-Cenozoic Maghrebian allochthonous units; (7) Pliocene-Early Pleistocene volcanics of the Hyblean Plateau; (8) Meso-Cenozoic foreland sequences of the western Hyblean Plateau; (9) Meso-Cenozoic foreland sequences of the eastern Hyblean Plateau; (10) buried front of the Maghrebian allochthonous units; (11) Pliocene-Early Pleistocene fault; (12) Middle Pleistocene extensional fault (Simeto and Scordia-Lentini Graben); (13) Late Quaternary and active tectonics and volcano-tectonics: strike-slip fault (a); normal fault of the "Siculo-Calabrian Rift Zone" (b), thrust (c), anticline (d) dry and eruptive fissure (e); (14[END_REF]. The deposits of the Simeto Graben are now covered by the alluvial plain of the Simeto and Dittaino Rivers (Catania Plain).

The Simeto Graben sequence is extensively exposed in the hangingwall of the thrust bordering the basin to the north. Along the Terreforti High (Fig. 3), Late Quaternary clastic deposits engrave several orders of terraced abrasion platforms, which have been modelled at the top of the Middle Pleistocene clay sequence. The terraced deposits represent regressive marine sequences (S 1 -S 7 in [START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF], upward grading into alluvial conglomerates (T 1 -T 7 in [START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF] see also [START_REF] Chester | The interaction of volcanic activity in Quaternary times upon the evolution of the Alcantara and Simeto rivers[END_REF][START_REF] Monaco | Tettonica pleistocenica nell'area a sud dell'Etna (Sicilia orientale)[END_REF]Monaco et al., 2000)(Tab. 1), that were accumulated during successive periods of highstand, under the control of active tectonic uplift and deformation (Fig. 3). The fluvial terraces along the left side of the Simeto River [START_REF] Chester | The interaction of volcanic activity in Quaternary times upon the evolution of the Alcantara and Simeto rivers[END_REF]) can be, thus, interpreted as cyclic features, strictly related to the eustatically A c c e p t e d M a n u s c r i p t controlled marine sequences distributed along the coastal slope of the town of Catania. This implies that the age of the marine deposits can be extended also to the correlative alluvial deposits [START_REF] Caputo | Late Quaternary activity of the Scorciabuoi Fault (southern Italy) as inferred from morphotectonic investigations and numerical modeling[END_REF].

The older terraced sequence (S 1 and T 1 in [START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF])(Tab. 1) rests on an abrasion platform undercutting 320-250 ka-old subaerial subalkaline lavas [START_REF] Gillot | The evolution of Mount Etna in the light of potassiumargon dating[END_REF] andis intruded by 200 ka-old sub-volcanic rocks [START_REF] Del Negro | The age of the neck of Motta Sant'Anastasia (Sicily) as revealed by a magnetic survey[END_REF], while the youngest sequence (S 7 and T 7 in [START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF])(Tab. 1) is largely covered by subaerial lava flows from the modern central volcano (< 35 ka-old; [START_REF] Gillot | The evolution of Mount Etna in the light of potassiumargon dating[END_REF]. Taking into account the relation between eustacy and terracing, the age of the terraced sequences is thus constrained within the OISs (Oxygen Isotope Timescale Stages) to have occurred since 250 ka B.P., as also demonstrated by occurrence of volcanic pebbles from the Etnean sub-alkaline lavas (320-250 ka old; [START_REF] Gillot | The evolution of Mount Etna in the light of potassiumargon dating[END_REF] in the coarse grained levels [START_REF] Monaco | Sudden change in the Late Quaternary tectonic regime in eastern Sicily: Evidences from geological and geomorphological features[END_REF].

In more detail, the lowermost five terraced sedimentary wedges and the correlative alluvial terraces (T 3 -T 7 in [START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF]) (Tab. 1) undercut the 180 ka-old alkaline volcanics and can be thus related to the five main highstands between the OIS 5.5 (125 ka) and 3.1 (40 ka)(Tab. 1). This indirect dating, based on geomorphic evidence, is also supported by the occurrence of volcanic pebbles from the ancient alkaline lavas in the coarse-grained levels [START_REF] Monaco | Tettonica pleistocenica nell'area a sud dell'Etna (Sicilia orientale)[END_REF]. This implies that the lowermost fluvial terrace (T 7 ; Tab. 1), flanking the courses of the Simeto and Dittaino Rivers, can be assigned to the OIS 3.1 (40 ka).

The Middle Pleistocene (240-200 ka) terraced deposits are severely deformed by the contractional structures developed by the positive tectonic inversion of the northern margin of the Simeto Graben as a result of the propagation of a dextral shear zone from the collision belt. These structures are mostly represented by E-W oriented ramp-thrusts associated with major NW-SE oriented dextral strike-slip fault zones (Fig. 3 and Tab. 1; [START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF]. In contrast, the Upper Pleistocene deposits, ranging in age from 125 to 40 ka, have been severely uplifted, as a consequence of the deformation along the extensional faults of the SCRZ, located in the Ionian off-A c c e p t e d M a n u s c r i p t shore region (Tab. 1; [START_REF] Monaco | Sudden change in the Late Quaternary tectonic regime in eastern Sicily: Evidences from geological and geomorphological features[END_REF][START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF]. During the last 40 ka, the syn-rift deformation has been accompanied by the propagation through the Catania Plain of a set of ENE-WSW oriented gentle folds (Contrada Agnelleria, Contrada Rotondella and Sigonella anticlines in Fig. 3; [START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF], that has deformed the most recent alluvial terrace (T 7 ; 40 ka)(Tab. 1).

5.

The most impressive contractional structure of the Catania region is represented by a N70 oriented, northwest-dipping ramp thrust. It involves terraced deposits as young as 80 ka [START_REF] Monaco | Sudden change in the Late Quaternary tectonic regime in eastern Sicily: Evidences from geological and geomorphological features[END_REF], as observed in the southern sector of the town of Catania, in a quarry (Fig. 4). In this area, the alluvial conglomerates belonging to terrace T 5 (OIS 5.1;Monaco et al., 2000) (Tab. 1) unconformably cover the 200 ka-old marine sandy sequence (S 2 )(see inset a of Fig. 4), which is severely tilted towards the south, at the eastern termination of the Terreforti High (Fig. 3; [START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF]. The syn-rift thrust-ramp is responsible for the local overthrusting of the marine terraced sequence (S 2 ; 200 ka) onto the alluvial conglomerates of the T 5 (80 ka)(see inset b of Fig. 4). The total vertical displacement along the structure has been estimated at about 10 m, considering the off-set of the 80 ka-old unconformity across the ramp. In the hanging-wall of the structure, the motion along the ramp has produced the propagation of a fan of minor slip planes that isolate a N70° oriented kink band, which has further tilted the S-dipping strata (Fig. 4).

The described thrust-ramp developed at the forelimb of the N70° oriented Sigonella Anticline (Fig. 3), a wide box-fold that extends towards the SW, as far as 13 km from the urban area of 5) shows a slight N-ward tilting of the layering, with the consequent S-ward narrowing of the Holocene alluvial deposits from 25 to 15 m, on a distance of about 300 m (dip=2°). The tomographic model crossing the buried crest of the Sigonella Anticline (profile 2 in Fig. 5) reveals a further narrowing (< 5 m) of the Holocene alluvial deposits. In the resistivity pseudo cross-section, the forelimb of the fold is also marked by the occurrence, at depth, of a NEdipping body that, characterised by higher resistivity values (20 to 35 Ohm•m), cutting at high angle the sub-horizontal layering. The geometry of this resistivity anomaly is consistent with the occurrence, as well as in the Catania area, of heavy fractured fault-rocks, associated with a NEdipping thrust-ramp. As similar to the Catania town section, the tomographic model across the structure suggests a very small vertical displacement that can be estimated at about 10 m, if A c c e p t e d M a n u s c r i p t measured on the offset of the base of the sandy levels, while it is reduced to about 5 m, considering only the offset of the base of the Holocene alluvial deposits.

In the northwestern sector of the Catania Plain and in the adjacent Terreforti High, further two Late Quaternary N70 oriented anticlines (Contrada Rotondella and Contrada Agnelleria anticlines of [START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF]) can be recognised (Fig. 3). These two gentle folds show very low amplitude (20 and 40 m) and long wavelength (3 to 4 km, respectively). The Contrada Agnelleria and Contrada Rotondella folds appear as slight flexures of the top-surface of the youngest terraced alluvial deposits (T 7 ; 40 ka) (Tab. 1), showing interlimb angles that approximate the 180° (Fig. 6).

The geometry of these folds shows almost uniform amplitude and length on surface of different age, where they deform the older terraced alluvial deposits (T 1 -T 6 ; Tab. 1). The morphology of the folds is completely absent in the Holocene floodplain, which shows a very low-gradient (0,3%) longitudinal profile, almost parallel to that of the active channel of the Simeto River.

Out of the folded area, the T 7 terrace (Tab. 1) is regularly dipping, at very low-angle (0.34°), towards the Ionian sea. In contrast, in the folded area, the T 7 terrace shows evident undulations, as it culminates to form ENE-WSW elongated ridges, corresponding with the well-rounded hinge zones of the anticlines, which are separated by wide U shaped topographic depressions that have developed where the terraced surfaces have been subsided at the core of the synclines (Fig. 6).

Approaching the crest of the Contrada Agnelleria and Contrada Rotondella anticlines, the T 7 terrace (Tab. 1) shows a slight divergence (1.34°) from the alluvial plain, coupled with convergence (0.40°) towards the intermediate synclines. The steeper and shorter NNW-dipping limbs suggest a NNW vergence of the folds. To the southeast of the Rotondella anticline, the 40 ka-old fluvial deposits have been covered by the Holocene deposits of the Catania Plain.

The syn-rift contractional deformation has also controlled the evolution of the Holocene fluvial morphology of the Catania Plain, which formed by the merging of the Simeto and Dittaino river floodplains (Fig. 3). The final tract of the Simeto River is characterised by an arcuated meandering course that, varying in direction from the N-S to the NW-SE, crosses at high-angle the

A c c e p t e d M a n u s c r i p t

Late Quaternary anticlines of the Catania region (Fig. 3). At a distance of about 25 km from the mouth, the active channel of the Simeto River as well as the top of the Holocene floodplain show a very low-gradient (about 0.3 %), crossing undisturbed the folds that deform the 40 ka-old alluvial terrace (T 7 ).

In the Catania Plain, a very impressive change in the sinuosity index (S) of the river channel occurs at the intersection with the Sigonella Anticline (from the reach 6 to 7; Fig. 7). This variation, even if partly obliterated by modification due to the human activities, has been estimated taking into account the topographic data from the 1867 1:50.000 topographic map, the 1:100.000 topography reported in the geological map of the end of the XIX century (Gemmellaro and Travaglia, 1885;Gemmellaro et al., 1885) and, finally, the 1960 edition of I.G.M.I. 1:25.000 topographic map. If measured on the XIX century topography, the variation of the sinuosity can be estimated from 1.75 to 3.75. A reduced change of the sinuosity (from 1.75 to 2.5) is imaged in the 1960 topography, as a consequence of the cut-offs of meanders that occurred in the last century.

To the west of the Simeto River, the Dittaino River is characterised by a prevalent NW-SE oriented meandering course, which flows from the external areas of the Maghrebian Orogen to the Catania Plain, crossing roughly perpendicular the NE-SW oriented termination of the Scicli Line.

At the intersection with the structure, which is located between the reach 2 and the reach 3 (Fig. 7), the Dittaino River shows an evident increase of the sinuosity, of about the 100% (from 1.5 to 3), that was attenuated at about 60%, due to the cut-off of meandersin the reach 3 during the last century. From the Scicli Line to the mouth, the sinuosity of the river remains almost uniform and constant in the time. Nevertheless, the Dittaino displays a large scale deflection, which bounds the western terminations of the Rotondella and Sigonella anticlines, as it enters in the Catania Plain (Fig. 7).

A c c e p t e d M a n u s c r i p t

Discussion and conclusions

A Late Quaternary fold and thrust system that developed along the flank of an active rift zone in the Catania region has been described. The Late Quaternary contractional structures as a whole generate a triangle zone, involving the Lower-Middle Pleistocene deposits of a former extensional basin (Simeto Graben) and deforming Upper Pleistocene and Holocene morphological features (Fig. 3). The south-verging active boundary of the triangle zone is represented by the ramp-thrust outcropping in the town of Catania, here referred to as Catania-Sigonella Thrust (CST), responsible for adjustments in the fluvial morphology along the Simeto River. At the intersection between the river and the structure, a dramatic change in the sinuosity is documented by several historical topographic maps (Fig. 8). The synoptic analysis of the 1:50.000 scale topographic map, edited in the 1867, of the topography reported in the 1:100.000 scale geological map of the end of the XIX century (Gemmellaro and Travaglia, 1885;Gemmellaro et al., 1885) and of I.G.M.I. 1:25.000 topographic map, edited in 1960, demonstrates that the difference in the sinuosity, between reaches 6 and 7 (Fig. 8), has in about 100 yrs progressively decreased from 2 to 0.75. The CST accumulated a total vertical displacement of about 10 m, half of which is related to Holocene motions as suggested by the clear control on the fluvial morphology. This implies that the structures developed in the last 20 ka, with an averaged vertical displacement-rate of about 0,5 mm/a. The opposite north-verging tectonic boundary of the active triangle zone can be located along the northern limb of the Agnelleria Fold, on the prolongation of the active NE-SW oriented splay of the Scicli Line, here designated as the Sferro-Paternò Thrust (SPT), displacement on which has caused the morphological responses recognised along the Dittaino River. The intersection between the river and the structure is also in this case marked by a well defined change in sinuosity, only partially attenuated in the last century. The recognised responses of the fluvial morphology, which have been modelled in the literature by [START_REF] Holbrook | Geomorphic and sedimentary of rivers to tectonic deformation : a brief review and critique of a tool for recognizing subtle epeirogenic deformation in modern and ancient setting[END_REF], document very recent longitudinal tilting of the active channels of the Simeto and Dittaino Rivers, at the intersection with Conversely, the undeformed Holocene morphological features, which are incised into the entire fold belt, and the geometry of the Late Pleistocene fluvial terraces suggest that the folding is restricted to the time interval from 40 to 10 ka. Considering the evolution of these contractional structures, the Late Pleistocene anticlines can be interpreted as fault-propagation folds affected by opposite vergences, that deactivated when the thrusts reached the surface. Thereafter, the entire triangle zone was uplifted, confined to the west by the large scale deflection of the Dittaino River.

The cumulative shortening across the entire triangle zone can be estimated taking into account the vertical displacement of the Late Quaternary terraces and the shortening due to the gentle folding (< 3 m). The terracing of the 40 ka-old fluvial deposits is confined at the hangingwall of the SPT. These old alluvial sediments were deposited during OIS 3.1, when the base level was about 40 m lower than Present [START_REF] Bassinot | The astronomical theory of climate and the age of the Brunhes-Matuyama magnetic reversal[END_REF]. Considering an average elevation of the terrace of about 30 m above the modern Simeto River, a cumulative uplift of about 70 m is implied. A large portion of this uplift can be related to the regional signal that has been estimated to be about 0.8 mm/a [START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF]. The residual amount of uplift (38 m) can be interpreted as the tectonic contribution of the SPT. Considering that the reverse faulting developed along a former extensional margin, a dip of 60° for the inverted fault-plane can be tentatively assumed. This would imply a conversion of the vertical motion into a shortening of about 22 m. Taking also into account the motion along the CST and the gentle fold deformation, a shortening of about 35 m has been accommodated by the entire triangle zone in the last 40 ka, at a deformation-rate of about 0.8-0.9 mm/a. A rough estimate of the depth of the sole-thrust of the triangle zone can be obtained comparing the amount of the relief produced at the surface by the Late Quaternary thrust motion with the related cumulative shortening. Assuming plane-strain deformation, a depth of the solethrust of about 7500 m is obtained by analysing the deformation of the 40 ka-old alluvial terrace (T 7 ) along the topographic profiles, drawn perpendicular to the trend of the contractional structures A c c e p t e d M a n u s c r i p t (Fig. 6). This implies, considering a depth of about 2000 m of the sole-thrust of the allochthonous units beneath the Catania Plain (T 1 in Fig. 3), that the relation between the Late Quaternary deformation of the triangle zone and the migration of the front of the orogen is unrealistic.

The geometry, the estimated depth of the basal detachment and the deformation-rate of the triangle zone are, however consistent with the Late Quaternary contractional deformation affecting the Hyblean crustal block, which is bounded to the west by the Scicli Line. In this region, the Early Pleistocene extensional features (e.g. Scordia-Lentini Graben) experienced a diffuse Late Quaternary positive tectonic inversion, associated with the inversion of the sense of motion along the Scicli Line (Catalano et al., 2008b). The authors have constrained the positive tectonic inversion of the Scordia-Lentini Graben, connected to the left-lateral motion along the southernmost segments of the Scicli Line, at about 850 ka. Recent studies [START_REF] Catalano | Rift-flank deformation in SE Sicily: evidence for an active Mantle upwelling[END_REF] have interpreted the inversion tectonics as related to a rift-flank deformation, affecting the eastern sectors of the Hyblean Plateau, located at the footwall of the SCRZ. The authors have related the surface inverted tectonics of the southern margin of the Scordia-Lentini Graben to the emergence of a NNW-verging crustal ramp that originated from a Mantle intrusion, which has been inferred at a depth of about 25 km.

It is worth noting that the described surface inverted structures have accommodated only a small amount (30%) of the sinistral strike-slip motion along the Scicli Line, thus suggesting that part of the contraction should be taken up along blind structures, located more to the north. The present study suggests that a NNW-ward propagation of the inverted contractional structures may have migrated for about 30 km in the last 850 ka (Fig. 8). The overall surface geometry of these younger inverted structures is compatible with a detachment, located at depth of about 7500 m, along the Permo-Triassic levels at the base of the sedimentary succession of the Hyblean Plateau [START_REF] Bianchi | Sicilia orientale: Profilo geologico Nebrodi-Iblei[END_REF], that ramps to form the SPT, since about 40 ka. In this context, the CST, which forms the opposite boundary of the triangle zone would represent a secondary back-thrust, characterised by a very low cumulative displacement. The measured shortening-rate (0.8-0.9 mm/a) along the inverted triangle zone of the Catania Plain is comparable with the residual left-lateral A c c e p t e d M a n u s c r i p t deformation-rate (≈1 mm/a), affecting the northernmost segments of the Scicli Line, to the north of the Scordia-Lentini Graben. According to this new kinematic model, the active contractional deformation of the Catania region can be related to the dynamics affecting the whole of southeastern Sicily, driven by mantle intrusion triggered by the rift-flank deformation along the SCRZ, rather than to the late evolution of the collision belt [START_REF] Catalano | Rift-flank deformation in SE Sicily: evidence for an active Mantle upwelling[END_REF].

This new kinematic and dynamic model has consequences for the regional seismotectonic picture, if the contractional deformation of the Catania region is to be related to motion along the northern segments of the Scicli Line and the southern inverted margins of the Scordia-Lentini Graben. These tectonic alignments represent major crustal structures that must be considered as potential seismogenic sources for locating historical earthquakes. As well, the measured deformation rates are significant and should be accounted for when evaluating the recurrence-time of seismicity in the area. [START_REF] Postpischl | Catalogo dei terremoti italiani dall[END_REF][START_REF] Boschi | Tectonic stresses and volcanism: in-situ stress measurements and neotectonic investigations in the Etna area (Italy)[END_REF][START_REF] Monaco | Tettonica pleistocenica nell'area a sud dell'Etna (Sicilia orientale)[END_REF][START_REF] Bianca | Quaternary normal faulting in southeastern Sicily (Italy): A seismic source for the 1693 large earthquake[END_REF]Monaco & Tortorici, 2000; Focal mechanisms are from [START_REF] Anderson | Active tectonics in the Adriatic region[END_REF][START_REF] Amato | Main shock and aftershocks of the December 13, 1990, Eastern Sicily earthquake[END_REF]Pondrelli et al., 2002;2004); ( 18 In the geological profiles, which has been drawn on the basis of surface data, supported by bore-hole logs and seismic lines [START_REF] Torelli | 5) Early-Middle Pleistocene deposits of the Gela Foredeep; (6) Undifferentiated Meso-Cenozoic Maghrebian allochthonous units; (7) Pliocene-Early Pleistocene volcanics of the Hyblean Plateau; (8) Meso-Cenozoic foreland sequences of the western Hyblean Plateau; (9) Meso-Cenozoic foreland sequences of the eastern Hyblean Plateau; (10) buried front of the Maghrebian allochthonous units; (11) Pliocene-Early Pleistocene fault; (12) Middle Pleistocene extensional fault (Simeto and Scordia-Lentini Graben); (13) Late Quaternary and active tectonics and volcano-tectonics: strike-slip fault (a); normal fault of the "Siculo-Calabrian Rift Zone" (b), thrust (c), anticline (d) dry and eruptive fissure (e); (14[END_REF], the tectonics related to distinct deformation stages are reported: T 1 ) sole thrust of the allochthonous nappe; T 2 ) Middle Pleistocene (> 850 ka) normal faults of the Simeto Graben; T 3 ) Middle Pleistocene (600-200 ka) ramp thrust; T 4 ) reactivated portion of T 1 (< 40 ka). A c c e p t e d M a n u s c r i p A c c e p t e d M a n u s c r i p t

Fig. 1 ;

 1 Fig. 1; Hollestein et al., 2003) the measured NW-oriented GPS-velocities vs. Eurasia are parallel to

  Catania. The anticline is well exposed in the area of the village of San Giorgio, on the northern margin of the Catania Plain, where it deforms 60 ka-old terraced alluvials (T 6 in[START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF]) (Tab. 1) as well as several older (125 to 80 ka) marine terraced deposits (S 3 -S 5 in[START_REF] Catalano | The relationship between Late Quaternary deformation and volcanism of Mt. Etna (eastern Sicily): new evidence from the sedimentary substratum in the Catania region[END_REF]) (Tab. 1). The marine strandline of the OIS 5.5 (125 ka) is involved in the inner flank of the box-fold, along which is regularly increasing in elevation from 120 to 150 m a.s.l.. The 100 kaold (OIS 5.3) marine strandline crosses the crest of the structure, showing a change in elevation from about 80 m a.s.l., at the base of the inner flank, to 120 m a.s.l., at the crest of the anticline. The

  s c r i p t younger marine strandlines, as well as the 60 ka-old alluvial deposits, are confined to the forelimb of the anticline, being deformed as far as 2 km from the crest of the structure. The morphological features of different ages that are involved in the structure, as a whole, suggest a wavelength of the fold of about 4 km with an estimated amplitude of about 40 m.Towards the southwest, the crest of the anticline has been buried by the Holocene deposits of the Catania Plain. The stratigraphy of the uppermost 60 m thick levels of the alluvial plain have been investigated, by a series of electrical resistivity surveys, across the southwestern prolongation of the Sigonella Anticline, near the confluence between the Simeto and Dittaino rivers. North and south of the structure (cross-sections 1 and 3 in Fig.5), the reconstructed resistivity pseudocrosssections have imaged three well defined superimposed electrostratigraphical units, each corresponding to as many sedimentary units: the shallower 15-25 m thick level (R= 8 -4 Ohm•m) is represented by the Holocene alluvial sands and conglomerates, the intermediate 35 m thick level (R= 10 -30 Ohm•m) corresponds to the Late Pleistocene sands and conglomerates (T 7 ; 40 ka), while the deeper very low-resistivity (R<4 Ohm•m) layer corresponds to the top of the Middle Pleistocene marly-clay levels of the Simeto Graben. The tomographic model to the north of the fold crest (profile 1 in Fig.

  s c r i p t the two tectonic boundaries of the triangle zone and, thus, directly related to recent motion along the SPT and the CST.

  16) maximum horizontal compression in the Catania Plain area; (17) main earthquake (a): in the map the major (M≥7) historical earthquakes are reported; the inset also shows the main instrumental events, with available focal mechanisms (b)(the location of epicentres are from

  ) GPS vector; (19) rift-flank volcanoes (see in the inset a). The inset a shows the deformation belts connected to the active Nubia-Eurasia convergent boundary and the SCRZ respectively. The inset b refers to the relative motion between the plates.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Schematic geological profile across the eastern portion of the Sicily collision belt.

  The stereoplots (Schmidt net, lower hemisphere projection) in the inset refers to measured bedding of marine terraces of the following Oxygen Isotope Stages (OISs): a) OIS 7 (240-200 ka); b) OIS 5.5 (125 ka); c) OIS 5.3 (100 ka); d) OIS 5.1 (80 ka); e) OIS 3.3 (60 ka); f) OIS 3.1 (40 ka).

Fig. 4 :

 4 Fig. 4: Photo showing the Late Quaternary (< 80 ka) thrust ramp (CST) exposed in the Catania

Fig. 5 :

 5 Fig. 5: Geoelectric tomographies in the Catania Plain, across the southwestern prolongation of the

Fig. 6 :

 6 Fig. 6: Digital elevation model (DEM) of the topography across the Sigonella, Rotondella and

Fig. 7 :

 7 Fig. 7: Evolution of the sinuosity along the Simeto and the Dittaino rivers, across the Late

Fig. 8 :

 8 Fig. 8: Evolution of the Late Quaternary positive tectonic inversion at the flank of the SCRZ, in the

  Figure 1

  Figure 3

  Figure 4

  Figure 5

  Figure 6

  Figure 7

  Figure 8

Table 1

 1


Yannick Gorsse

Angelo Iollo

Haysam Telib

Lisl Weynans

A simple second order cartesian scheme for compressible Euler flows

We present a finite-volume scheme for compressible Euler flows where the grid is cartesian and it does not fit to the body. The scheme, based on the definition of an ad hoc Riemann problem at solid boundaries, is simple to implement and it is formally second order accurate. Error convergence rates with respect to several exact test cases are investigated and examples of flow solutions in one, two and three dimensions are presented.

Introduction

The computation of flows in complex geometries is a crucial issue to perform realistic simulations of physical or biological applications like for instance biolocomotion (fish swimming or insect flight), turbomachines, windmills... To this end several class of methods exist. Here we are concerned with immersed boundary methods, i.e., integration schemes where the grid does not fit the geometry. These methods have been widely developed in the last 15 years, though the first methods were designed earlier (see for example [START_REF] Peskin | The fluid dynamics of heart valves: experimental, theoretical and computational methods[END_REF], [START_REF] Berger | An adaptive cartesian mesh algorithm for the euler equations in arbitrary geometries[END_REF], [START_REF] Berger | Stable boundary conditions for cartesian grid calculations[END_REF]). The general idea behind immersed boundary methods is to take into account the boundary conditions by a modification of the equations to solve, either at the continuous level or at the discrete one, rather than by the use of an adapted mesh. The main advantages of using these approaches, compared to methods using body-conforming grids, are that they are easily parallelizable and allow the use of powerful line-iterative techniques. They also avoid to deal with grid generation and grid adaptation, a prohibitive task when the boundaries are moving.

A recent through review of immersed boundary methods is provided by Mittal and Iaccarino [START_REF] Mittal | Immersed boundary methods[END_REF]. They distinguish two different kinds of immersed boundary methods. In the first category are methods that handle the presence of boundaries by adding a continuous forcing directly to the flows equations, like in the original immersed boundary method ( [START_REF] Peskin | The fluid dynamics of heart valves: experimental, theoretical and computational methods[END_REF]). The boundary is represented by lagrangian markers exerting a forcing in the momentum equations. Also the penalty method ( [START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF], [START_REF] Arquis | Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide-milieux poreux: application la convection naturelle[END_REF]) and its recent developments ( [START_REF] Liu | A brinkman penalization method for compressible flows in complex geometries[END_REF], [START_REF] Boiron | A high-resolution penalization method for large mach number ows in the presence of obstacles[END_REF], [START_REF] Bergmann | Modeling and simulation of fish-like swimming[END_REF]), belong to this category: the solid is modeled as a porous medium with vanishing porosity, typically leading to a singular volume force in the momentum equations.

The methods in the second category are usually refered as "cartesian grid methods", where the forcing accounting for the presence of boundaries is performed at the discrete level. The ghost-cell approach, inspired by the ghost-fluid method developed by Fedkiw ( [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF], [14]) for multiphase flows belongs to this category. The first developments ( [START_REF] Majumdar | Rans solver with adaptive structured boundary non conforming grids[END_REF], [START_REF] Fadlun | Combined immersed boundary finite difference methods for three dimensional complex flow simulations[END_REF]) were followed by many other extensions [START_REF] Tseng | A ghost cell immersed boundary method for flow in complex geometry[END_REF], [START_REF] Ghias | A sharp interface immersed boundary method for compressible viscous flows[END_REF], [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF]. Ghost cells are cells in the solid with at least one neighbour in the fluid. The values on these ghost cells are extrapolated from the values in the fluid in order to impose the appropriate boundary conditions at the interface. The sub-mesh penalty method introduced by Sarthou et al. [START_REF] Sarthou | The sub-mesh penalty method[END_REF] is also related to this family. The immersed interface method is another approach, developed by Leveque and Li [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF] and extended to flow problems in [START_REF] Linnick | A high order interface method for simulating unsteady incompressible flows on irregular domains[END_REF]. It is based on Taylor expansions of the solution on each side of the interface. The "cutcell" approach belongs also to this class, and was mainly developed for compressible flows, see for example [START_REF] Pember | An adaptive cartesian grid method for unsteady compressible flows in irregular regions[END_REF], [START_REF] Ye | An accurate cartesian grid method for viscous incompressible flows with complex immersed boundaries[END_REF], [START_REF] Miller | A conservative three-dimensionnal eulerian method for coupled solid-fluid shock capturing[END_REF], [START_REF] Colella | A cartesian grid embedded boundary method for hyperbolic conservation laws[END_REF], [START_REF] Hu | A conservative interface method for compressible flows[END_REF]. These methods preserve conservation properties near the boundary.

Parallel to this classification, the issue of the accuracy of the immersed boundary methods, especially near the immersed interface, is worth being raised. For example, methods such as the original immersed boundary method from Peskin, or the penalty method, are known to be order one accurate in space. For Navier-Stokes equations, either compressible or incompressible, several types of methods provide an order two accuracy in space, global as well as local (near the interface): Tseng and Ferziger [START_REF] Tseng | A ghost cell immersed boundary method for flow in complex geometry[END_REF], testing their ghost-cell method on an uniform flow past a cylinder, report an order two accuracy in the L ∞ norm in the domain and on the boundary, provided that the values of the ghost cells are computed with a quadratic extrapolation. Mittal et al. [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF] and Ghias et al. [START_REF] Ghias | A sharp interface immersed boundary method for compressible viscous flows[END_REF], using a ghost-cell method, also obtained an order two accuracy for the same kind of test case. Ye et al [START_REF] Ye | An accurate cartesian grid method for viscous incompressible flows with complex immersed boundaries[END_REF], using a cut-cell method, reported an order two accuracy for the L 1 and L 2 norms of the error on the velocity in the whole domain for a two-dimensional Stokes flow. Kim et al. [START_REF] Kim | An immersed boundary finite volume method for simulations of flow in complex geometries[END_REF], with a finite volume approach, for decaying vortices obtained an order two accuracy for the L ∞ norm of the velocity error on the immersed boundary.

However, in the case of the compressible Euler equations, it seems more difficult to get order two accuracy near the boundary. In practice, only a few cartesian methods obtain local second order accuracy in the case of Euler flows.

In this sense, by investigating the semi-discrete approximation of hyperbolic and parabolic equations, Abarbanel et al. [START_REF] Abarbanel | On error bounds of finite difference approximations to partial differential equations -temporal behaviour and rate of convergence[END_REF] show that if the order of the space approximation at the boundary is lower than that in the interior domain, the overall order of convergence in the hyperbolic case will be lower than that corresponding to the parabolic case. Pember et al. [START_REF] Pember | An adaptive cartesian grid method for unsteady compressible flows in irregular regions[END_REF] proposed a cut cell method based on a Godunov algorithm and a correction step for cells at the boundary inspired by volume of fluid approaches, with a redistribution step to maintain conservativity. They tested the accuracy of their method on a Prantl-Meyer expansion wave. They reported for the entropy an order two of convergence for the L 1 norm on the whole domain, and an order one for the cells at the interface (and smaller convergence rates for the stagnation enthalpy). Forrer and Jeltsch [START_REF] Forrer | A high order boundary treatment for cartesian grid methods[END_REF], with a ghost-cell type method reported similar results on the same test case. Coirier and Powell [START_REF] Coirier | An accuracy assessment of cartesian mesh approaches for the euler equations[END_REF] had similar convergence results too for the Ringleb flow, with another cut cell method. Colella et al. [START_REF] Colella | A cartesian grid embedded boundary method for hyperbolic conservation laws[END_REF] obtained also the same kind of accuracy results for the case of a wave propagating in a straight circular channel, with another cutcell method generalizing the work of Pember et al [START_REF] Pember | An adaptive cartesian grid method for unsteady compressible flows in irregular regions[END_REF]. Hu et al. [START_REF] Hu | A conservative interface method for compressible flows[END_REF] devised a cut cell method based on a level-set description of the interface and a mixing procedure to merge small fluid cells with their neighbours. For a moving wall problem, they obtained an order two accuracy for the global error, and a slightly superior to one error on the boundary. Sjogreen et al. [START_REF] Sjogreen | A cartesian embedded boundary method for hyperbolic conservation laws[END_REF] presented a ghost-cell method based on a second order interpolation in the normal direction. They studied the conservativity properties of their scheme and showed a 1.6 order for mass loss. Chaudhuri et al. [START_REF] Chaudhuri | On the use of immersed boundary methods for shock/obstacle interactions[END_REF] coupled a 5th order WENO scheme to an immersed boundary method where the fluid variables at the ghost cells are obtained using a local quadratic interpolation. They showed accuracy between one and two for shock/obstacle interaction problems. Dadone and Grossman [START_REF] Dadone | Ghost-cell method for analysis of inviscid three dimensionnal flows on cartesian grids[END_REF] proposed a ghost-cell method based on an isoentropic pressure gradient correction at the boundary and obtained an order two accuracy for the normal velocity on the surface of a spheroid.

At this stage let us notice that ghost-cell methods used for compressible non-viscous flows are in fact more sophisticated than the ones developed for Navier-Stokes equations. For example, Forrer and Jeltsch [START_REF] Forrer | A high order boundary treatment for cartesian grid methods[END_REF] computed integrals of piecewise constant functions on polygons corresponding to the reflection of the computational fluid cells along the boundary. Dadone and Grossman [START_REF] Dadone | Ghost-cell method for analysis of inviscid three dimensionnal flows on cartesian grids[END_REF] used an isentropic model near the interface, and computed a corrected pressure value accounting for the curvature of the interface. In the context of cut-cell methods, the shape of the cut cells must be carefully identified. This task is non-trivial in two dimensions, and can become very tedious in three dimensions. Small cell effects can deteriorate the stability of the method, and strategies developed to cope with this problem increase the complexity of the programming ( [START_REF] Berger | An adaptive cartesian mesh algorithm for the euler equations in arbitrary geometries[END_REF], [START_REF] Berger | Stable boundary conditions for cartesian grid calculations[END_REF] and [START_REF] Berger | h-box methods for the approximation of hyperbolic conservation laws on irregular grids[END_REF]).

In this paper we present a simple globally second order scheme inspired by ghostcell approaches to solve compressible inviscid flows. In the fluid domain, away from the boundary, we use a classical finite-volume method based on an approximate Riemann solver. At the cells located on the boundary, we solve an ad hoc Riemann problem taking into account the relevant boundary condition for the convective fluxes by an appropriate definition of the contact discontinuity speed. To avoid pressure oscillations near the boundary, we weight the boundary condition with a fluid extrapolation, as a function of the angle between the normal to the boundary and the cartesian mesh. These ideas can be adapted to reach higher order accuracy. However, here we limit the discussion to second order accuracy as our objective is to device a method that can easily be implemented in existing codes and that is suitable for massive parallelization.

The paper is organized as follows: in section 2 we describe the finite volume scheme that we use to solve the flow equations in the fluid domain, away from the interface. In section 3 the scheme at the boundary is detailed. In section 4 we present numerical tests to validate the expected order of convergence.

Resolution in the fluid domain

We briefly describe how we solve the flow equations in the fluid domain, that is, on grid cells containing only fluid and whose neighbours contain only fluid too. The method consists in a classical finite-volume scheme based on a MUSCL ( [START_REF] Leer | Towards the ultimate conservative difference scheme i. the quest of monotonicity[END_REF]) reconstruction with slope limiters.

Governing equations

The compressible Euler equations are:

∂ρ ∂t + ∇ • ρu = 0 ∂ρu ∂t + ∇ • (ρu ⊗ u + pI) = 0 ∂E ∂t + ∇ • ((E + p)u) = 0
where the symbols have usual meaning. For a perfect gas

E = p γ -1 + 1 2 ρu 2 , p = ρRT and c = γRT 2.

Discretization in space

We focus for the moment on a two-dimensional setting. The extension to three dimensions of the interior domain scheme is straightforward. Let i and j be integers and consider the rectangular lattice generated by i and j, with spacing ∆ x and ∆ y in the x and y direction, respectively. Let W be the conservative variables and F x (W ), F y (W ) the convective flux vector in the x and y direction, respectively. By averaging the governing equations over any cell of the rectangular lattice we have

dW i j dt + 1 ∆ x (F x i+1/2 j -F x i-1/2 j ) + 1 ∆ y (F y i j+1/2 -F y i j-1/2 ) = 0
where W i j is the average value of the conservative variables on the cell considered, F x i+1/2 j the average numerical flux in the x direction taken on the right cell side, and similarly for the other sides.

The average numerical convective fluxes at cell interfaces are approximated here using the Osher numerical flux function [START_REF] Osher | Upwind difference schemes for hyperbolic systems of conservation laws[END_REF], so for example

F x i+1/2 j ≈ F x -+ W + W - A - x (W )dW (1)
where W -and W + are the values of W at the left and the right of the appropriate interface retrieved by a MUSCL ( [START_REF] Leer | Towards the ultimate conservative difference scheme i. the quest of monotonicity[END_REF]) reconstruction of the primitive variables,

F x -= F x (W -), A - x (W ) is the negative part of A x (W ) = ∂F x /∂W
and the integral is taken on a specific path assuming a locally isoentropic Riemann problem. The choice of a specific Riemann solver to compute the numerical flux function is inessential to the following developments.

Time integration

A second order Runge-Kutta scheme is used for the time integration. We denote by W the cell averages of the conservative unknowns at a given time stage. The right-hand side of the semi-discrete equations is a non-linear vector function of the conservative variable averages, -F(W). According to this notation, given W n , the second order Runge-Kutta scheme that we employ in the following is:

W (1) = W n -∆tF(W n ) W n+1 = W n -∆t 2 F(W n ) + F(W (1) ) (2)

A second order impermeability condition

The canonical finite-volume scheme described in the previous section is modified as follows for fluid cells that have one or more neighbouring cells separated by an interface.

The boundary condition on this interface derives from the impermeability assumption, i.e., given normal velocity to the boundary (zero for a steady wall, but non-zero for a moving wall). We are concerned with recovering second order accuracy on the impermeability condition, as for convection dominated flows on slender bodies, i.e. thin boundary layers, this is expected to significantly affect the computational results near the solid walls at high Reynolds numbers.

Level set method

In order to improve accuracy at the solid walls crossing the grid cells we need additional geometric information. This information, mainly the distance from the wall and the wall normal, is provided by the distance function. The level set method, introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF], is used to implicitly represent the interface of solid in the computational domain. We refer the interested reader to [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF], [START_REF] Sethian | implementation, and application of level set and fast marching methods for advancing fronts[END_REF] and [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] for recent reviews of this method. The zero isoline of the level set function represents the boundary Σ of the immersed body. The level set function is defined by:

ϕ(x) = dist Σ (x) outside of the solid -dist Σ (x) inside of the solid (3)
A useful property of the level set function is:

n(x) = ∇ϕ(x) (4)
where n(x) is the outward normal vector of the isoline of ϕ passing on x. In particular, this allows to compute the values of the normal to the interface, represented by the isoline ϕ = 0. The computation of the distance function follows state of the art numerical techniques [START_REF] Shan Jiang | Efficient implementation of weighted eno schemes[END_REF].

Euler equations in one dimension

To make the ideas clear, let us start from a simple one-dimensional configuration. The typical situation for a grid that does not fit the body is shown in Fig. 1. The plan is to modify the numerical flux at the cell interface nearest to the solid boundary, in order to impose the boundary condition at the actual fluid-solid interface location with a given accuracy. For a steady body, we impose u b = u(x b ) = 0 at the boundary point

x b where ϕ(x b ) = 0. i+1 i-1 i x b u i-1 u i s 1 s 2 s 3 φ i Figure 1:
Mesh points and zero level set. The red segment represents ϕ = 0; hence the interface lies between the center of cell i (fluid) and the center of cell i + 1 (solid). Slopes s 1 , s 2 and s 3 are used to compute the velocity at i + 1/2. Slope s 3 is only used when a slope limiter is applied.

Let u * be the contact discontinuity speed resulting from the solution of the Riemann problem defined at the interface between cell i and cell i + 1. We define a fictitious fluid state to the right of i + 1/2 such that u * takes into account, at the desired degree of accuracy, the boundary condition u b in x b . In particular, the value of u * at x i+1/2 is determined by the following formula:

u * = u b + 1 2 -d s b , (5)
where d = ϕ i ∆x , u b is the velocity of the obstacle (= 0 for a steady body) and s b is a slope defined by:

s b = d s 1 + (1 -d)s 2 , (6)
where

s 1 = u b -u i d and s 2 = u b -u i-1 1 + d
are the slopes represented in Fig. 1. This formulation has been chosen to avoid stability issues in the slope computation when

x b is close to x i since s b = u b -u i + 1 -d 1 + d (u b -u i-1 ) . (7)
If a slope limiter is needed, the limited slope is defined by s l b = minmod(s b , s 3 ), where s 3 = u i -u i-1 . An higher polynomial reconstruction to retrieve s b does not significantly modify the convergence rates shown in the following.

At x i+1/2 the left fluid state of the Riemann problem U -= (u -, p -, c -) is computed as usual with the MUSCL reconstruction. The right state is U + = (-u -+ 2u * , p -, c -), so that u * is the contact discontinuity speed of the resulting Riemann problem. The left and right states of the variables p and c are identical in order to impose the correct wave reflection.

The scheme is non conservative at the numerical interface x i+1/2 . However, the loss of conservativity at this point is negligible compared to all the other points and the shocks are correctly resolved as shown in the following.

Euler equations in two dimensions

In two dimensions the flow equations are solved by computing independently the flux in each direction. Therefore, we apply in each direction the method employed in one dimension in order to accurately enforce the boundary condition. When the level set function changes sign between two cells, we modify the numerical fluxes at the interface between these cells. In particular, we compute a value of the contact discontinuity velocity at the cell side that will take into account the boundary condition at a point belonging to the solid-fluid interface.

The interface point naturally appearing in this context is the intersection between the interface (ϕ = 0) and the segment connecting the two cell centers concerned by the sign variation of ϕ (for example the points A or B on Fig. 2). For the numerical flux computation, a fictitious state is created for instance between the cells (i, j) and (i + 1, j) on Fig. 2. This fictitious state will be used to compute the numerical flux at this cell interface. However, the boundary condition that we have to impose now is u A .n A = 0, where u A is the speed of the fluid at the boundary, and n A the outward normal vector of the body.

When n A is parallel to the cell side normal n cell we recover the one-dimensional case. However, when the scalar product between the normal to the physical boundary and the normal to the mesh side is close to zero, see Fig. 3 for example, we assume that the boundary condition only weakly affects the numerical flux, as it would happen for a body fitted mesh. The state corresponding to a cell in the solid is therefore modified compared to the one-dimensional case in order to take into account this requirement thanks to a convex interpolation.

We compute as before the left state primitive variables u -, p -, c -relative to the Riemann problem at the concerned numerical interface by a standard MUSCL reconstruction. According to Fig. 2 and considering the (i + 1/2, j)-flux, the right state will be U + = (-u -+ 2u w , v w , p w , c w ), where: The speeds u * and v * : The normal component of the contact discontinuity speed u * is calculated as in the one-dimensional case. With reference to Fig. 2, the level set function changes sign between x i,j and x i+1,j at point A. Let the physical boundary normal point to the fluid side. Thanks to the level set function ϕ, the normal of the level set is defined at each cell center. If we assume that the directions of the levelset normal vectors n i,j and n i+1,j are parallel, we find that the normalized distance between x i,j and A is

•      u w v w p w c w      = α      u * v * p - c -      + (1 -α)      u f v f p f c f      ; • α = n A .n cell ; • (u f , v f , p f , c f ) is an extrapolated fluid state detailed hereafter. and U -= (u -, v w , p -, c -). B A (i,j) (i+1,j) φ=0 n A n cell
d = |ϕ i,j | |ϕ i,j | + |ϕ i+1,j | (8)
and the normal vector n A is computed by

n A = n i,j + d (n i+1,j -n i,j ) (9)
where n i,j is a fourth order centered finite-difference approximation of ∇ϕ at point (i, j), or the exact normal if the geometry is analytically known. We determine the value of the contact discontinuity speed u * , relative to a Riemann problem defined in the direction normal to the cell side through x i+1/2,j , consistent at second order accuracy with u A • n A = 0 in A. The vector u * is determined as follows:

u * .n A = u * n = u A .n A + 1 2 -d s n A u * .τ A = u * τ = u -.τ A      ⇒ u * = u * n n x + u * τ τ x u * n n y + u * τ τ y
Where u A is the velocity of the obstacle (= 0 for a steady body), n A = (n x , n y ) t and τ A = (τ x , τ y ) t are respectively the normal and tangential vectors to the boundary at point A and the slope s n A is defined as in 1D:

s n A = u A .n A -u i .n A + 1 -d 1 + d (u A .n A -u i-1 .n A ) . ( 10
)
If a slope limiter is needed, the limited slope is defined by s n,l A = minmod(s n A , s n 3 ), where

s n 3 = u i .n A -u i-1 .n A .
Extrapolated fluid state: (u f , v f , p f , c f ) is computed as a linear extrapolation of the fluid variables to the point x i+1/2 j using the variables and the slopes of the closest upstream fluid cell. For example, considering the velocity direction is west to east, the extrapolated variables are determined with the variables and slopes of the closest upwind fluid cell with respect to u (starting cell of the arrow on Fig. 3), in order to respect the domain of dependence of the entropy equation. For a (i, j + 1/2)-flux, the same procedure is applied considering the closest upwind fluid cell with respect to v.

Euler equations in 3D

The method described in the previous sections is easily extended in three dimensions. The interior cells, having all neighbours in the fluid, are treated by the classical finitevolume method described in Section 2. For the numerical flux at the mesh interface across two points where the distance function changes sign, exactly the same procedure described in two dimensions is applied. First of all the intersection point A between the solid boundary and the segment connecting the cell centers concerned is found. At this point the normal is determined thanks to the level set function and u * n is computed thanks to a polynomial extrapolation in the direction orthogonal to the numerical interface considered and taking into account the boundary condition at A.

Two vectors τ 1 and τ 2 are defined in the tangent plane through A. As before, the speeds in the tangential directions u * τ 1 and u * τ 2 are extrapolated from the fluid side. The extrapolated fluid state is retrieved in the same as in two dimensions: we extrapolate starting from the closest upwind fluid cell. When several choices are possible, we consider the fluid cell whose velocity component in the extrapolation direction is higher.

Finally, the same correction based on the scalar product between the normal to the solid interface and the normal to the cell side is employed for the convex interpolation described before.

Numerical results

We present a set of numerical illustrations of increasing complexity ranging from onedimensional to non-symmetrical three-dimensional configurations. In particular, we compute convergence rates to reference solutions of well known test cases. We define the L p error as

N i |e i | p /N 1/p
where e i is the error relative to the i-th cell and N is the number of computational cells; the L ∞ error is defined as max(|e i |) over the grid. The convergence rates are computed as the linear regression of the error data.

Analysis of the accuracy in 1D

In the following, the accuracy of the method is evaluated on three test cases involving a reflection on a solid wall, as a function of the relative position of the reflecting wall and the computational boundary. The test cases concern a rarefaction wave, a continuous pulse and a shock wave. The exact solutions are represented on Fig. 4. We assume that the wall is located at

x b = x i + ψ(x i+1 -x i ), ψ ∈ [0, 1]
. The accuracy of the present method is compared with the case of a reflexion on a solid wall exactly located at i+1/2 (ψ = 1/2), where the boundary condition can be imposed without approximations. The exact solution for the reflection of the rarefaction wave is obtained using Riemann invariants. The velocity to the left of the rarefaction is -0.35. To the right, where the gas is steady, the pressure is 1 and the speed of sound is 1. Concerning the continuous pulse, the velocity inlet boundary condition is a time-dependent portion of a sinusoid of amplitude 0.01. The inlet pressure is 1. The reference solution is computed using 10000 discretization points, with the wall located at a grid interface (ψ = 1/2). The exact solution for the reflection of the shock is calculated with the Rankine-Hugoniot conditions and an exact Riemann problem at the wall. The shockwave has a speed corresponding to a Mach number of 1.22. The computations are initialized at the time corresponding to initial time in the figures and the errors are calculated at the time corresponding to final time, with respect to the exact or reference solutions (see Fig. 4).

The error in several norms versus ψ is computed on 128 grid points and is represented on Fig. 5. On the plots we compare the error of our scheme with the error of a standard scheme for which ψ = 1/2. This corresponds to exactly imposing the impermeability boundary condition in the computation of the numerical flux. We denote with "optimal" the error corresponding to this ideal situation.

Overall, for the three test cases, the errors in all norms are close to the optimal one. Let us first consider the continous solution cases. We opted for a numerical set up where the relative position of the initial condition and the wall is fixed. We remark that for increasing ψ the error slightly increases in L 1 and L 2 norms. This is coherent with the fact that we compute a boundary condition using a larger distance for extrapolation. The L ∞ norm is less regular because the grid point of the largest error may vary with ψ. It is however, on average, of the order of the optimal error. The case ψ = 1 corresponds to the case ψ = 0 and therefore the errors are identical.

For the shock case, the initial position of the discontinuity can be determined at most with an accuracy of the size of a grid cell at the discrete level. For this reason we opted for a wall location that varies with respect to the shock position as a function of ψ. Therefore the numerical solution corresponding to ψ = 0 does not exactly match to that obtained for ψ = 1. The error norms can be discontinuous for the same reason: negligible differences in the shock position can lead to finite errors as function of the relative position of the exact solution, the approximate solution and the discretization points. This is of course particularly evident in the L ∞ norm.

Ringleb flow

In this section, we will compare in some cases our results to two other existing methods that are comparable in terms of complexity of implementation. The first one is a simple symmetry technique, and the second one the ghost-cell CCST method, see [START_REF] Dadone | Ghost-cell method for analysis of inviscid three dimensionnal flows on cartesian grids[END_REF], that relies on a local isoentropic flow model at the wall. For both methods, two rows of ghost points are computed, and slopes are calculated independently in each region (fluid or solid). With respect to Fig. 6, the first step is to determine the value of all the variables U P = (u n P , u τ P , p P , c P ) at point P , u n P and u τ P being respectively the velocity component normal and tangential to the wall at point P . U P is determined with a bilinear interpolation of the four surrounding cell centers. Then, the value on the ghost point will be:

•      -u n P u τ P p P c P     
for the symmetry technique

•        -u n P u τ P + 2γ γ-1 p P ρ P -p G ρ G + u n P -u n G p P -ρ P u τ P R ∆n c P p G p P 1 γ       
for the ghost cell CCST technique called GCM in [START_REF] Dadone | Ghost-cell method for analysis of inviscid three dimensionnal flows on cartesian grids[END_REF]. R is the wall curvature and ∆n the distance between P and G.

The Ringleb flow refers to an exact solution of Euler equations. The solution is obtained with the hodograph method, see [START_REF] Shapiro | The Dynamics and Thermodynamics of Compressible Fluid Flow[END_REF].

The exact solution is formulated in (θ, V ) variables with u = V cos θ, v = V sin θ and V = u 2 + v 2 . The stream function is given by: Ψ = sin θ V . The streamlines equations are:

x = 1 2ρ 1 V 2 -2Ψ 2 + L 2 , y = sin θ cos θ ρV 2 with : L = - 1 2 ln 1 + c 1 -c - 1 c - 1 3c 3 - 1 5c 5 , c 2 = 1 - γ -1 2 V 2 , ρ = c 5 .
In our test case, the computational domain is [-0.5; -0.1]×[0; 0.6] and we numerically solve the flow between the streamlines Ψ 1 = 0.8 and Ψ 2 = 0.9. The inlet and outlet boundary condition are supersonic for y = 0 and y = 0.6 respectively. The convergence orders are calculated in L 2 and L ∞ norms on four different grids 32 × 48, 64 × 96, 128 × 192 and 256 × 384. We also perform a similar convergence study using a sequence of body fitted meshes and the same finite-volume scheme in the interior domain. The meshes are composed of streamlines and potential lines and an example is represented on Fig. 7.

The results for the L 2 and L ∞ norm are compared to the symmetry technique, the ghost-cell CCST method [START_REF] Dadone | Ghost-cell method for analysis of inviscid three dimensionnal flows on cartesian grids[END_REF] and the case with body-fitted mesh and standard wall boundary conditions in Fig. 8 and Fig. 9. The symmetry technique requires a special interpolation if one or more cell centers surrounding point P are in the solid part, see Fig. 6. To avoid any bias, at grid points not belonging to the fluid we take U P = U ex (x P , y P ), with U ex (x P , y P ) the exact solution of the Ringleb flow at point (x P , y P ). For the ghost cell CCST method, we proceed to a bilinear interpolation for all cases, the ghost points being a smooth extrapolation of the fluid region. The overall results show that the classical symmetry scheme is first order accurate in the L ∞ norm, even if the ghost points are based on the reflection of the exact solution. The other schemes have overall comparable accuracy, although the amplitude of the error is lower for the present scheme compared to the ghost-cell CCST method.

Subsonic circular cylinder

We consider a circular cylinder of radius 1, centered at (0, 0), in a computational domain of [-3, 3] × [0, 2]. The boundary conditions at the left and right borders correspond to an unperturbed flow at Mach 0.2; at the bottom and top boundaries we impose symmetry conditions. The finest grid is 1200 × 400. To compute the convergence rates, we have successively reduced by half the grid points in each space direction. The pressure and Mach number fields obtained on the finest grid are shown in Fig. 10(a). For pressure we plot 29 isocontours equally distributed in [0.89, 1.03]. The Mach number isocontours are 19 and they are equally distributed in [0., 0.45]. The pressure field is almost perfectly symmetric, whereas the Mach field shows the usual entropy layer. The convergence rates of the entropy error are presented in Fig. 10(b). They are even higher than expected thanks to the fact that the solution is approximately linear in space before and after the cylinder.

Transonic NACA0012 airfoil

The simulation of the flow over a transonic NACA0012 airfoil is performed. We choose a case with M ∞ = 0.8 and angle of attack α = 1.25 • to test our method. This test case is classical. However, The precise solution is difficult to obtain even though the overall pattern of the flow may be recovered with a low grid resolution (see a recent and complete discussion in [START_REF] Vassberg | In pursuit of grid convergence for twodimensional euler solutions[END_REF]). Here the computation is performed on a 30c × 20c domain, where c is the airfoil chord and with ∆ x = ∆ y = c Further resolution was needed for the lift and drag coefficients. The results obtained with ∆ x = ∆ y = c 400 in [-c, 2c] × [-c, c] are detailed and compared to [START_REF] Vassberg | In pursuit of grid convergence for twodimensional euler solutions[END_REF] on table 4.4. The reference results are obtained with state of the art schemes on body fitted orthogonal meshes with over 16 million points in the whole domain and more than 4000 points on the airfoil. To be remarked that the last significant digits of the drag coefficient are strongly influenced by the particular closure of the NACA0012 at the trailing edge.

Min [START_REF] Vassberg | In pursuit of grid convergence for twodimensional euler solutions[END_REF] Max [START_REF] Vassberg | In pursuit of grid convergence for twodimensional euler solutions[END_REF] Present

C l 0.

Blunt body

This example is relative to a Mach 2 flow around a circular cylinder of radius 1 centered in (0, 0). Because of the symmetry of the problem with respect to axis y = 0, the computational domain is restricted to [-3 , 0] × [0 , 6]. A reference simulation, relative to a polar body fitted grid, is performed to compare to the results obtained on the cartesian mesh. For the reference simulation the computational domain is r ∈ [START_REF] Abarbanel | On error bounds of finite difference approximations to partial differential equations -temporal behaviour and rate of convergence[END_REF][START_REF] Berger | Stable boundary conditions for cartesian grid calculations[END_REF] and θ ∈ [ π 2 , π]. In Fig. 12 the comparison between the results obtained on the cartesian grid are contrasted to those obtained using the polar grid. The shock position and the isolines are basically the same for both methods and coincide on the horizontal axis as the grid is refined.

Supersonic wedge

A more complicated shock flow pattern can be investigated considering a supersonic wedge. For given upstream Mach number, and below a critical wedge angle, a straight oblique shock develops at the wedge corner. In this configuration there exist a supersonic uniform flow solution downstream the shock. Using the present scheme with a 300 × 200 grid, the upstream Mach set at 4 and a 20 degrees wedge, we computed a shock angle of 32.49, whereas the exact value is 32.46.

Above a maximum flow deflection angle the oblique shock wave is no longer attached to the corner and is replaced by a detached bow shock. We consider a steady 25 degree wedge with an initial datum that consists in two discontinuous states connected through the Rankine-Huginot conditions for a shock moving with Mach=1.7. The grid is 1200 × 800. The discontinuity is placed ahead of the wedge angle. The speed of the shock and the wedge deflection are such that this configuration corresponds to a super-critical case. In Fig. 13 we present the Mach and pressure fields. we remark a contact discontinuity resulting from the bow shock interaction with the shock to the right. The pressure isocontours do not exhibit any perturbation at the contact discontinutity. These results are in good accordance with the literature (see for example Fig. 17.9, pag. 590 [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics -A practical introduction -2nd edition[END_REF]).

Conclusion

In this paper we have presented a cartesian method to solve compressible flows in complex domains with second order accuracy in smooth regions. Results show that the scheme is overall locally and globally second order accurate. This method is based on a classical finite volume approach, but the values used to compute the fluxes at the cell interfaces near the solid boundary are determined solving a Riemann problem where the left and right states are modified so to satisfy the boundary conditions with second order accuracy. The method is therefore non-conservative at solid boundaries, but shocks are shown to be correctly captured. This method is particularly simple to implement, as it doesn't require any special cell reconstruction at the solid-wall interface. The extension to three-dimensional cases is natural as the same procedure at the boundary is repeated in each direction. Forthcoming work will concern the extension of the present approach to multi-physics problems.

Figure 2 :

 2 Figure 2: Example of geometric configuration at the interface.A is the interface point located between (i, j) and (i + 1, j). The flux on cell interface (i + 1/2, j) is modified to enforce the boundary condition on A.

Figure 3 :

 3 Figure 3: Cells and stencils used for the extrapolation state, considering u * n > 0.

Figure 4 :

 4 Figure 4: (x, t) diagrams of exact solutions in terms of velocity isocontours for the three test cases.

Figure 5 :

 5 Figure 5: Present scheme error and optimal error.

Figure 6 :

 6 Figure 6: Illustration of the determination of ghost-cell variables.

  5 -0.45 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 (b) Body fitted mesh.

Figure 7 :

 7 Figure 7: Comparison between the cartesian mesh with the isoline ϕ = 0 and the body fitted mesh.

(

  cell CCST ; order 1.9 Present scheme ; order 2.0 Body fitted mesh ; order 2.0 cell CCST ; order 2.0 Present scheme ; order 2.0 Body fitted mesh ; order 2.0 (b) Sound speed

Figure 8 :

 8 Figure 8: Comparison of the L 2 accuracy of the present scheme with several methods. The convergence orders are detailed in the legend.

Figure 9 :

 9 Figure 9: Comparison of the L ∞ accuracy of the present scheme with several methods. The convergence orders are detailed in the legend.

  100 . A numerical zoom is then performed on [-c, 2c] × [-c, c] in order to achieve further resolution on the airfoil: the fully converged solution obtained on the 30c × 20c domain is interpolated on the boundary points of the [-c, 2c]×[-c, c] domain where ∆ x = ∆ y = c 200 . Pressure, mach, density contours and the pressure coefficient thus obtained are shown on Fig. 11.

  (a) Upper half plane: pressure; lower: mach

Figure 10 :

 10 Figure 10: Flow field and convergence rates for entropy.

Figure 11 :

 11 Figure 11: NACA0012, pressure, mach and density contours. Pressure coefficient on the airfoil surface compared to reference results.

  For pressure we plot 40 isocontours equally distributed in [0.97, 4.82]. The Mach number isocontours are 40 and they are equally distributed in [0., 0.91], the time is 0.318. In the Mach field

Figure 12 :

 12 Figure 12: 20 pressure isolines between 0.6 and 4, ∆ x = ∆ y = 1 20 on the left and 1 40 on the right. The results with the body fitted mesh are on upper half of the figure, and with the cartesian mesh on the lower half.

Figure 13 :

 13 Figure 13: Flow field past a super-critical wedge.

4. 7

 7 Mach 10 shock over three spheresThe computation of a planar shock reflecting over three spheres is performed. The spheres are located at (0,0,0), (-0.75, 1.4,-1.4), (0,-2,0) with radii 1, 0.4, 0.3, respec-tively. The size of the domain is [-2.5; 2.5] 3 . The numerical computation is performed on a 256 3 mesh. Four snapshots are shown on Fig.14. The interactions of the bow shocks give rise to complicated flow structures in the wake of the spheres.

Figure 14 :

 14 Figure 14: Mach 10 planar shock reflecting on a 3D sphere. 20 isopressure surfaces.

Table 1 :

 1 Drag and lift coefficient.

		348226045	0.356208937	0.34951
	C d	0.022453440	0.022684938	0.02056