Vincent Mouchi
email: vmouchi@gmail.com

Camille Godbillot

Catherine Dupont

Marc-Antoine Vella

Vianney Forest

Alexey Ulianov

Franck Lartaud

Marc De Rafélis

Laurent Emmanuel

Eric P Verrecchia

Provenance study of oyster shells by LA-ICP-MS

Provenance determination of archaeological remains is a valuable tool for reconstruction of past exchange networks. Among these materials, oyster shells are ubiquitous in sites from all ACCEPTED MANUSCRIPT 2 prehistorical and historical periods. Thus, they seem to be as promising candidates for provenance identification as they include chemical elements from the environment in their shells, which implies that an elemental fingerprint of the region of origin can be recorded in the shell composition. In this study, we present elemental measurements from 15 groups of modern and archaeological shells from 13 continental localities in mainland France and the island of Corsica (western Mediterranean Sea).

; therefore, they are not relevant to the geographical identification of the specimen's origin beyond a single region.

INTRODUCTION

Reconstruction of exchange networks in archaeology is paramount to study large-scale medium-to long-term economics, as well as cultural and technological transfers, since routes were vectors of ideas and cultural exchanges. Numerous studies have attempted to reconstruct the provenance of marbles [START_REF] Herz | Classical marble: Geochemistry, Technology, Trade[END_REF][START_REF] Matthews | The establishment of a data base of neutron activation analyses of white marble[END_REF][START_REF] Mrozek-Wysocka | Ancient marbles: Provenance determination by archaeometric study[END_REF], obsidian [START_REF] Robin | Identification and characterization of two new obsidian subsources in the Nemrut volcano (Eastern Anatolia, Turkey): The Sıcaksu and Kayacık obsidian[END_REF][START_REF] Kuzmin | Determination of the source for prehistoric obsidian artifacts from the lower reaches of Kolyma River, Northeastern Siberia, Russia, and its wider implications[END_REF], variscite [START_REF] Querré | Témoin d'échanges au Néolithique le long de la façade atlantique : la parure en variscite des tombes de l'ouest de la France[END_REF], ceramics [START_REF] Michelaki | On establishing ceramic chemical groups: exploring the influence of data analysis methods and the role of the elements chosen in analysis[END_REF], and ore [START_REF] Vogl | On the ore provenance of the Trojan silver artefacts[END_REF]. Yet, only a few studies have focused on shells and shell beads from chemical analyses [START_REF] Shackleton | Neolithic trade routes re-aligned by oxygen isotope analyses[END_REF][START_REF] Eerkens | Provenance analysis of Olivella biplicata shell beads from the California and Oregon Coast by stable isotope fingerprinting[END_REF].

In Mesolithic and Neolithic sites from the European Atlantic seaboard, food wastes composed of oyster shells are exclusively found on the coastline (Dupont, 2016) as the hunter-gatherer populations consumed the soft body and discarded the shell on site. More recent sites from Antiquity comprise oyster shells used for consumption in numerous localities, even distant from marine shores [START_REF] Bardot-Cambot | Consommer dans les campagnes de la Gaule Romaine[END_REF]. Indeed, during this period, oysters were consumed by the elite as delicacies Other techniques have been developed in order to identify the provenance of shell remains, such as geochemical analyses, not only for archaeology [START_REF] Eerkens | Provenance analysis of Olivella biplicata shell beads from the California and Oregon Coast by stable isotope fingerprinting[END_REF][START_REF] Eerkens | Stable isotope provenance analysis of Olivella shell beads from the Los Angeles Basin and San Nicolas Island[END_REF]Mouchi et al., 2020a), but also for quality control of modern cultured specimens prior to market distribution [START_REF] Bennion | Trace element fingerprinting of blue mussel (Mytilus edulis) shells and soft tissues successfully reveals harvesting locations[END_REF][START_REF] Morrison | Spatio-temporal trace element fingerprinting of king scallops (Pecten maximus) reveals harvesting period and location[END_REF]. In addition, geochemistry has been used to study connectivity of separated populations of mollusc species that colonize distant localities through larval dispersion [START_REF] Becker | Spatial and temporal variation in trace elemental fingerprints of mytilid mussel shells: A precursor to invertebrate larval tracking[END_REF][START_REF] Carson | Population connectivity of the Olympia oyster in Southern California[END_REF][START_REF] Fodrie | Connectivity clues from short-term variability in settlement and geochemical tags of mytilid mussels[END_REF][START_REF] Sorte | Elemental fingerprinting of mussel shells to predict population sources and redistribution potential in the Gulf of Maine[END_REF].

Molluscs build their shells in calcium carbonate by sampling chemical elements from their surrounding waters, and the shell's chemical composition reflects (at least partially) that of the environment [START_REF] Urey | Measurements of paleotemperatures and temperatures of the Upper Cretaceous of England, Denmark, and the southeastern United States[END_REF]. Most studies have used oxygen and carbon stable isotope ratios ( 18 O and  13 C, respectively) from shells as a provenance proxy [START_REF] Eerkens | Provenance analysis of Olivella biplicata shell beads from the California and Oregon Coast by stable isotope fingerprinting[END_REF][START_REF] Eerkens | Stable isotope provenance analysis of Olivella shell beads from the Los Angeles Basin and San Nicolas Island[END_REF][START_REF] Bajnóczi | Tracing the source of Late Neolithic Spondylus shell ornaments by stable isotope geochemistry and cathodoluminescence microscopy[END_REF][START_REF] Milano | Oxygen and carbon stable isotopes of Mytilus galloprovincialis Lamarck, 1819 shells as environmental and provenance proxies[END_REF][START_REF] Apolinarska | Can stable isotopes of carbon and oxygen be used to determine the origin of freshwater shells used in Neolithic ornaments from Central Europe[END_REF], but the use of this approach appears to be restricted to differentiate between marine and brackish/estuarine environments. Indeed,  18 O and  13 C are mainly representative of local biogeochemical conditions rather than the geographical location (temperature and salinity for  18 O and food source for  13 C; [START_REF] Craig | The measurement of oxygen isotope palaeotemperatures[END_REF][START_REF] Sharp | Principles of stable isotope geochemistry[END_REF]Lartaud et al., 2010a). Temporal homogeneity for the regional fingerprint cannot be insured as seasonal temperatures have an impact on  18 O, leading to the overlap of amplitude values with remote areas exhibiting similar seasonal temperature contrasts (Mouchi et al., 2020b). Additionally, it was recently demonstrated that oysters do not mineralize their shell in equilibrium with seawater during their first year of growth, making unlikely the use of  18 O as a proxy for local physicochemical conditions during this life stage [START_REF] Huyghe | Oxygen isotope disequilibrium in the juvenile portion of oyster shells biases seawater temperature reconstructions[END_REF]. Overall, elemental analyses tend to provide better results in provenance determination than carbon and oxygen isotopic compositions [START_REF] Zhao | Effects of sediment, seawater, and season on multi-element fingerprints of Manila clam (Ruditapes philippinarum) for authenticity identification[END_REF][START_REF] Ricardo | Spatio-temporal variability of trace elements fingerprints in cockle (Cerastoderma edule) shells and its relevance for tracing geographic origin[END_REF][START_REF] Bennion | Trace element fingerprinting of blue mussel (Mytilus edulis) shells and soft tissues successfully reveals harvesting locations[END_REF][START_REF] Morrison | Spatio-temporal trace element fingerprinting of king scallops (Pecten maximus) reveals harvesting period and location[END_REF]. Nevertheless, such studies have identified issues related to 'vital effects': the metabolism acts as a filter between the living environment and the shell, thus modifying the shell's final composition. It has indeed been shown (Mouchi et al., 2020a) that concentrations in rare earth elements (REE) can be used to identify the region of provenance for the cupped oyster Crassostrea gigas. Still, these elements are not discriminant for the flat oyster, Ostrea edulis, which is the only oyster species found in European archaeological sites, as the cupped oyster Crassostrea angulata was imported to Portugal from Asia during the 17 th century. Consequently, other chemical elements need to be investigated as provenance proxies for O. edulis.

In the present study, we use a new multi-elemental approach based on laser-ablation inductively-coupled plasma mass spectrometry (LA-ICP-MS) performed on 68 modern and archaeological oyster specimens (C. gigas and O. edulis) from 15 groups of different localities and ages. We aim at using this approach to identify potential discriminant elements for provenance studies applied to both oyster species found in Europe. This technique was selected as it is fast and commonly used in carbonate geochemistry studies [START_REF] Durham | Rapid determination of oyster lifespans and growth rates using LA-ICP-MS line scans of shell Mg/Ca ratios[END_REF].

MATERIALS AND METHODS

Collection sites and specimens

Samples were collected from multiple modern localities and archaeological sites in continental France and the island of Corsica (Figure 1; Table 1). As elemental composition of the rivers (suggested to influence the composition of seawater and shells at the output) reflects the composition 103 of substrate weathered on land, a brief description of the geological formations of watersheds are 104 indicated with site descriptions in Supplementary Information 1. Specimens are divided into groups 105 of origin, according to their species ('Cg' and 'Oe' for C. gigas and O. edulis, respectively), their 106 locality (see Table 1), and their age ('M' for modern and 'A' for ancient specimens). Although the 107 species Crassostrea gigas has been recently renamed Magallana gigas [START_REF] Salvi | Molecular taxonomy in 2D: a novel ITS2 rRNA sequence structure approach guides the description of the oysters' subfamily Saccostreinae and the genus Magallana (Bivalvia: Ostreidae)[END_REF], 108 this genus change is debated by numerous experts [START_REF] Bayne | The proposed dropping of the genus Crassostrea for all Pacific cupped oysters and its replacement by a new genus Magallana: A dissenting view[END_REF], so we will continue to use the 109 name Crassostrea gigas throughout the paper. 110 [START_REF] Dechambenoy | Apport de la thermographie infrarouge aéroportée à la connaissance de la dynamique superficielle des estuaires (système Charente-Seudre-Anse de l'Aiguillon)[END_REF]. 5: Tès modern site. 6: Leucate modern site. 7: Corsica island sites, comprising the modern site of Diana lagoon (A) and the Antiquity site of Aléria (B). Additional archaeological groups are the Medieval shells found at La Malène of confirmed Mediterranean origin and two groups unearthed at Cybèle (Lyons) of unknown origin.

Sample selection and preparation

The preparation protocol is given in Mouchi et al. (2020a). Multiple specimens were selected from each group based on the apparent preservation of the umbo region of the shell. Epibionts, when present, were mechanically removed. The umbo of the left valve was embedded in Huntsman Araldite 2020 epoxy resin, cut through using a Buehler low-speed saw to expose the preserved internal structures along the growth axis, and made into polished 750 µm-thick sections. As the umbo records the entire growth of the shell in a condensed area, geochemical analyses on such sections have the benefit of being deprived of surface contaminations of chemical or organic origin [START_REF] Richardson | The age determination and growth rate of the European flat oyster, Ostrea edulis, in British waters determined from acetate peels of umbo growth lines[END_REF][START_REF] Kirby | Stable isotope sclerochronology of Pleistocene and recent oyster shells (Crassostrea virginica)[END_REF]. Indeed, the umbo is a dense area, protected from the exterior environment within the bivalve, which generally presents a much better preservation compared to the rest of the shell. In addition, oyster shells are made of lamellae that can be intertwined with empty spaces that can be filled with organic and sedimentary contaminants, which could strongly interfere with the quality of the measurements.

Each section was observed under cathodoluminescence to ensure the pristine state of the surface and to outline the seasonal calibration of the shell based on the intensity of luminescence (Lartaud et al., 2010b;[START_REF] Mouchi | Reconstructing environments of collection site from archaeological bivalve shells: case study from oysters (Lyon, France)[END_REF]Huygue et al., 2019). Following the growth direction, areas of strong luminescence correspond to parts of the shell built in summer periods, while areas in dull luminescence have grown during winter periods. Observations were completed using a Cathodyne-OPEA cold cathode at ISTeP, Sorbonne Université (Paris, France) with 15-20 kV, 200-400 µA.mm -2 and a pressure of 0.05 Torr. Shells from Leucate (Cg_Leu_M) were the only specimens that did not exhibit cathodoluminescence, which made the seasonal calibration impossible. The seasonal calibration from morphology proposed by [START_REF] Kirby | Stable isotope sclerochronology of Pleistocene and recent oyster shells (Crassostrea virginica)[END_REF] was also unsuccessful due to the absence of the necessary curved surface of the umbo. All specimens from the other groups were seasonally calibrated.

Geochemical analyses by LA-ICP-MS

Measurements were performed using an Element XR (ThermoScientific) ICP-MS coupled with a RESOlution 193 nm ArF excimer ablation system equipped with an S155 two-volume ablation cell (Australian Scientific Instruments) at ISTE (University of Lausanne, Switzerland). The analyses were carried out with a pulse repetition rate of 20 Hz and an on-sample energy density of 4 J.cm -2 .

Sample areas were cleaned from potential surface contamination (during polishing) by conducting preablation. Spot size was 200 µm. Measurements were conducted on both summer and winter zones, based on the seasonal calibration. Measurements on Leucate specimens (Cg_Leu_M) were positioned randomly instead. Each specimen (shell) was measured several times in order to avoid bias due to internal variability. We carefully avoided performing measurements in suspicious areas regarding preservation, based on cathodoluminescence observation. We performed at least two adjacent measurements for each selected season area, with generally six to eight measurements per shell (sometimes less, e.g. four measurements for some Oe_MO_M specimens, due to lower growth rates restricting the surface available for measurements). For most specimens, three to four successive season areas were sampled. Prior and following each 15-sample series, NIST SRM 612 was measured for external calibration. Accuracy of the analyses was checked against measurements of the BCR-2 basalt reference material from USGS according to the GeoReM preferred values [START_REF] Jochum | GeoReM: A new Geochemical database for reference materials and isotopic standards[END_REF].

Relative standard deviation for all measured elements range from 1.2 % to 5.8 % (Supplementary Information 2). Measured elements were La, Ce, Pr, Nd, Sm, Ti, Cr, Cu, Zn, Sr, Ba, Pb, U, Y, and Ca as the internal standard. Data reduction was performed using the LAMTRACE software [START_REF] Jackson | LAMTRACE data reduction software for LA-ICP-MS. In: Laser ablation ICP-MS in the Earth sciences: Current practices and outstanding issues[END_REF]. A total of 472 measurements were acquired.

Data processing

All data processing was conducted using the Matlab software (MathWorks, www.mathworks.com, v. R2017a). As none of the 14 measured elements presented a normal distribution (Kolmogorov-Smirnov test) due to a right-skewed distribution (towards larger element abundances), it was necessary to transform the data to normality for statistical processing, which was performed by using the cubic root transform [START_REF] Chen | Power transformations to induce normality and their applications[END_REF]. Hierarchical cluster analyses were used to estimate seasonal differences (using the seasonal age models based on the [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF] with the exact Euclidean method. This method allows high-dimensional data points to be imbedded in low dimensions (typically two dimensions) in a way that respects similarities between points: distant points in a high-dimensional space correspond to distant embedded lowdimensional points, and nearby data points in the high-dimensional space correspond to nearby embedded low-dimensional points (MathWorks, www.mathworks.com, v. R2017a). In other words, each point represented by t-SNE corresponds to the overall composition of one measurement of all 14 elements, two close points having a similar composition while two distant points reflect different compositions. To prevent the undesirable over-influence of some elements showing wide variance, we performed t-SNE on normalized data. The Matlab script for t-SNE processing is available as Supplementary Information 3. As we focus on comparing the concentrations between individual shells from multiple localities using the same (LA-ICP-MS) analytical method, the host phase (organic or mineral) of specific elements and their resulting co-variations within a measurement are not discussed here.

RESULTS

All measurement results are presented in Figure 2a. The three elements with the largest coefficients of variation are Ce, La and Zn, followed by U, Ba, and Ti (Fig. 2b). All other elements appear to have a low variance in the dataset. Except for Sr, Ba, Zn, La, and Ce, all measured elements have concentrations below 1 ppm. Boxplots by element for all the shell groups are provided as Supplementary Information 4. None of the elements can be used to differentiate between the origins of the various groups on its own. Although some O. edulis groups are well defined without any measurement overlap from other groups (i.e., Oe_Ale_A and Oe_Lan_A), some groups of this species are more difficult to discriminate geochemically (Figure 4a). Indeed, when only considering groups originating from localities on the English Channel and the Atlantic Ocean coastlines (i.e., open seas with strong tidal influence), numerous measurements cluster without discrimination between the locality of origin (dashed ellipse in Figure 4b). Only the Landévennec (Oe_Lan_A) group is absent from this common cluster.

Inversely, the Oe_BdV_M group measurements are found solely in this common cluster. For the Oe_RsB_M, Oe_MO_M, Oe_Bro_A and Oe_BrT_A groups, other measurements are found in 'secondary' clusters in addition to their presence in the common cluster, which are specific to the group. Some specimens from these groups have a unique fingerprint corresponding to one or the other cluster (common or secondary), while other specimens have fingerprints corresponding to both clusters along successive measurements, following shell growth (Figure 5). Amongst them, Oe_MO_M and Oe_Bro_A (both from Charente-Maritime) share this 'secondary cluster', while Broue tower (Oe_BrT_A; from the Charente-Maritime region as well) has a 'secondary cluster' isolated from the others (Figure 4c). It is also worth noting that specimens from the Oe_BrT_A group, which have been collected from two distinct temporal units in the site (see Supplementary Information 1), present measurements in the common cluster as well as in the secondary cluster for specimens from both stratigraphic units (sometimes along the same specimen), confirming that the local elemental fingerprint is identical over centuries. However, the two specimens from Landévennec (Oe_Lan_A), belonging to two different temporal periods, are discriminated from one another by t-SNE (Figure 3)

in the form of two groups of measurements in the Landévennec cluster. When considering the two remaining groups of confirmed origin, both from the Mediterranean Sea (an enclosed basin with a narrow opening to the Atlantic Ocean through the Strait of Gibraltar), each group is defined by a single cluster without any overlap (in green and orange in Figure 4d). The measurements from the groups of unknown origin overlap at least partially with the Oe_Mal_A group.

Data collected from specimens bred at Baie des Veys and Marennes-Oléron were used to study the differences between seasons for both species based on dendrograms (Supplementary Information 5). Seasons cannot be reliably discriminated, as indicated by the cluster analyses that show both seasonal measurements mixed together. For Cg_BdV_M, the cophenetic correlation coefficients are 0.80 and 0.52 for average and Ward methods, respectively. In the same site, for Oe_BdV_M, coefficients are 0.89 and 0.77 for average and Ward methods, respectively. Similar cophenetic correlation coefficients are found for Cg_MO_M with 0.70 and 0.62 for average and Ward methods, respectively, and for Oe_MO_M with 0.89 and 0.77 for average and Ward methods, respectively. The identification of seasons for all the measurements (except Leucate specimens; see section 2.2) in the t-SNE distribution is also reported in Supplementary Information 5.

DISCUSSION

Large-scale provenance discrimination

It was previously reported that REE and Y cannot be used for O. edulis as a provenance proxy due to vital effects, unlike C. gigas (Mouchi et al., 2020a). Consequently, other elements were investigated in this study in order to provide a potential regional fingerprint for this species and to find a way to efficiently compare measurements from O. edulis and C. gigas in terms of sourcing. The results show that REE are not the only elements whose incorporation differ between these two species.

Indeed, the Cg_BdV_M and Oe_BdV_M groups, which originate from the same modern locality (and reared simultaneously on the same culture table in Baie des Veys, Normandy), exhibit obvious different compositions (Figure 3). Observations are similar for Cg_MO_M and Oe_MO_M from Marennes-Oléron, Charente-Maritime. However, it seems possible, when considering only one of the two species, to discriminate the composition depending on their region of origin, if not the locality.

For C. gigas specimens, the shells from the Mediterranean Sea (Cg_Leu_M and Cg_Dia_M), seem to be separated by t-SNE from the Atlantic Ocean/English Channel groups along the second dimension (approximately around unit 25 in Figure 3). This separation appears to be controlled by a lower abundance in light REE in the Mediterranean groups compared to the Atlantic groups (Supplementary Information 4). This low abundance in light REE can either be specific to the Mediterranean shells or to the relatively small watersheds from these two localities (the smallest of the entire dataset). Although most of the Baie des Veys measurements are discriminated by t-SNE from other groups, some measurements share their composition with a part of the Tès group (Cg_Tes_M), which itself is partially connected to the Marennes-Oléron group (Cg_MO_M). This observation tends to indicate that specimens from the English Channel (Baie des Veys) and the Atlantic Ocean (Marennes-Oléron and Tès) exhibit some similarities in their geochemical composition.

For O. edulis, the common cluster (dashed ellipse in Figure 4b) gathers measurements from specimens from the Atlantic Ocean and English Channel coastlines, covering roughly 700 km of coastline (Figure 1) from Charente-Maritime in the south (Oe_MO_M, Oe_Bro_A, Oe_BrT_A) to Normandy in the north (Oe_BdV_M). As the English Channel and the Atlantic Ocean are connected across a wide area, the overall composition of seawater should share some similarities. On the contrary, the Mediterranean Sea is a semi-enclosed basin and is recharged by the Atlantic Ocean through the narrow Strait of Gibraltar. A difference in composition is expected in shells from Mediterranean localities. In addition, the near absence of tidal range in this region compared to the Atlantic Ocean can cause differences in fluid mixing and elemental incorporation in the shells.

Therefore, it appears possible to discriminate, at a large scale, between specimens from the Atlantic Ocean/English Channel and the Mediterranean Sea.

Fine-scale provenance discrimination

The common cluster composed of samples from Oe_BdV_M, Oe_RsB_M, Oe_MO_M, Oe_Bro_A, and Oe_BrT_A shells (Baie des Veys: Normandy; Riec-sur-Belon: Brittany; Marennes-Oléron, Brouage and Broue tower: Charente) represents a homogenous fingerprint. These groups are from watersheds with quite distinct geological formations (limestones, basalts, etc…) but share a large seawater body (i.e., the Atlantic Ocean). However, we observe that each group in this cluster has also a secondary cluster outside the common one, which seems group-specific, except for the Baie des Veys (Oe_BdV_M; Figure 4b). Within these groups, some specimens have all their measurements within a single cluster (either the common cluster or the secondary cluster of the corresponding group), while other specimens have some measurements in both clusters (Figure 5). As these specimens have not been transferred from one location to another during growth, it appears that these changes in elemental incorporation occurred over the growth of the organism in the same site. As these changes occur multiple times over the growth period (Figure 5), it cannot be caused by ontogeny. This cyclical pattern is not linked to seasonality (this was checked by the seasonal calibration using cathodoluminescence; see t-SNE figure in Supplementary Information 5).

Consequently, a different factor is at play. Temporal resolution of the LA-ICP-MS measurements is approximately one month, based on a tidal cathodoluminescence signal [START_REF] Mouchi | Chemical labelling of oyster shells used for time-calibrated high resolution Mg/Ca ratios: a tool for past estimation of seasonal temperature variations[END_REF]; see monthly cycles on Figure 5b) and in vivo labelling (Lartaud et al., 2010b; Supplementary Information 6). Growth rates can however be different according to the specimen and the period considered, and some measurements may represent a slightly shorter or longer period [START_REF] Huyghe | New insights into oyster high-resolution hinge growth patterns[END_REF]. Since all the corresponding groups originate from the English Channel and the Atlantic Ocean shorelines, tidal influence can be a pertinent factor. Indeed, all the groups from the Mediterranean Sea (Oe_Mal_A, Oe_Ale_A), which exhibits a very weak tidal range (only a few centimetres), are excluded from the common cluster (Figure 4d). In localities with a broad tidal range, the influence of the chemistry of nearby rivers could fluctuate widely, increasing during low tide and, on the contrary, being diluted with seawater at high tide. If this were the case, the secondary clusters would represent the correct local fingerprint, while the common cluster would correspond to the fingerprint of Atlantic seawater. Manganese is known to fluctuate in oyster shells from ebb to flood currents [START_REF] Huyghe | New insights into oyster high-resolution hinge growth patterns[END_REF]. Other chemical elements (e.g. those measured here) could present similar fluctuations.

Switching between these two compositions, which are responsible for the presence of two distinct fingerprints, could be caused by the relative proportion (by volume in the analysed shell) of carbonate formed during ebb and flow currents. The proportion may vary, depending on growth rates. The case of the Baie des Veys shells (Oe_BdV_M group) is singular as it has no secondary cluster (Figure 4b).

Interestingly, this site is also the locality with the most open marine settings from all modern groups.

Although some shells of other groups did not have measurements corresponding to secondary clusters, the fact that no measurements from the 38 performed on five specimens are outside this common cluster would tend to imply that there is only one fingerprint in this site. This observation is supported by the marine water mass current of the Baie des Veys, which flows towards the direction of the rivers' outputs from the locality where the oysters were collected (Figure 1 This assumption is supported by a higher and more stable salinity in Baie des Veys (approximately 33 all year) than in Marennes-Oléron (27 to 34; Lartaud et al., 2010c).

Temporal changes in the local elemental fingerprint have already been reported in mollusc shells, such as Cerastoderma edule and Mytilus edulis from Ireland [START_REF] Ricardo | Spatio-temporal variability of trace elements fingerprints in cockle (Cerastoderma edule) shells and its relevance for tracing geographic origin[END_REF]Bennion et al., 2019). These differences occur within months and observations in this study confirm that it is also the case for O. edulis, although the same is not true for C. gigas. These differences appear to switch back and forth between two fingerprints (Figure 5), one being common to multiple groups of shells of various origins in the French Atlantic Ocean and English Channel coastlines.

The case of the Charente-Maritime groups (Marennes-Oléron: Oe_MO_M; Brouage:

Oe_Bro_A; and Broue tower: Oe_BrT_A) is interesting, as two distinct secondary clusters are visible in the t-SNE plot (Figure 4c). One of these groups of shells has been collected in the present day (Marennes-Oléron) with the exact location known, while the two other groups (Brouage and Broue tower) have been found in archaeological sites on land, with evidence that these locations were on the shoreline at the corresponding ages of the specimens. The fingerprint of the modern group is shared by the Brouage group dated from the 16 th century CE (Oe_Bro_A). On the contrary, the Broue tower shells from the 11 th and the 14 th centuries CE (Oe_BrT_A) present a distinct fingerprint. It may be possible that the oysters from this group were living in a locality in the region under the influence of another freshwater flow with a different water composition. Indeed, the Marennes-Oléron bay has two major tributary rivers, which cross different geological formations, although the main rock type is limestone [START_REF] Bourgueil | Notice explicative, Carte géol[END_REF][START_REF] Bourgueil | Notice explicative, Carte géol[END_REF][START_REF] Platel | Notice explicative, Carte géol[END_REF][START_REF] Platel | Notice explicative, Carte géol[END_REF][START_REF] Bambier | Notice explicative, Carte géol[END_REF][START_REF] Hanztpergue | Notice explicative, Carte géol[END_REF][START_REF] Mourier | Notice explicative, Carte géol[END_REF]. The site of Brouage, from which were unearthed the shells of the Oe_Bro_A group, is located at the same level in the bay as the modern Marennes-Oléron group (Oe_MO_M; [START_REF] Champagne | La citadelle de Brouage et la dynamique paléoenvironnementale du marais charentais : l'apport de la malacologie et de la palynologie[END_REF]; Figure 1), with which it shares its secondary cluster. The site of the Broue tower, where the shells from the Oe_BrT_A group were discovered, is located further south, near the current output of the Seudre river (Figure 1; [START_REF] Normand | Broue (Saint-Sornin -Charente-Maritime) : un site élitaire au coeur des marais charentais. Rapport intermédiaire de fouille programmée triennalecampagne 2019[END_REF]. As seawater currents cross the bay from north to south, the oysters from the Oe_MO_M and Oe_Bro_A groups must have been under the influence of the Charente River (north side of the bay), but not the Seudre River. On the contrary, specimens from the Oe_BrT_A group must have been mainly under the influence of the Seudre River. The measurements show that, in such contexts, the identification of the region can be narrowed down to the watershed level. In addition to this, although temporal differences in local elemental fingerprints have been reported in other regions [START_REF] Ricardo | Spatio-temporal variability of trace elements fingerprints in cockle (Cerastoderma edule) shells and its relevance for tracing geographic origin[END_REF][START_REF] Bennion | Trace element fingerprinting of blue mussel (Mytilus edulis) shells and soft tissues successfully reveals harvesting locations[END_REF], the fact that both modern (Oe_MO_M) and ancient (Renaissance, Oe_Bro_A) shells share a common cluster and a local cluster over the course of several centuries indicate that multiple measurements over the growth direction of shells can allow the local elemental fingerprint to be recognized, whatever the period of collection. This interpretation is also supported by the Oe_BrT_A specimens, which present fingerprints corresponding to both the common cluster and the same (locality-specific) secondary cluster, whether specimens are from the 11 th or the 14 th c. CE.

The only group of shells from the English Channel/Atlantic Ocean coastline that have no measurements corresponding to the common cluster is Oe_Lan_A. The Medieval site at Landévennec, from which these shells were unearthed, is located in the direct proximity of the Aulne River (Figure 1). Although the salinity is close to marine (generally around 32 except during flooding events where salinity drops to 15; [START_REF] Monbet | Bilan des connaissances océanographiques en rade de Brest. Rapport technique[END_REF], the river rate of flow (annual mean of 24 m 3 .s -1 ; [START_REF] Auffret | Dynamique sédimentaire de la marge continentale celtique[END_REF] For "Atlantic Ocean" C. gigas specimens, the Baie des Veys group (Cg_BdV_M; although situated in the English Channel, with a wide access to the Ocean) appears to be the most distinct by t-SNE between these groups, with only some measurements overlapping with the Tès group. This is interesting, as O. edulis shells from the same locality (Oe_BdV_M) had the peculiar characteristic of representing the composition of a cluster common to other Atlantic groups. The Tès (Cg_Tes_M) and Marennes-Oléron (Cg_MO_M) groups share their distribution in the t-SNE graph, with the Marennes-Oléron distribution entirely comprised into a restricted portion of the Tès distribution. The large distribution by t-SNE of Tès measurements compared to Marennes-Oléron can be explained by the annual range of salinity (approx. 15 and 5, respectively for Tès and Marennes-Oléron; Lartaud et al., 2010c). These two localities (Arcachon Basin and Charente-Maritime) have their watersheds crossing different geological formations with mainly Cenozoic river deposits for Tès [START_REF] Dubreuilh | Could occupation duration be related to the diversity of faunal remains in Mesolithic shell middens along the European Atlantic seaboard?[END_REF] and mainly Mesozoic limestone for Marennes-Oléron (Bourgueil and[START_REF] Bourgueil | Notice explicative, Carte géol[END_REF][START_REF] Bourgueil | Notice explicative, Carte géol[END_REF][START_REF] Platel | Notice explicative, Carte géol[END_REF][START_REF] Platel | Notice explicative, Carte géol[END_REF]. Still, the deposit cover of the respective watersheds does not appear to be sufficient to entirely discriminate these groups.

Attempting to identify the provenance of the Lyons/Cybèle shells

The Oe_Mal_A group from La Malène belongs to a specific locality, which is unknown, but at least confirmed from the Mediterranean coastline, probably from Languedoc coastline due to the associated remains from other taxa specific to Roman-influenced Languedoc customs (Bardot-Cambot and [START_REF] Forest | Les paysans consommateurs et les vestiges fauniques archéologiques : une fausse évidence[END_REF][START_REF] Forest | Les paysans consommateurs et les vestiges fauniques archéologiques : une fausse évidence[END_REF][START_REF] Forest | Étude archéozoologique préliminaire : conchyliologie. La Malène -Les Piboulèdes (Lozère)[END_REF][START_REF] Mureau | Consommation et exploitation des ressources animales dans l'est du Massif central et le Languedoc de la fin de l'Antiquité tardive au haut Moyen Âge[END_REF]. Measurements from this group are discriminated from the other O. edulis groups of known origins (all from the Atlantic Ocean/English Channel coastlines) and these groups are relatively close to one another. We note that the Oe_X1_A and Oe_X2_A group distributions are partially overlapping (Figure 4d). It seems that the shells from these groups of unknown origin, and found at the Antiquity Cybèle site (Lyons, over 200 km away from the nearest shore), share part of their fingerprint, which could imply that all these specimens lived in the same geographic region, yet not the exact same locality. Therefore, although geochemistry confirms that the Cybèle groups initially determined by morphometric parameters (Bardot-Cambot, 2013) do exist, our data tend to discard the hypothesis that one group originated from the Mediterranean Sea and the other from the Atlantic Ocean. Although additional sites from both Atlantic and Mediterranean shorelines would be required to affirm this, it seems here that populations from these groups lived a short distance away from each other, and the fingerprint of the Oe_X1_A group seems to correspond to that of La Malène specimens (Oe_Mal_A). If true, additional specimens in other localities on the Mediterranean shoreline should be analysed to identify their locality of origin, but at the moment, we were not able to collect wild modern specimens.

Quality of the discrimination method

As our dataset is constituted of 14 elements per measurement, attempts were made to reduce this number for ease of analytical and data processing work. Some elements were found to be of increased importance compared to others to adequately separate measurements from different groups using t-SNE. By successive removal of one element from the dataset, Pb, Ba and U are the most effective at successfully discriminating O. edulis groups from the Mediterranean Sea and the Atlantic Ocean (Supplementary Information 7). Rare Earth Elements (La, Ce, Pr, Nd, Sm), as found by Mouchi et al. (2020a), were most useful for discriminating C. gigas groups. All other elements were found to contribute equally to discriminating specific localities. In the end, the removal of any of the measured elements systematically reduces the quality of the discriminating method.

Cross-validation of the clusters was attempted on O. edulis measurements (without the groups of unknown origin, i.e. Oe_X1_A and Oe_X2_A) using discriminant analysis. We obtained a crossvalidated classification error of 24.3% (out of 210 measurements). However, when we removed all measurements from the common cluster (putatively corresponding to the open marine fingerprint, including all measurements from the Oe_BdV_M group) from the comparative dataset, we obtained a cross-validated classification error of 0.8% (out of 123 measurements from 7 groups). Becker et al. in bays present some similarities to the regional open marine fingerprint. Our observations concur with these conclusions, albeit with additional precision.

CONCLUSION

Previous attempts to discover the geographical origin of archaeological oyster shells based on geochemistry presented various degrees of success due to the use of isotopic ratios that are strongly influenced by local environmental conditions, which can be similar over extremely remote areas. In this study, elemental fingerprinting of oyster shells from 15 different groups was used, highlighting several points. First, different elemental incorporation behaviours prevent the direct comparison of different species. Second, cluster analysis from the t-SNE method suggests that a regional fingerprint seems to exist and to be distinguishable from large water bodies (seas and oceans). Third, in environments with restricted influence of oceanic water, a more precise fingerprint is identified and the locality (at the watershed scale) can be found. In addition, as previous studies reported, none of the elements is able to differentiate between the origins of various groups on its own. Last but not least, this technique was successful to identify the provenance of oyster shells of known and previously unknown origin from archaeological sites. Therefore, this technique appears adequate and promising to pinpoint the provenance of oyster shells. It remains necessary, however, to have the elemental fingerprint of a multitude of coastal sites as a database, so that the fingerprint of unknown specimens can be recognized. respectively. The authors would like to acknowledge three anonymous reviewers for their proposed improvements of the manuscript and Karin Verrecchia for editing the English.

Data availability

The data used in this study have been deposited in the Zenodo data repository (http://doi.org/10.5281/zenodo.4681223; [START_REF] Mouchi | Data for the publication by Mouchi et al[END_REF].

Declarations of interest: none.

  cathodoluminescence signal) at Baie des Veys on 38 measurements from C. gigas specimens (Cg_BdV_M; n=5) and 38 measurements from O. edulis specimens (Oe_BdV_M; n=5), and at Marennes-Oléron on 23 measurements from C. gigas specimens (Cg_MO_M; n=3) and 16 measurements from O. edulis specimens (Oe_MO_M; n=3). Two methods were tested: unweighted average distance and "Ward" inner squared distance. The cophenetic correlation coefficient between the distances obtained from the cluster tree and the original distances (in the multivariate space) is an indicator of the quality of the estimation (of the hierarchical cluster analysis) to faithfully represent dissimilarities between observations being compared. Comparison and classification between all groups were performed using t-SNE (t-distributed Stochastic Neighbour Embedding; van der Maaten

Figure 2 :

 2 Figure 2: Distribution of the 472 measurements for each analyte. a: Boxplot with values indicated in cubic root for best approaching the normality in the distribution and for clarity. Note that the ordinate axis is in logscale. b: Coefficient of variation for each analyte.

Figure 3 :

 3 Figure 3: Visualization of shell groups' partitioning using t-SNE (exact method using Euclidean distance) based on 472 measurements. Crassostrea gigas shell measurements are indicated with dots. Light blue symbols represent measurements from shells originating from Normandy, dark blue symbols represent Brittany, red symbols represent Charente, purple symbols represent Arcachon Basin, and green and orange symbols represent the Mediterranean Sea, including Corsica. Two groups of Antiquity shells of unknown origin are indicated by yellow and black circles.

Figure 4 :

 4 Figure 4: Visualization of O. edulis shell groups' partitioning using t-SNE showing specific groups. a: t-SNE from Figure 3 with measurements from C. gigas shells shaded. b: Same as a showing only groups of confirmed English Channel or Atlantic Ocean coastlines. The dashed ellipse indicates a cluster common to measurements from several groups. c: Same as a showing only the three Charente groups. d: Same as a showing only the groups from the Mediterranean coastline and the two groups of unknown origin.

Figure 5 :

 5 Figure 5: Successive measurements along shell growth indicating different fingerprints on specimen 1297-2 from Broue tower (Oe_BrT_A), namely "common cluster" and "secondary cluster" on Fig. 4b. a: Cathodoluminescence view of the last two years of growth (corresponding to the area analysed by LA-ICP-MS) on the umbo region and positions of the measurements. Red and blue dots correspond to summer and winter measurements, respectively. Note that the first measurement on the left is out of the picture frame. b: Intensity of luminescence (in arbitrary units) along the umbo following the line indicated on a, showing the areas of light (corresponding to summers; "S") and dull (corresponding to winters; "W") luminescence. Smoothing of this signal is shown overlapping in red, and highlights the presence of a monthly tidal cycle (numbers).

;

  [START_REF] Lazure | Manche -Mer du Nord -Etat physique et chimique, caractéristiques physiques, courantologie. Caractéristiques et Etat Ecologique[END_REF], hence limiting the influence of the composition of the freshwater coming from the land.[START_REF] Becker | Spatial and temporal variation in trace elemental fingerprints of mytilid mussel shells: A precursor to invertebrate larval tracking[END_REF] measured the elemental fingerprint of young mussel shells in several localities in San Diego County to study the larval dispersal along the coastline. They were able to observe distinct fingerprints for specimens from estuarine locations, but could not clearly discriminate open marine locations. We suspect this is what is observed here, with a 'marine' fingerprint common to several groups (being identical to that of an open marine locality, i.e. Baie des Veys), and a secondary fingerprint, specific to the locality, in environments where continental inputs are sufficiently strong.

(

  2005) were able to retrieve an open marine origin from San Diego County mussel shell composition with a 5% error (out of 41 measurements), which corresponded to misclassification to open bays. However, although these specimens had been collected from six open marine localities along the Pacific Ocean coastline of San Diego, the correct identification of the site had a 44 % error. Misclassification of specimens from two open bays (between both bays and all open marine localities as a whole) presented errors of 33 and 62 % (out of 6 and 8 measurements, respectively) due to a misclassification of measurements to open marine environments. These results indicated that the open marine fingerprint from the mussel shell composition is similar along 50 km of coastline (range of the collection sites of the study), and that elemental fingerprinting of specimens from intertidal localities

  carried out by the support of the TracOstrea project, funded by the Collaborative Research Project "Les Marais Charentais au Moyen-Âge et à l'époque moderne" and the University of Rennes 1. The authors would like to thank Frédéric Delbès and Mikaël Guiavarch for the work they performed on the preparation of the thin sections. We thank Michel Ropert from the marine station of Port-en-Bessin (IFREMER), Philippe Geairon from the marine station of La Tremblade (IFREMER) and Danièle Maurer from the marine station of Arcachon (IFREMER) for their precious help in the breeding of oysters at Baie des Veys, Marennes-Oléron and Tès localities,

Table 1 :

 1 Number of specimens and measurements per group of shells. Sample collection sites. Modern specimen groups are indicated by circles while archaeological specimen groups are presented as stars. Colour coding corresponds to that used in the t-SNE analysis on the following figures. 1: Baie des Veys modern site. Arrow: direction of the English Channel current near the shore[START_REF] Lazure | Manche -Mer du Nord -Etat physique et chimique, caractéristiques physiques, courantologie. Caractéristiques et Etat Ecologique[END_REF]. 2:

	111
	Group	Location	Coordinates	ID on	Age	Number of	Number of
				Figure		specimens	measurements
				1
	Cg_BdV_M	Baie des Veys	49°23.11 N,	1	Modern	5
		(Normandy)	01°06.05 W
	Oe_BdV_M	Baie des Veys	49°23.11 N,	1	Modern	5
		(Normandy)	01°06.05 W
	Oe_Lan_A	Landévennec (Brittany)	48°17.42 N,	2	8 th -15 th c. CE	2
			04°16.00 W
	Oe_RsB_M	Riec-sur-Bélon (Brittany)	47°48.55 N,	3	Modern	3
			03°43.08 W
	Cg_MO_M	Marennes-Oléron	45°52.23 N,	4A	Modern	3
		(Charente-Maritime)	01°10.60 W
	Oe_MO_M	Marennes-Oléron	45°52.23 N,	4A	Modern	3
		(Charente-Maritime)	01°10.60 W
	Oe_Bro_A	Brouage (Charente-	55°85.05 N,	4B	16 th c. CE	4
		Maritime)	01°07.44 W
	Oe_BrT_A	Broue tower (Charente-	45°47.67 N,	4C	11 th -14 th c.	6
		Maritime)	00°58.57 W		CE
	Cg_Tes_M	Tès (Arcachon Basin)	44°40.01 N,	5	Modern	8
			01°08.18 W
	Cg_Leu_M	Leucate (Aude)	42°52.48 N,	6	Modern	5
			03°01.50 E
	Cg_Dia_M	Diana lagoon (Corsica	42°08.00 N,	7A	Modern	2
		Island)	09°32.17 E
	Oe_Ale_A	Aléria (Corsica Island)	42°06.50 N,	7B	Antiquity	2
			09°30.72 E
	Oe_X1_A	Unknown	Unknown	-	1 st c. CE	6
	Oe_X2_A	Unknown	Unknown	-	1 st c. CE	8
	Oe_Mal_A	Undefined (Mediterranean	Unknown	-	6 th c. CE	6
		Sea)

  is probably high enough to induce a strong influence in terms of chemical composition of the local seawater. That way, even at high tide, the composition of shells, derived from the local seawater, would still be substantially influenced from the Aulne River. The two sub-groups from this cluster, with each corresponding to one specimen of different age (8 th and 15 th centuries CE), are probably due to a different degree of the influence of freshwater input. As several centuries separate these specimens, it is probable that the local geomorphology and water level changed, causing relocation of the supply of oysters from the inhabitants. Nevertheless, the close similar fingerprint of these specimens implies that they originate from nearby locations.For C. gigas groups, all measurements from the Cg_Leu_M group of shells are clustered away from the other groups in Figure3, with highly negative values on the 1 st dimension axis and highly positive values on the 2 nd dimension axis. This observation indicates that the elemental compositions of these shells are substantially different than the shells from the other regions. The fingerprint from the Diana lagoon (Cg_Dia_M) is also easily discriminated from the others. Compared to the other regions, these localities are characterized by small watersheds, uniquely consisting in limestone. As one can expect a smaller surface exposed to weathering than those of the other studied localities, the abundance of some elements in the shell should be lower at Leucate and Diana lagoon. Indeed, both groups are characterized by low abundances of light REE compared to the other groups, but the other

	analysed elements do not present such a systematic lower content (Supplementary Information 4). It
	is also possible that, as both these localities are in restricted lagoons, salinity can be at play. Indeed,
	Diana Lagoon (Mediterranean Sea) is known to have a high salinity (37.3 annual mean; Orsini et al.,
	2001) and Leucate lagoon has a seasonally dominated salinity (26 in winter and 42 in summer due to
	evaporation; Andrisoa et al., 2019), whereas Baie des Veys and Marennes-Oléron (Atlantic coastline)
	have lower salinities (33.4 and 32.4 annual mean, respectively; Lartaud et al., 2010c).


R Abgrall

A Larat

M Ricchiuto

C Tavé

A simple construction of very high order non oscillatory compact schemes on unstructured meshes

 have been constructed very high order residual distribution schemes for scalar problems. They were using triangle unstructured meshes. However, the construction was quite involved and was not very flexible. Here, following [1], we develop a systematic way of constructing very high order non oscillatory schemes for such meshes. Applications to scalar and systems problems are given.

Introduction

We are interested in the approximation of the following model problem

λ • ∇u = f x ∈ Ω u = 0 x ∈ Γ - (1)
where Ω ⊂ R d is a polygonal set, Γ -is the inflow boundary

Γ -= {x ∈ ∂Ω, λ • n x < 0}
and n x is the local normal at the point x ∈ ∂Ω.

We consider a conformal triangulation T h which elements K are triangles, quads in 2D or tetrahedrons/hex in 3D. More general elements could, in principle, be considered. The parameter h denotes the maximum value of all the diameters h K of the circumscribed circle/sphere to the elements K of T h . We also assume that the meshes are regular. The mesh T h are assumed to be adapted to [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF], i.e. Γ -is a collection of edges/faces of T h . The space V p h is the set of continuous functions that, on each element K, are polynomials of degree p that vanishes on Γ -.

The equation ( 1) is discretized by a variational formulation of the type : Let us give p, q ∈ N ⋆ , find u h ∈ V p h such that for all v h ∈ V q h , a(u h , v h ) = ℓ(v h ).

(

The examples we are interested in are the SUPG schemes [START_REF] Hughes | A new finite element formulation for CFD: I. symmetric forms of the compressible Euler and Navier-Figure 8: 4-state Riemann problem, comparison stabilized, unstabilized solution on a coars mesh Stokes equations and the second law of thermodynamics[END_REF][START_REF] Hughes | A new finite element formulation for CFD: III. the generalized streamline operator for multidimensional advectivediffusive systems[END_REF] and the stabilized residual distribution schemes [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF] which are in general non linear schemes, even for linear problems.

In the SUPG schemes, we take p = q ∈ N ⋆ and the relation [START_REF] Abgrall | Construction of very high order residual schemes for unstructured meshes[END_REF] writes

Ω v h λ • ∇u h -f (x) dx + h Ω λ • ∇v h λ • ∇u h -f (x) dx = 0 (3)
i.e.

a(v h , u h ) = Ω v h λ • ∇u h dx + h Ω λ • ∇v h λ • ∇u h dx ℓ(v h ) = Ω v h f (x)dx + h Ω λ • ∇v h f (x)dx (4)
In the second example, the Residual Distribution schemes (RD schemes for short), we also take q = p ∈ N ⋆ but the formulation is completely different in order to account, for example, of a maximum principle. These schemes are described in section 2.

The solutions of (2) are obtained by an iterative scheme. The convergence of the iterative procedure is important for two reasons 1. The uniqueness of the solutions of ( 2) is essential to have a well posed problem, and one wishes to obtain a good approximation of the solution of (2).

2. One can show, and we recall this later in the text, that if the the problem (2) is not solved with enough precision, the formal accuracy of the scheme (2) is lost.

In this respect, the SUPG scheme is dissipative, and coercive in a proper norm, so that existence and uniqueness is guarantied. However, in the case of the RD schemes, this may be no longer true, at least for their unstabilized version. Indeed, a RD scheme can be constructed so that it is positivity preserving but in general, the solution of (2) may not be unique, see [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF].

Coming back to the SUPG scheme, we see that the forms a and b are the sum of two terms, the forms

a ′ (u h , v h ) = Ω v h λ • ∇u h dx and ℓ ′ (v h ) = Ω v h f dx (5)
that define the Galerkin formulation of (1). The problem (2) with a ′ and ℓ ′ is known to be very unstable. It is stabilized by adding a dissipative term q and a linear form b to keep the consistancy of the scheme,

q(u h , v h ) = h K D K (v h , u h ), D K (v h , u h ) = K λ • ∇v h λ • ∇u h dx b(v h ) = h Ω λ • ∇v h f (x)dx (6)
The exact evaluation of D K may be quite costly in practice. If in the case of p = 1, the terms λ • ∇u h λ • ∇v h can be evaluated with second order accuracy with only one point (the centroid). For p = 2, the components of ∇v h are of degree one, and an exact quadrature formula (for a constant velocity) λ is obtained with 3 quadrature points (the mid-points of the edges of the triangle) of a 5 point formula as indicated in Table 1.

When p = 3, λ • ∇u h λ • ∇v h
is of degree 4 and 7 quadrature points are needed. This can be seen from Table 1 where weights and quadrature points are displayed for triangular elements. For Q k elements, the situation is worse. The question we are interested in this paper is the following. Given a scheme of the type (2),

Error O(h k ) Weights Λ 1 Λ 2 Λ 3 3 1/3 1/2 1/2 0 1/3 0 1/2 1/2 1/3 1/
a(u h , v h ) = ℓ(v h )
what are the requirements about the forms q and b such that the scheme

a(u h , v h ) + q(u h , v h ) = ℓ(v h ) + b(v h ) (7)
is well posed, has provable error estimates in a well behaved norm ? How can q and b be chosen such that the evaluation of these terms is as simple as possible with the minimal number of operations ?

The schemes we are interested in, like (3) or the RD scheme, share several formal properties in common. Namely, 1. if u is a smooth solution of (1), then for any v h ∈ V q h , we have

a(u, v h ) = ℓ(v h ) (8a)
Moreover, if u h denotes the solution of scheme, we have

a(u -u h , v h ) = 0. ( 8b
)
Note that this property, which is well known for the SUPG scheme, is also true for the RDS scheme, even if the RDS scheme is non linear.

2. From this, if u h denotes now the interpolant of the exact solution u of (1), then the equivalent equation of the scheme is

a(u h , v h ) -ℓ(v h ) = O(h p+d ) (9)
from which we deduce the formal order of accuracy. Of course, in the case of the SUPG scheme, things can be made more rigorous.

These properties must remain intact.

In the first section, we explain in detail what is a RD scheme. The SUPG schemes are particular cases. The second section is devoted to the describe and discuss natural necessary conditions. The third section is devoted to examples and numerical results.

Examples of "unstabilised" schemes

The example of the Galerkin formulation of ( 1) is well known so we skip it. We give some details on the RD schemes that are less known.

We consider a conformal mesh, the generic element is denoted by K. The degrees of freedom are denoted by x σ . In the case of a P 1 interpolant, they are just the vertices of the mesh. For a P 2 interpolant, we have to add the midedge points, etc. Obvious generalization can be described for other continuous elements such as the P k or Q k elements.

In order to construct a RD scheme for (1), on has first to construct "residuals" Φ K σ such that the two conditions are met :

1. Compact stenceil condition : Φ K σ (u h ) := Φ K σ only depends on the values of u at the degrees of freedom in K,

Conservation condition : Φ

K σ are such that σ∈K Φ K σ = ∂K λ • n u h dx - K f (x)dx := Φ K
This is a conservation constraint.

The function u h has to be solution of for any σ,

K,σ∈K Φ K σ = 0. ( 10
)
As said previously, the SUPG schemes are examples of RD schemes, since they are exactly [START_REF] Lerat | Residual-based compact schemes for multidimensional hyperbolic systems of conservation laws[END_REF] with

Φ K σ = K ϕ σ λ • ∇u h -f dx + h K λ • ∇ϕ σ λ • ∇u h -f dx.
Of course these conditions (conservation and a compact stencil) are not enough to provide a working scheme in term of stability and accuracy. Here we foccus on the L ∞ stability and the residual property (8a)-(8b) which ensure formal accuracy. These two additional constraints are achieved by the following procedures 1 .

One starts from a monotone scheme, say the Lax-Friedrichs one,

Φ K,LxF σ = 1 N K Φ K + α K σ ′ ∈K (u σ -u σ ′ )
which is only first order, N K represents the number of degree of freedom in K.

Then we define

x σ := Φ K σ Φ K , they sum up to unity thanks to the conservation relation and

β K σ := x + σ σ ′ ∈K x + σ ′ . ( 11
)
There is no problem in the definition of β K σ since

σ ′ ∈K x + σ ′ ≥ σ ′ ∈K x σ ′ = 1.
The RD scheme is then defined by [START_REF] Lerat | Residual-based compact schemes for multidimensional hyperbolic systems of conservation laws[END_REF] with

Φ K σ = β K σ Φ K . ( 12
)
The solution of (18)-( 12) is sought for by an iterative method. The simplest one is

u n+1 σ = u n σ -ω σ K,σ∈K Φ K σ , for all σ (13)
with u 0 σ = 0 for example, and one hopes that u σ = lim n→+∞ u n σ . Thanks to the definition of β K σ , one can see that the sequence {u n σ } n,σ satisfy a maximum principle provided a CFL-like condition

0 ≤ ω σ max K,σ∈K |K| max σ ′ ∈K max x∈K ||∇ϕ σ ′ (x)|| ≤ 1.
Note that sharper estimates can be given, but this is not the point here.

The variational formulation of ( 10)-( 12) is easily obtained. If one multiply (10) by v(x σ ) and sum over all the degrees of freedom, one obtains

0 = K Ψ K with Ψ K = σ∈K v σ β T σ K λ • ∇u h -f dx = K v h λ • ∇u h -f dx + σ∈K v σ K β K σ -ϕ σ λ • ∇u h -f dx
Since by definition

σ∈T β T σ = 1 = σ∈T ϕ σ , we have σ∈K v σ K β K σ -ϕ σ λ • ∇u h -f dx = 1 N K -1)! σ,σ ′ (v σ -v σ ′ ) K γ σ,σ ′ -ψ σ,σ ′ λ • ∇u h -f dx = h N K -1 ! σ,σ ′ θ σσ ′ v σ -v σ ′ || x σ x σ ′ || K γ σ,σ ′ -ψ σ,σ ′ λ • ∇u h -f dx with θ σσ ′ = || x σ x σ ′ ||/h which is bounded since the mesh is regular, γ σ,σ ′ = β T σ -β T σ ′ and ψ σ,σ ′ = ϕ σ -ϕ σ ′ . The form a is a(u h , v h ) = K v h λ • ∇u h dx + h N K -1 ! σ,σ ′ θ σσ ′ v σ -v σ ′ || x σ x σ ′ || K γ σ,σ ′ -ψ σ,σ ′ λ • ∇u h dx (14) and ℓ is ℓ(, v h ) = K v h f dx + h N K -1 ! σ,σ ′ θ σσ ′ v σ -v σ ′ || x σ x σ ′ || K γ σ,σ ′ -ψ σ,σ ′ f dx (15)
which have the same structure as [START_REF] Corre | A residual-based compact scheme for the unsteady compressible Navier-Stokes equations[END_REF]. The problem of this scheme is that even though the iteration (13) is L ∞ stable, it does not converge in general. The same conclusion holds for more involved iterative scheme and the reason is that (10)-( 12) is not well posed except for very special situations.

An example is given for a second order (hence P 1 interpolation) using the local Lax Friedrichs scheme (12) on

-y ∂u ∂x + x ∂u ∂y = 0 (x, y) ∈ [0, 1] 2 u(x, 0) =    -sin π x -0.7 0.6 if x ∈ [0.1, 0.7] 0 else (16)
The convence history and a solution is given on Figure 1 The solution of

Solution

Convergence history Figure 1: Solution of (1) with the scheme (12) with β defined by [START_REF] Struijs | Fluctuation Splitting Schemes for the 2D Euler equations[END_REF].

is obviously not a second order accurate approximation of (16). The next section is devoted to describe a simple modification of the scheme. This problem has already been solved in [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF] for second order schemes, we show how to extend the method in a simple and efficient way.

Construction and discussion

Let us consider the problem

a(u h , v h ) = ℓ(v h ) (17a)
with a and ℓ given by

a(u h , v h ) = K∈T h a K (u h , v h ) a K (u h , v h ) = K v h λ • ∇u h dx + hb K (u h , v h ) ℓ(v h ) = K∈T h ℓ K (v h ) ℓ K (v h ) = K v h f dx + hl K (v h ) (17b)
We assume that a, ℓ, a K and ℓ K satisfy the following assumptions :

Assumption 3.1. 1. a K and l K are linear in v h .
2. if u is the solution of (1),

a(u, v h ) = ℓ(v h ) for any v h ∈ V h and a(u -u h , v h ) = 0
for any v h ∈ V h . More precisely, because of the structure of the forms a and b, we assume that for any K, and any

v h ∈ V p h , a K (u, v h ) = l K (v h ) and h a K (u -u h , v h ) = O(h p+d ).
These assumtions are true for the SUPG and RD schemes. Moreover, for these two schemes, we have the conservation constraint a K (u h , 1) = l K (1) = 0 for any K.

Remark 3.1 (About the linearity assumption). The problem (1) is linear.

All what is said here can be extended to the non linear case, and the linearity assumption still holds.

The scheme writes in the RD form (14). To see this, we consider the list of degrees of freedom {x σ }. For piecewise linear interpolant and triangular elements or Q 1 interpolant, they are just the vertices of the mesh. For quadratic interpolant and triangular meshes, they are the vertices of the mesh and the mid-edges points, etc. The Lagrange interpolant of degree p associated to a given degree of freedom x σ is denoted as ϕ p σ . We have

1. ϕ p σ (x σ ′ ) = δ σ ′ σ , 2. ϕ p
σ is continuous, 3. for any element T , the restriction of ϕ σ to T is a polynomial of degree d.

By definition, any u h ∈ V p h can be written as

u h = σ u h (x σ )ϕ p σ ,
and v h ∈ V q h can be written as

v h = σ v h (x σ )ϕ q σ ,
so that the scheme can be reformulated as finding u h ∈ V p h such that for any σ, we have

T ∈T h K ϕ p σ λ • ∇u h dx + h K a K (ϕ σ , u h ) + T ⊂T h T f (x)ϕ q σ (x)dx + h K l K (ϕ q ) = 0 (18)
Such a scheme is rewritten as a RD scheme with the residual

Φ T σ = K ϕ p σ λ • ∇u h -f dx + h K a K (ϕ σ , u h ) -h K l K (ϕ σ ) (19)
the conservation constraint is automatically satisfied because σ∈T ϕ p σ = 1 and using the linearity with respect to v h , σ∈K

Φ K σ = K λ • ∇u h -f dx + h K a K u h , σ∈K ϕ σ -l K σ∈K ϕ σ = K λ • ∇u h -f dx = Φ K
Formally, the scheme is accurate because if that if w h ∈ V p h is the interpolant of the exact solution of (1) assumed to be smooth, then Φ T σ (w h ) = O(h p+d ). In fact, for RD scheme, this property is true because Φ

T (w h ) = O(h p+d ) since Φ T (w h ) = ∂T λ • nw h dx - K f dx = ∂T λ • n w h (x) -u(x) dx = O(h d-1 ) × O(h p+1 ) = O(h p+d ).
Using this remark, the relation ( 9) follows for regular meshes.

A first scheme of the type ( 7) can be obtained by perturbing the residual (19) as

Φ T σ ⋆ = Φ T σ + ψ T σ (20)
with the constraint σ∈T ψ T σ = 0 to ensure conservation. The formal accuracy property is also conserved if

ψ T σ (w h ) -f = O(h p+d ) (21) whenever λ • ∇w h = O(h p ).
A first example is obviously given by

ψ T σ = θh K λ • ∇ϕ q σ λ • ∇u h -f dx ( 22
)
where θ is chosen such that

Ω ξ • ∇u λ • ∇u dx + θ Ω λ • ∇u 2 dx ≥ 0. ( 23
)
Under this condition, the iterative scheme ( 13) is convergent when n → +∞. In ( 22), ψ T σ is evaluated by a quadrature formula of exact order,

K λ•∇ϕ q σ λ•∇u h dx = x quad ω quad λ•∇ϕ q σ (x quad ) λ•∇u h (x quad )-f (x quad ))
(24) The question is now : given a formula of the type

ψ T σ = θ T h K |T | x quad ω quad λ•∇ϕ q σ (x quad ) λ•∇u h (x quad -f (x quad )) , ( 25
)
what are the requirements on the points x quad and the weights ω quad , so that we still have the inequality (23) and the accuracy condition (21) ?

Accuracy constraint. Assuming that the solution of ( 1) is smooth enough, we have on T

w h -u = O(h p+1 ) and ∇ w h -u = O(h p ),
and for a regular mesh

∇ϕ σ = O(h -1 )
so that for any x quad ,

λ • ∇ϕ q σ (x quad ) λ • ∇w h (x quad -f (x quad )) = O(h p-1 ) so that ψ T σ = h × O(h d ) × O(h p-1 ) = O(h p+d ).
The formal accuracy is automatically guarantied.

Constraints on the weights and the points x quad . In order to have the inequality (23), a necessary condition is that the quadratic form

q K (v h ) := x quad ω quad λ • ∇v h (x quad )
2 must be positive definite whenever the polynomial λ • ∇v h = 0.

A sufficient condition is for all x quad , ω quad > 0

λ • ∇v h (x quad ) = 0 for each x quad , then λ • ∇v h = 0. (26)
Under these conditions, there exist constants C 1,q and C 2,q such that

C 1,q q K (v h ) ≤ h K K λ • ∇v h 2 dx ≤ C 2,q q K (v h ) (27)
because P q is a finite dimensional space, hence

Q(v h ) = K q K (v h ) defines a norm on V h which is equivalent to the norm v h → K λ • ∇v h 2 dx.
We have shown the following result Proposition 3.2. If a and ℓ are defined by (17b), under the assumptions 3.1 and provided that the conditions (26) hold, for each element K, there exists θ K,0 > 0 such that the scheme (17a) for θ K > θ K,0 . is well posed and .

Proof. The scheme writes in variational formulation : find

u h ∈ V h such that for all v h ∈ V h , a ′ (u h , v h ) = ℓ(v h ) with a ′ = a + b, ℓ ′ = ℓ where b(u h , v h ) = K∈T h σ∈K v σ ψ K σ
with ψ K σ defined by (25) and Using the scalar product

u h , v h = K∈T h |K| σ∈K u σ v σ ,
The iterative scheme (13) writes,

u n+1 = u n -ω Au n -F ) with Au h , v h := a(u h , v h ), F, v h = ℓ(v h ).
The scheme is convergent if

||Id -ωA|| < 1.
A necessary condition is that for any v h ,

-2 Av h , v h + ω||Av h || ≤ 0
for some ω > 0. This condition needs

a(v h , v h ) = Av h , v h > 0 (28)
for any v h . Coming back to the problem,

a(v h , v h ) = ∂Ω + λ • nv 2 h dl + K h K a K (v h , v h ) + θ H q K (v h ) ,
hence a necessary condition for haveing (28) is that on any K we have

a K (v h , v h ) + θ H q K (v h ) > 0
From assumption 3.1, we see that

Ker q K = {v h ∈ P K , q K (v h ) = 0} ⊂ Ker a K = {v h ∈ P K , a K (v h ) = 0},
so that which means that, since q K is positive definite, the scheme is contractant provided that

θ H > θ K,0 = min 0, - sup v h ∈P q a(v h , v h ) inf v h ∈P q ,v h ∈KerqK q K (v h ) ∈ R
for each K.

Examples and numerical illustrations 4.1 Accuracy study

We apply the method on a simple linear problem, namely

∂u ∂y = 0 (x, y) ∈ [0, 1] 2 u(x, y) = sin(πx) 2 x = 0 (29)
for which the solution is simply u(x, y) = sin(π x 2 + y 2 ) 2 .

We have run the formaly second order scheme, third order and fourth order schemes with respectively 1, 3, 6 "quadrature points" in (25).

For the third order scheme, the "quadrature" points are simply the vertices of the triangle. For the third order case, we have chosen the vertices and the mid point edges. The weights are 1, 1/3 and 1/6 respectively. The results, in term of accuracy, are independant of choices of the "quadrature" points, provided that (24) defines a positive definite quadrature form. The constant θ in ( 25 2: Accurary results for (29). The third order accurate scheme uses the vertices in (25). The fourth order scheme uses the vertices and the mid-points (25). In each case, θ = 1. order of accuracy is met in each case. If now we repeat the same experiment with a smaller number of "quadrature" points, the accuracy is degraded and the The "fourth" order accurate scheme uses the vertices in (25). In each case, θ = 1.

h = 1/N L

points 1 point

Figure 2: Isolines of the solution of (29) when 1 point or 3 points are used in (25). The baseline scheme is formaly third order. All the degrees of freedom are represented.

The non linear case

The second example is the solution of the problem

1 2 ∂u 2 ∂x + ∂u ∂y = 0 if x ∈ [0, 1] 2 u(x, y) = 1.5x -0.5 when y = 0 or x ∈ {0, 1} (30)
The solution consists in a compression merging into a shock which foot is located at (0.5, 0.75). Several schemes are tested. We only represent the solutions obtained by the formaly third order scheme since the behavior for the fourth order one is the same. The "quadrature" points are again the vertices of the elements with and without the centroid depending on if we take 3 or 4 points. On Figure 3, we represent the isolines of the scheme when θ is set to 0, 1 or

θ = |T | vertices |k i | min 1, |T | vertices |k i | ∂T n x u 2 2 + n y u dl + ǫ . (31)
In (31), if n i x and n i y are the components of the inward normal opposite the vertex i in the triangle, k i = n i x u 2

2 + n i y u and ǫ = 10 -10 . Once again, the same conclusions hold : 3 points are necessary to get accurate results. We compare the solutions depending on which option is chose (3/4 quad points, the choice of θ). To do this, we make cross-section at y = 0.25, i.e. in the fan, and y = 0.75, i.e. in the discontinuity. The results of Figure 4 show that if θ = 0, the oscillations visible in Figure 3 are not a manifestation of an instability, the scheme is overcompressive. When θ = 1 or is chosen as (31), there is no difference in the solution, whatever the ρ u v p This integral is evaluated by

Φ T = ∂T F (V ) • n -(ξ, ν) • n dl + T V (ξ, ν)dξdν.
Again we see the improvement obtained by adding the term (25). The scheme is very robust and non oscillatory, despites the interaction between many waves.

Conclusion

In this paper, we have discussed a simple way to construct simple and accurate very high order residual distribution schemes. A theoretical discussion is provided which is confirmed by numerical experiments on scalar problems and the Euler equations. We have foccused on schemes like Residual Distribution schemes, bu we believe however that the method we present in this paper can be adapted to other type of schemes. Note also that it shares common features with the work of Corre and Lerat, see [START_REF] Lerat | A residual-based compact scheme for the compressible Navier-Stokes equations[END_REF][START_REF] Lerat | Residual-based compact schemes for multidimensional hyperbolic systems of conservation laws[END_REF][START_REF] Corre | A residual-based compact scheme for the unsteady compressible Navier-Stokes equations[END_REF] for example.

Figure 3 :

 3 Figure3: Results obtained for problem (30) with various choices of θ and quadrature points. In the case of 4 quadrature points we have chosen the centroid (weight -27/48), and the points of coordinates (0.6, 0.2, 0.2), (0.2,0.6,0.2), (0.2, 0.2, 0.6) with weights 25/48. In the case of 3 quadrature points, they are simply the vertices of the triangle with the weights 1/3.

Figure 5 :

 5 Figure 5: Supersonic ject, third order solution

Figure 6 :

 6 Figure 6: With and without dissipation, density isolines.

Table 1 :

 1 Examples of quadrature points and weights for triangles.

	2	0	1/2

  ) is set to unity. The results are displayed in table 2. We see that the expected

	h = 1/N	L 2	rate	L ∞ rate
	25	0.50493 10 -2	-	0.30340 10 -1	-
	50	0.14684 10 -2 1.78	0.12726 10 -1 1.25
	75	0.74684 10 -3 1.66	0.82311 10 -2 1.07
	100	0.41019 10 -3 2.08	0.52882 10 -2 1.54
		Second order accurate results
	25	0.32612 10 -4	-	0.15748 10 -3	-
	50	0.48741 10 -5 2.742 0.31276 10 -4 2.33
	75	0.13334 10 -5 3.19	0.11363 10 -4 2.49
	100	0.66019 10 -6 2.44	0.46897 10 -5 3.07
		Third order accurate results
	25	0.20860 10 -5	-	0.12811 10 -4	-
	50	0.17001 10 -6 3.61	0.17880 10 -5 2.84
	75	0.27027 10 -7 4.53	0.26772 10 -6 4.68
	100	0.91462 10 -8 3.76	0.95526 10 -7 3.58
		Fourth order accurate results
	Table

Table 3 :

 3 Accurary results for (29). The "third" order accurate scheme uses the gravity center in (25). In each case, θ = 1. results are only first order accurate or the scheme is only consistant, see tables 3 and 4. This can also be seen visually on Figure2.

		2	rate	L ∞	rate
	25	0.25122 10 -1	-	0.42887	-
	50	0.12935 10 -1 0.9577 0.39237	0.1283
	100	0.83978 10 -2 0.6232 0.43656 -0.1540

Table 4 :

 4 Accurary results for (29).

Note that other RD scheme exist, they do not satisfy a L ∞ stability property. An example is the SUPG scheme, another one is the LDA scheme, see[START_REF] Deconinck | Compact advection schemes on unstructured meshes[END_REF][START_REF] Struijs | Fluctuation Splitting Schemes for the 2D Euler equations[END_REF].

Acknowledgements

The work of A. Larat is supported by the European projetc ADIGMA on the development of innovative solution algorithms for aerodynamic simulations. The work of C. Tavé was supported by CNES.

The first author would like to take this opportunity to thank Professor Alain Lerat for the major impact he had on his scientific life. I was lucky to have him as a professor during my master studies, and he made me discover that computing could be fascinating. After this, there was no possible return to pure mathematics, my original background.

number of quadrature points. In Figure 4 Some difference appears for the cross section at y = 0.75. The choice (31) appears to be a good compromise.

add the addition of the term (25), the scheme is no longer formaly monotonicity preserving, but the Figure 4 indicate that no undershoot nor overshoot are created. The same figure also indicate that the choice (31) is the best compromise between accuracy and stability. The effect of this term is that when the solution is smooth, θ ≃ 1 while θ ≃ 0 in the discontinuity.

Euler equations

The last examples that we show are for the Euler equations. Details about the scheme can be found in [START_REF] Abgrall | Construction of very high order residual schemes for unstructured meshes[END_REF] in particular about the way equation ( 12) is implemented. The method has been implemented only for P 2 element so far again with 3 "quadrature" points. The first example is a supersonic jet with M = 2.4 on the bottom and M = 4.4 on the top. The solution, see Figure 5, is made of a shock wave followed by a contact and a fan. On Figure 6, we show the effect of adding and removing the term (25). We can also see the increase of accuracy. On Figure [START_REF] Hughes | A new finite element formulation for CFD: III. the generalized streamline operator for multidimensional advectivediffusive systems[END_REF], we have run the first order, second order and third order RD schemes with the same number of degrees of freedom, namely the vertices and the mid-points of the mesh. A last example is a 4 state shock tube problem (configuration 12 of [START_REF] Kurganov | Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers[END_REF]). This case is time dependant, but we can compute the solution at time t since the solution is self-similar,

The case has been chosen that the boundary condition can easily be computed analyticaly. The scheme is the same as before, but we modify the definition of 201 × 201