Hélène Laguerre

Grégory Raymond

Patrick Plan

Nadia Améziane

Xavier Bailly

Patrick Le Chevalier
email: patrick.lechevalier@univ-brest.fr

First description of embryonic and larval development, juvenile growth of the black sea-cucumber Holothuria forskali (Echinodermata: Holothuroidea), a new species for aquaculture in the north-eastern Atlantic

Keywords: Holothuria forskali, sea cucumber, aquaculture, spawning, larvae, embryonic development, metamorphosis, juvenile growth

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

Introduction

Holothurians (Echinodermata: Holothuroidea), known as "sea cucumbers" or "beche-de-mer", are popular luxury food in Asia, called "trepang", and represent an economic important market as high value seafood products [START_REF] Battaglene | Survival and growth of cultured juvenile sea cucumbers, Holothuria scabra[END_REF][START_REF] Chen | Present status and prospects of sea cucumber industry in China[END_REF][START_REF] Conand | Present status of world sea cucumber resources and utilization: An international overview[END_REF]Purcell et al., 2012a;[START_REF] Purcell | Sea cucumber fisheries: Global analysis of stocks, management measures and drivers of overfishing[END_REF]. They are also used in traditional medicine for their nutritional virtues [START_REF] Chen | Present status and prospects of sea cucumber industry in China[END_REF]Vergara & Rodriguez, 2016) and in medicine, pharmacology and in cosmetics for their bioactive compounds with anti tumoral, anti-oxidant and anti-inflammatory properties actually explored and valued [START_REF] Janakiram | Sea Cucumbers Metabolites as Potent Anti-Cancer Agents[END_REF][START_REF] Kiew | Jewel of the seabed: sea cucumbers as nutritional and drug candidates[END_REF][START_REF] Kumar | Antifungal Activity in Triterpene Glycosides from the Sea Cucumber Actinopyga lecanora[END_REF][START_REF] Mena-Bueno | Sea cucumbers with an anti-inflammatory effect on endothelial cells and subcutaneous but not on epicardial adipose tissue[END_REF][START_REF] Santos | Sea cucumber Holothuria forskali, a new resource for aquaculture? Reproductive biology and nutraceutical approach[END_REF].

The increasing demand of commercialized species and the uncontrolled exploitation of wild stocks in the last decades have led to a collapse of natural resources, classifying some species in threatened species (Anderson et al., Due to decrease in natural stocks of tropical species, European species of Holothuroidea are targeted by the Asian Market [START_REF] Gonzalez-Wanguemert | The fast development of sea cucumber fisheries in the Mediterranean and NE Atlantic waters: From a new marine resource to its over-exploitation[END_REF][START_REF] Hamel | Population status, fisheries and trade of sea cucumbers in temperate areas of the Northern Hemisphere. In Sea cucumbers: A global review of fisheries and trade[END_REF]; Atlantic and Mediterranean species are now subject to an important fishing, inducing illegal fishing and effects on biometry and genetic structure [START_REF] Dominguez-Godino | Profitability of the new target sea cucumber species from Mediterranean and NE Atlantic: Holothuria polii[END_REF][START_REF] Gonzalez-Wanguemert | Effects of fishery protection on biometry and genetic structure of two target sea cucumber species from the Mediterranean Sea[END_REF][START_REF] Gonzalez-Wanguemert | Sea cucumbers as new marine resource in Europe[END_REF]. Only few species of European species have been described for hatchery and cultivating procedures, in a pilot-scale, H. arguinensis [START_REF] Dominguez-Godino | A new species for sea cucumber ranching and aquaculture: Breeding and rearing of Holothuria arguinensis[END_REF], H. tubulosa (Rakaj et al., 2018), H. mammata, (Dominguez-Godino and[START_REF] Gonzalez-Wanguemert | The fast development of sea cucumber fisheries in the Mediterranean and NE Atlantic waters: From a new marine resource to its over-exploitation[END_REF] and H. polii [START_REF] Rakaj | Artificial reproduction of Holothuria polii: A new candidate for aquaculture[END_REF].

H. forskali (Delle Chiaje, 1823) is a native unexploited species from Brittany Coast (France) and can be found in the North East Atlantic Area and Mediterranean Sea, in rocky sea beds, from 1 to 100 meters depth. This species is difficult to maintain in captivity because of its behavior in case of stress. It ejects -as a defensive strategy -cuvierian tubules and sometimes even several organs that can lead to death, even if H. forskali, like the other Holothurians and Echinoderms, are able to regenerate most of its organs. H. forskali is a deposit feeder as the majority of Holothurians, and ingest sediment containing fresh or detrital organic matter [START_REF] Purcell | Ecological roles of exploited sea cucumbers[END_REF][START_REF] Massin | Observations écologiques sur Holothuria tubulosa, H. poli et H. forskali (Echinodermata-Holothuroidea) et comportement alimentaire de H. tubulosa[END_REF][START_REF] Roberts | Feeding and digestive strategies in deposit-feeding holothurians[END_REF]. This gonochoric species has an annual reproduction cycle with a synchronous spawning for male and female in spring in the Atlantic Area [START_REF] Despalatovic | New data about reproduction of the holothurian Holothuria forskali (Echinodermata) living in geographically different places[END_REF][START_REF] Santos | Sea cucumber Holothuria forskali, a new resource for aquaculture? Reproductive biology and nutraceutical approach[END_REF][START_REF] Tuwo | Reproductive biology of the holothurian Holothuria forskali (Echinodermata)[END_REF][START_REF] Tuwo | La fécondité de trois holothuries tempérées à développement pélagique. 8th Echinoderm Conference[END_REF]). However, there is no exhaustive references on reproduction behavior, embryonic and larval development of H. forskali.

We report herein for the first time the reproductive behavior, embryonic and larval development, and juvenile growth of H. forskali in captivity.

Material and methods

Broodstock collection and maintenance

Adults of H. forskali (n=36, mean weight = 322.78 ± 68.27g) were collected at Brezellec (West Coast of Brittany, France: 48°4'15.15"N and 4°39'45.29"W) by scuba diving at 5 to 12 meters depth on 5th of April 2018 (fig. 1, the geographic map was obtained from the French national website www.geoportail.gouv.fr) and the seawater temperature was 12°C. Broodstocks were transferred directly in a 30 L box to the hatchery of the marine experimental station of Agrocampus Ouest in Beg-Meil (Fouesnant, France). Broodstocks were placed in a 250 L tank and were maintained in flowing 5 µm filtered seawater at 14°C prior to thermal induction of spawning. Animals were carefully handled to avoid expulsion of cuvierian tubules and organs such as gonads.

Macroscopic and microscopic examination of the gonadal tubules

Before induction of spawning, gonads were removed from the animals by dissection to determine the maturity stages: immature (I), resting (II), maturating (III), pre-spawning (IV) and post-spawning (V) [START_REF] Ramofafia | Reproductive biology of the commercial sea cucumber Holothuria fuscogilva in the Solomon Islands[END_REF][START_REF] Ramofafia | Development of three commercial sea cucumbers Holothuria scabra, Holothuria fuscogilva and Actinopyga mauritiana: Larval structure and growth[END_REF][START_REF] Santos | Sea cucumber Holothuria forskali, a new resource for aquaculture? Reproductive biology and nutraceutical approach[END_REF][START_REF] Tuwo | Reproductive biology of the holothurian Holothuria forskali (Echinodermata)[END_REF]. Macroscopic and microscopic observations were performed on gonadal tubules, in order to ensure the maturity stage: primary (T1), secondary (T2), growing and fecund (T3, T4), and spend tubules (T5) [START_REF] Tuwo | Reproductive biology of the holothurian Holothuria forskali (Echinodermata)[END_REF].

Induction of spawning by thermal stimulation and fertilization

Induction of spawning was performed by thermal stimulation on the 6 th of April 2018. Water temperature was quickly raised to 19°C by continuous flow of 1µm filtered seawater, and was maintained at 19°C until spawning. Photos and videos were taken during spawning to describe spawning behavior. Fertilization in the tanks was performed and monitored.

Embryonic and larval breeding

When spawning occurred, fertilized eggs were removed from the tank with 80µm sieve, and washed with seawater. The seawater used for embryonic and larval rearing was sand filtered, then filtered with 1µm filter bag and sterilized by UV system. Eggs were observed and incubated at 18°C (±1°C) for 12 hours in 10 L flat bottom trays in closed-water systems, before larval breeding, at a density of approximately 10 embryos/ml. Embryos were sampled and were observed every hour to describe cell divisions, and measured under a microscope equipped with a micrometric eyepiece.

After incubation, eggs were transferred to hatchery 35 L and 100 L aerating tanks and maintained in closed-water systems at 17°C (±1°C), at a density of 5-10 larvae/ml, with natural photoperiod. Larvae were sampled and observed daily to describe the larval development (size, color, organs formation). Seawater of hatchery tanks was renewed 3 times a week: twice with partial seawater change (50%) and once with full change (100%), using a 120µm mesh siever. Water temperature was daily monitored, and 5 ml was collected to observed and describe larval development. Between the late Gastrula and the early Auricularia stage, larvae were fed daily with a mixture of the microalgae Isochrysis galbana and Chaetoceros muelleri (1:1), at a final concentration from 5.10 3 to 10 4 cell.mL -1 . From the late Auricularia stage until Doliolaria stage, they were daily fed with I. galbana and Skeletonema costatum (1:1). Microalgae were cultivated in Conway Growth Medium, successively in 500 mL, 2 L, 10 L Erlenmeyers and 310 L cylinder, using the method of the batch culture of microalgae. Salinity was maintained at 35 PSU.

Metamorphosis and juvenile breeding

At Pentactula stage, larvae were transferred in a micro-nursery with airlift systems, in 300 mm diameters sievers containing Plexiglas plates, previously striped with a wire brush and colonized by a fouling of periphytic microalgae. Plates of slightly roughened Plexiglas, which were vertically stacked in the water column with each sheet separated by a spacer. Tanks were placed in black during metamorphosis -settlement at ambient temperature, 18°C (±1°C). Microalgae S. costatum and C. muelleri (2:1) were daily added 2 days after settlement. When juveniles were observed in tanks, fish micropellets (Carp fish feed from Alltech Coppens, The Netherlands) were distributed once a week. Approximate measurements were made with a ruler under the binocular magnifier.

Results

Gonadal morphology in mature broodstock

Gonads of mature broostock appeared as a well-developed tree constituted by a multitude of branched tubules, converging towards the gonadal basis and attached to the body wall. The gonadal basis was divided into primary tubules, and then bifurcated in secondary ramifications and then fecund tubules during the spawning period. Gonadal tubules in mature female was dark orange to red and white to pale pink in mature male (fig. 2 A-B). In April and during the spawning period, mature gonads were observed, corresponding to stage IV: mature and pre-spawning gonads [START_REF] Santos | Sea cucumber Holothuria forskali, a new resource for aquaculture? Reproductive biology and nutraceutical approach[END_REF][START_REF] Tuwo | Reproductive biology of the holothurian Holothuria forskali (Echinodermata)[END_REF]. Gonads occupied a large place in the body cavity and tubules were more than 10 centimeters long and 1 to 2 mm in diameter, corresponding to class T4 fecund tubules [START_REF] Tuwo | Reproductive biology of the holothurian Holothuria forskali (Echinodermata)[END_REF]. Male mature tubules were filled with agglomerated spermatozoa, composed by a tail and a head (fig. 2C). Microscopic examination of female tubules shown mature and densely packed ovaries (fig. 2D-E), composed by mature oocytes (165µm in diameter) with their germinal vesicle (50-55µm in diameter) and enclosed with their follicle. The newly released oocytes had lost their follicle coat and were surrounded by a 20 µm-thick translucent jelly layer. The nucleus (5µm in diameter) was well visible in the germinative vesicle (fig. 2F). A,D,E,F) and in male (B, C) of H. forskali. A: Macroscopic aspect of mature gonadal tubule in female (T4 fecund tubules, stage IV pre-spawning). B: Macroscopic aspect of mature gonadal tubule in male (T4 fecund tubules, stage IV pre-spawning). C: active spermatozoa (spz) with a focus on spermatozoa morphology, tail (tl) and head (hd) in the black box. D-E: Microscopic morphology of ovarian tubules and oocyte: mature vitellogenic oocytes (vo) with their follicle (f), and germinal vesicle (vg). F: view on mature oocyte with a well-visible nucleus (nc) in the germinative vesicle (vg), and the translucent jelly layer (tjl) around the oocyte.

Fig. 2. Macro-and microscopic morphology of pre-spawning gonadal tubules in female (

Spawning behavior

Spawning started in the night of the 8 th of April 2018, after a 2 day-thermal stimulation. Prior to spawning, male and female climbed on tank walls, to the upper water level, and then extended and erected the anterior region of their body, ventral face in front of the surface (fig. 3A). More than half of their body sought the standing position during spawning. Male and female ejected their gametes by the gonopore (fig. 3B), localized dorsally above the mouth (arrow, fig. 3C). Male ejected continuous streams of sperms and female streams of orange mature oocytes (fig. 3C, fig. 4A). Male and female performed oscillations to help the dispersal of the gametes. Animals were spawning successively during 3 days, mainly during the night, starting at the end of the day, corresponding to the decrease of photoperiod. Eggs were recovered, washed and placed in incubation.

Embryonic and larval development

The average diameter of the mature eggs released by females was 165.46 ± 4.29 µm (n=15). Because fertilization was not controlled by adding sperm, duration of fertilization can be only estimated at 30 minutes. Embryonic development of H. forskali started with the elevation of the fertilization envelop (fig. 4,table 1) and the expulsion of the polar bodies (fig. 4B). The embryonic development of H. forskali was radial holoblastic. The first cleavage appeared 2 hours after fertilization and generated 2 equal and hemispheric blastomeres (fig. 4C,table 1). The two blastomeres were symmetrically divided in the following cleavages (fig. 4D-F), egg diameters evolved from 180 to 195 µm until the morula stage (table 1). The blastula (mean 193.33 ± 7.99 µm) stage was observed 24 hours after fertilization, exhibiting a central blastocoele and a peripheral ciliary cover able to ensure continuous rotary movements (fig. 4G). The blastula stage was the first mobile stage in H. forskali embryonic development. The early gastrula stage (mean 242.50 ± 8.2 µm) was reached between 36 hours and 2 days after fertilization by invagination and elongation of the blastula wall and evolved to the late gastrula (mean size 296.67 ± 9.43 µm) 3 to 4 days after fertilization (fig. 4H-I, table 1). Table 1. Embryonic and larval development and size (mean ± SD; N > 15; *N = 5; size range and mean in µm) of H. forskali at 17°C. Time is represented in weeks (w), days (d), hours (hr) and minutes (min), time after fertilization (AF). The common larval stages in Holothuroidea species are auricularia, doliolaria and pentactula, before metamorphosis and settlement. The auricularia stage hatched 4 days after fertilization (table 1). The early auricularia was characterized by a functional, well visible mouth, esophagus, cloaca (fig. 5A-B) and a ciliated band used for feeding and locomotion. At this stage, the larvae was able to feed. The mean length and width were respectively 639.87 ± 40.93 µm and 394.08 ± 89.30 µm (table 1). During the auricularia stages, the larvae continuously elongated and their size evolved from 300 in early stages to 1450 µm in late auricularia (table 1). The mid auricularia was obtained 7 days after fertilization (fig. 5C), axohydrocoel was forming, the intestine increased in size and larvae began to accumulate hyaline spheres, to obtained the late auricularia stage (fig. D-E), with a mean length and width of 1060 ± 138.83 µm and 790.29 ± 111.25µm respectively, the maximum length registered was 1417 µm. The late auricularia was reached after 25 to 35 days post-fertilization, accumulated 10 well-formed hyaline spheres and the edges of the larvae were invaginating. After 4 weeks, the late auricularia progressively began to compact, round and shrink down to obtain a mean length of 772.33 ± 48.61 µm at the non-feeding doliolaria stage reached 36 days after fertilization. Dolioriaria larvae performed rotary movements with their 5 transverse ciliary bands, hyaline spheres were aligned (fig. 5G-H). Quickly, the larvae evolved towards the pentactula stage less than a week later (fig. 5I,table 1), the ultimate stage before metamorphosis, and changed in behavior from a planktonic larvae to a benthic larvae. Ciliated bands disappeared the 5 tentacles were well-visible and the larvae moved to the bottom of the tanks to enhance metamorphosis (fig. 5I). During metamorphosis, aerating was decreased, feeding was stopped and larvae were transferred in settlement tanks.

Metamorphosis and juvenile development

Metamorphosis and settlement of H. forskali began 42 days after fertilization. Duration of settlement still remained undetermined. First juveniles visible to the naked eye were observed 2 months after metamorphosis. More than 200 juveniles were found in the settlement tank. The body of 120-day juvenile was elongated, yellowish and translucent colored (fig. 6A-B). The dermis was covered with ossicles of varying shapes and with papillae along the dorsal body wall; ambulacral podia appeared in the ventral face and helped for movement on substrate. In the oral pole, 10 wellramified tentacles were present. The digestive tract was visible through the dermis and exhibit peristaltic movements. Most of the juveniles had a length between 0.5 and 1.5 mm (mean size 8 mm) and presented a large heterogeneity between specimens. The 150 day-juveniles began to darken and were less translucent (fig. 6C-D); ossicles, papillae and podia were well-developed. The average length was comprised between 5 and 30 mm (mean size 12 mm). The 250-day juveniles were small replicas of adults (fig. E-H), the body color were darkening, a part of the juveniles were dark brown (fig. E-F-H); the length was between 0.5 and 4 cm. The 300-day juveniles sized had a large range size comprised between 1 cm and 6 cm (fig. 6I) and presented an important heterogeneity. The growth of juveniles was approximately 0.1mm per day from the settlement to 300 days post-fertilization (fig. 7). This heterogeneity of juvenile growth may result from the feeding system and otherwise due to the difficulty of length measurements because of the retraction capacity of the body and the fragility of the juveniles during manipulations.

Discussion

The purpose of this study on H. forskali was to describe spawning, embryonic and larval development, metamorphosis, juvenile growth, and associated breeding times in hatchery, prior to aquaculture development.

Reproductive behavior and spawning

Broodstocks of H. forskali were collected in West Brittany in April, during the breeding season as previously described in the Atlantic Region and more precisely in Glenan Archipelago [START_REF] Santos | Sea cucumber Holothuria forskali, a new resource for aquaculture? Reproductive biology and nutraceutical approach[END_REF][START_REF] Tuwo | Reproductive biology of the holothurian Holothuria forskali (Echinodermata)[END_REF] and confirmed by gonadal examination of the collected broodstock. The maturation cycle of H. forskali differed according to the site, as described in the Adriatic Sea [START_REF] Despalatovic | New data about reproduction of the holothurian Holothuria forskali (Echinodermata) living in geographically different places[END_REF] or in the Atlantic Ocean [START_REF] Tuwo | Reproductive biology of the holothurian Holothuria forskali (Echinodermata)[END_REF][START_REF] Santos | Sea cucumber Holothuria forskali, a new resource for aquaculture? Reproductive biology and nutraceutical approach[END_REF]; and theses variations were also recorded between different populations of H. arguinensis and H. mammata in the Iberian Peninsula [START_REF] Marquet | Sea cucumbers, Holothuria arguinensis and H. mammata, from the southern Iberian Peninsula: variation in reproductive activity between populations from different habitats[END_REF].

Induction of spawning is an important challenge to optimize spawning and fertilization rates in invertebrate hatchery. In hatchery culture of bivalves, different techniques for spawning induction are used: thermal shocks, thermal induction, mechanical shocks, biochemical induction, and lethal methods such as gonad laceration [START_REF] Bourne | A Manual for Scallop Culture in British Columbia[END_REF][START_REF] Helm | Experiments in the hatchery rearing of Pacific oyster larvae (Crassostrea gigas Thunberg)[END_REF][START_REF] Helm | Hatchery culture of bivalves -a practical manual[END_REF][START_REF] Loosanoff | Rearing of Bivalve Mollusks[END_REF]Utting & Pencer, 1991). In Holothuroidea, these techniques have been used and the response to the type of stimulation could be different according to the species [START_REF] Rakaj | Spawning and rearing of Holothuria tubulosa: A new candidate for aquaculture in the Mediterranean region[END_REF]. Thermal stimulation demonstrated the best results during the spawning season, both for temperate species from North-Eastern and Mediterranean regions such as H. polii, H. tubulosa, H. mammata and H. arguinensis, (Dominguez-Godino et al., 2015;[START_REF] Dominguez-Godino | Breeding and larval development of Holothuria mammata, a new target species for aquaculture[END_REF][START_REF] Rakaj | Spawning and rearing of Holothuria tubulosa: A new candidate for aquaculture in the Mediterranean region[END_REF][START_REF] Rakaj | Artificial reproduction of Holothuria polii: A new candidate for aquaculture[END_REF] and for tropical species of the Indo-Pacific Region such as H. scabra [START_REF] Kumara | Artificial breeding and larval rearing of three tropical sea cucumber species -Holothuria scabra, Pseudocolochirus violaceus and Colochirus quadrangularis -in Sri Lanka[END_REF][START_REF] Mazlan | Spawning induction and larval rearing of the sea cucumber Holothuria scabra in Malaysia[END_REF]. Others techniques have shown very contrasting results such as mechanical shocks by dry stimulation for the tropical species [START_REF] Hu | Spawning, larval development and juvenile growth of the sea cucumber Stichopus horrens[END_REF][START_REF] Huang | Spawning, larval development and juvenile growth of the tropical sea cucumber Holothuria leucospilota[END_REF], thermal shocks [START_REF] Battaglene | Spawning induction of three tropical sea cucumbers, Holothuria scabra, H. fuscogilva and Actinopyga mauritiana[END_REF][START_REF] Morgan | Spawning of the temperate sea cucumber Australostichopus mollis Levin[END_REF][START_REF] Rakaj | Spawning and rearing of Holothuria tubulosa: A new candidate for aquaculture in the Mediterranean region[END_REF] or phytoplankton induction [START_REF] Dominguez-Godino | Breeding and larval development of Holothuria mammata, a new target species for aquaculture[END_REF][START_REF] Hu | Spawning, larval development and juvenile growth of the sea cucumber Stichopus horrens[END_REF][START_REF] Mazlan | Spawning induction and larval rearing of the sea cucumber Holothuria scabra in Malaysia[END_REF]. Otherwise, artificial induction of oocyte maturation by chemical inducers or by maturation inducing biochemical fractions from echinoderms shown high efficiency on Holothuroidea [START_REF] Léonet | A new method to induce oocyte maturation in holothuroids (Echinodermata)[END_REF][START_REF] Mercier | Perivisceral coelomic fluid as a mediator of spawning induction in tropical holothurians[END_REF] and on H. forskali [START_REF] Léonet | A new method to induce oocyte maturation in holothuroids (Echinodermata)[END_REF], but this technique requires the sacrifice of a large number of mature animals.

Spontaneous and synchronous spawning of H. forskali was observed for both males and females in a delay of 2 days after the thermal stimulation. Therefore, the release of the gametes occurred mainly during the night. The amount of produced and spawned gametes was high, which is a common reproductive strategy for marine animals characterized by an external fecundation and a planktonic embryonic development. Reproductive behavior was similar to those observed in other species, thereby broodstock erected their body, below to the surface of the tank, and ejected sperms or eggs by the gonopore.

This preliminary study demonstrates the necessity to acquire data on the maturation cycle and spawning season of H. forskali in its natural habitat. Additional works on spawning induction and fertilization control have to be processed, to optimized hatchery procedures. Minimally invasive methods and non-lethal methods would be preferred to limit the sacrifice of broodstock. Moreover, to enhance a new and sustainable aquaculture, broostock selection and genetic diversity will be an important point.

Embryonic and larval development

For this first complete description of embryonic and larval development of H. forskali, the larval rearing duration was 40 days (table 2, fig. 8). The embryonic development was holoblastic and the first larval stage, auricularia, hatched after 4 days and lasted more than 30 days (table 2). Pentactula and metamorphosis were obtained after 40 days, this stage was longer than the other European species and Indo-pacific species (12-15 days for H. scabra, [START_REF] Ramofafia | Reproductive biology of the commercial sea cucumber Holothuria fuscogilva in the Solomon Islands[END_REF]15-17 days for Apostichopus japonica, [START_REF] Ito | Technical development in seed production of the Japanese sea cucumber, Stichopus japonicus[END_REF][START_REF] Soliman | Spontaneous captive breeding and larval development in the green and red variants of the Japanese sea cucumber Apostichopus japonicus (Selenka, 1867)[END_REF][START_REF] Yang | The sea cucumber Apostichopus japonicus: History, biology and aquaculture[END_REF] (table 2). However, in other cultivated Echinoderms, larval life duration can be variable, as in the edible seaurchin Strongylocentrotus intermedius whose duration was 20 days [START_REF] Lawrence | Large-scale production of sea urchin (Strongylocentrotus intermedius) seed in a hatchery in China[END_REF], in comparison with Tripneustes depressus, whose duration can reached 54 days [START_REF] Sonnenholzner-Varas | Breeding, larval development, and growth of juveniles of the edible sea urchin Tripneustes depressus: A new target species for aquaculture in Ecuador[END_REF].

The breeding temperature may explain these results in part. In invertebrate hatchery, the temperature plays an important role and can decreased breeding time, for the bivalve Crassotrea gigas, the larval life duration decrease from 27-32 days at 18°C, to 22 days at 21°C, and 18 days at 28°C [START_REF] Ben Kheder | Impact of temperature on larval development and evolution of physiological indices in Crassostrea gigas[END_REF]. For comparison with temperate holothurians, breeding temperature is 24°C for H. polii [START_REF] Rakaj | Artificial reproduction of Holothuria polii: A new candidate for aquaculture[END_REF]), H. arguinensis (Dominguez-Godino et al., 2015) and H. tubulosa [START_REF] Rakaj | Spawning and rearing of Holothuria tubulosa: A new candidate for aquaculture in the Mediterranean region[END_REF]) (table 2). Generally, an increase of temperature reduced breeding time [START_REF] Asha | Effect of feed concentrations on larval growth, survival and development of Holothuria (Theelothuria) spinifera Theel[END_REF][START_REF] Asha | Effects of temperature, salinity and pH on larval growth, survival and development of the sea cucumber Holothuria spinifera Theel[END_REF][START_REF] Hamel | Early development, settlement, growth, and spatial distribution of the sea cucumber Cucumaria frondosa (Echinodermata: Holothuroidea)[END_REF], therefore, larval development of temperate species is longer than tropical species [START_REF] Chen | Larval development of the sea cucumber, Actinopyga echinites (Echinodermata: Holothuroidea)[END_REF][START_REF] Ramofafia | Development of three commercial sea cucumbers Holothuria scabra, Holothuria fuscogilva and Actinopyga mauritiana: Larval structure and growth[END_REF].

Another way to reduce larval life duration is to reduce larval density. In this present study, a high larval density has been chosen (5 to 10 larvae/mL). However, in other species, an average density of 1-2 larvae per mL is commonly used [START_REF] Dominguez-Godino | A new species for sea cucumber ranching and aquaculture: Breeding and rearing of Holothuria arguinensis[END_REF][START_REF] Dominguez-Godino | Breeding and larval development of Holothuria mammata, a new target species for aquaculture[END_REF][START_REF] Rakaj | Artificial reproduction of Holothuria polii: A new candidate for aquaculture[END_REF]. In this experiment, the larvae density was high during all the auricularia stage and a high larval density can increase the larval duration, the risk of pathogen development and diseases, and finally mortality [START_REF] Kennedy | Potential drivers of virulence evolution in aquaculture[END_REF]. Moreover, the high density decrease the food availability [START_REF] Battaglene | Culture of the tropical sea cucumbers for the purpose of stock restoration and enhancement[END_REF][START_REF] Ramofafia | Development of three commercial sea cucumbers Holothuria scabra, Holothuria fuscogilva and Actinopyga mauritiana: Larval structure and growth[END_REF].

The food sequencing and the food availability, from the first feeding larval stage, auricularia, to the non-feeding stage, doliolaria, are crucial points which have consequences on duration, shape, and development of larvae [START_REF] Morgan | Spawning of the temperate sea cucumber Australostichopus mollis Levin[END_REF][START_REF] Morgan | The effect of food availability on early growth, development and survival of the sea cucumber Holothuria scabra (Echinodermata: Holothuroidea)[END_REF][START_REF] Knauer | Growth and Survival of Larval Sandfish, Holothuria scabra (Echinodermata: Holothuroidea), Fed Different Microalgae[END_REF][START_REF] Militz | Successful large-scale hatchery culture of sandfish (Holothuria scabra) using micro-algae concentrates as a larval food source[END_REF] and metamorphosis [START_REF] Ito | Induction of larval metamorphosis in the sea cucumber Stichopus japonicus by periphitic diatoms[END_REF][START_REF] Yang | The sea cucumber Apostichopus japonicus: History, biology and aquaculture[END_REF]. The phytoplankton species Chaetoceros muelleri and Isochrysis galbana are commonly used in hatchery and for holothurians breeding [START_REF] Hu | Larval development and juvenile growth of the sea cucumber Stichopus sp. (Curry fish)[END_REF][START_REF] Hu | Spawning, larval development and juvenile growth of the sea cucumber Stichopus horrens[END_REF][START_REF] Mazlan | Spawning induction and larval rearing of the sea cucumber Holothuria scabra in Malaysia[END_REF][START_REF] Morgan | The effect of food availability on early growth, development and survival of the sea cucumber Holothuria scabra (Echinodermata: Holothuroidea)[END_REF][START_REF] Rakaj | Spawning and rearing of Holothuria tubulosa: A new candidate for aquaculture in the Mediterranean region[END_REF]. Feeding larvae with C. muelleri can increase growth, survival and metamorphosis [START_REF] Knauer | Growth and Survival of Larval Sandfish, Holothuria scabra (Echinodermata: Holothuroidea), Fed Different Microalgae[END_REF]. The use of microalgae concentrates can be explored to diversify food [START_REF] Duy | Use of micro-algae concentrates for hatchery culture of sandfish, Holothuria scabra[END_REF][START_REF] Hair | Potential breakthrough in hatchery culture of Sandfish Holothuria scabra by using algal concentrate as food[END_REF].

Additional experiments will be performed to monitor survival, in order to measure zoootechnical efficiency and to optimize protocols.

Metamorphosis and juveniles growth

Metamorphosis represents a critical point in hatchery, with mortality associated. Succeed the metamorphosis stage depends on parameters such as water quality, substrate, settlement material, food, hydrology, or light intensity (Altaramino et al., 2017;[START_REF] Battaglene | Survival and growth of cultured juvenile sea cucumbers, Holothuria scabra[END_REF][START_REF] Ramofafia | Reproductive biology of the commercial sea cucumber Holothuria fuscogilva in the Solomon Islands[END_REF]Gunay et al., 2015;[START_REF] Indriana | Survival and growth of sandfish Holothuria scabra juveniles reared in various tanks and rearing conditions[END_REF][START_REF] Ito | Induction of larval metamorphosis in the sea cucumber Stichopus japonicus by periphitic diatoms[END_REF][START_REF] Mercier | Settlement preferences and early migration of the tropical sea cucumber Holothuria scabra[END_REF] Growth and feeding methods of sea cucumbers are well-known for cultivated species, but remains elusive for this species [START_REF] Navarro | Feeding and movement patterns of the sea cucumber Holothuria sanctori[END_REF]Robert et al., 2000;[START_REF] Uthicke | Sediment bioturbation and impact of feeding activity of Holothuria (Halodeima) Atra and Stichopus Chloronotus, two sediment feeding Holothurians, at Lizard Island, Great Barrier Reef[END_REF]. Metamorphosis and settlement were effective, which represent a technical challenge. For H. forskali, the use of Plexiglas plates covered with a microalgae fouling was successful. This type of support is common in aquaculture and is used for holothurians settlement [START_REF] Ito | Induction of larval metamorphosis in the sea cucumber Stichopus japonicus by periphitic diatoms[END_REF][START_REF] Morgan | The effect of food availability on early growth, development and survival of the sea cucumber Holothuria scabra (Echinodermata: Holothuroidea)[END_REF]). An additional work on fouling will have to be done, on periphitic diatoms, to optimize settlement rate and to homogenize juvenile growth [START_REF] Ito | Induction of larval metamorphosis in the sea cucumber Stichopus japonicus by periphitic diatoms[END_REF][START_REF] Morgan | The effect of food availability on early growth, development and survival of the sea cucumber Holothuria scabra (Echinodermata: Holothuroidea)[END_REF], on substrate, on density and on the moment to transfer the pentactula larvae in the settlement tanks [START_REF] Battaglene | Survival and growth of cultured juvenile sea cucumbers, Holothuria scabra[END_REF]. Because of it efficiency in shellfish hatchery, settlement on fine crushed shell chips could to be explored. Juveniles have been observed only after 60 days for this first rearing of H. forskali. The mean juvenile size of H. forskali was 0.8 cm after 120 days and more than 4 cm after 300 days (fig. 7,table 2). In the first 200 days, the juvenile growth followed a polynomial equation (y = 0,0004x 2 -0,0022x + 1,3443 with R² = 0.9672), described in other species such as Isostichopus fuscus [START_REF] Mercier | Aquaculture of the Galapagos sea cucumber, Isostichopus fuscus[END_REF]. The growth rate of approximately 0.1 mm per day (fig. 7), lower than the tropical cultivated species Stichopus horrens (0,6mm per day, [START_REF] Hu | Spawning, larval development and juvenile growth of the sea cucumber Stichopus horrens[END_REF] and 0.2-0.3 mm per day after 60 days of culture for H. leucospilota [START_REF] Huang | Spawning, larval development and juvenile growth of the tropical sea cucumber Holothuria leucospilota[END_REF]. These growth values were lower than commercialized species but similar to temperate species (table 2).

Conclusion

Taking into account over fishing of sea cucumbers in many places such as Indo-Pacific, Mediterranean regions and in the Overseas Collectivity of Saint-Pierre-et-Miquelon, European species are a new target for the Asian Market. Moreover, sea cucumber play an important ecological role in marine benthic ecosystems by taking part in bioturbation process and nutrient recycling. H. forskali presents an interesting candidate as deposit feeder in integrated multitrophic aquaculture, and thus offering many opportunities to develop a sustainable aquaculture. In this context, the development of aquaculture is in fact a guarantor for the sustainable management of a new fisheries resource.

This first description of H. forskali embryonic, larval development and juvenile growth aimed to enhance a new aquaculture to promote diversification of production and potential export market.

Even if this work was carried out in a lab-scale, these preliminary results are very encouraging for the domestication of this species, and 5 cm long juveniles have been produced after 250 days of breeding. This experiment is promising and further research is needed to optimized zootechnical procedures and then to create a new innovative and profitable aquaculture.

Acknowledgement and funding information

This work was funded by European maritime and fisheries Fund in France (HOLOFARM project, FEAMP -Innovative Aquaculture) and the Brittany region (BioHol'eau project). Hélène Laguerre was supported by a "Quimper-Bretagne Occidentale" grant for PhD thesis. Thanks to Morgane Nedelec and Julien Delaporte-Blanc for zootechnical maintenance, Yannick Fleury for work supervising and Camille Jegou for sea-cucumbers sampling by scuba-diving.

Authors' contributions

HL designed the aquaculture experiments and analysed data. HL wrote the first draft of this manuscript and PLC reviewed it. All authors read and approved the final manuscript. [START_REF] Ramofafia | Reproductive biology of the commercial sea cucumber Holothuria fuscogilva in the Solomon Islands[END_REF][START_REF] Ramofafia | Development of three commercial sea cucumbers Holothuria scabra, Holothuria fuscogilva and Actinopyga mauritiana: Larval structure and growth[END_REF][START_REF] Yan | Studies on embryonic and larval development of sea cucumber (Holothuria scabra)[END_REF][START_REF] Knauer | Growth and Survival of Larval Sandfish, Holothuria scabra (Echinodermata: Holothuroidea), Fed Different Microalgae[END_REF][START_REF] Battaglene | Survival and growth of cultured juvenile sea cucumbers, Holothuria scabra[END_REF][START_REF] Altamirano | Substrate preference for burying and feeding of sandfish Holothuria scabra juveniles[END_REF][START_REF] Kumara | Artificial breeding and larval rearing of three tropical sea cucumber species -Holothuria scabra, Pseudocolochirus violaceus and Colochirus quadrangularis -in Sri Lanka[END_REF][START_REF] Mazlan | Spawning induction and larval rearing of the sea cucumber Holothuria scabra in Malaysia[END_REF] Holothuria fuscogilva 151,7 ± 5,2 40-70 hr 420-730 mortality 25-27 mortality Indo-Pacific -warm waters [START_REF] Ramofafia | Reproductive biology of the commercial sea cucumber Holothuria fuscogilva in the Solomon Islands[END_REF][START_REF] Ramofafia | Development of three commercial sea cucumbers Holothuria scabra, Holothuria fuscogilva and Actinopyga mauritiana: Larval structure and growth[END_REF] Actinopyga mauritiana 109,9 ± 0,99 40-70 hr 300-740 mortality 25-27 mortality Indian Ocean, Red Sea -warm waters [START_REF] Ramofafia | Reproductive biology of the commercial sea cucumber Holothuria fuscogilva in the Solomon Islands[END_REF][START_REF] Ramofafia | Development of three commercial sea cucumbers Holothuria scabra, Holothuria fuscogilva and Actinopyga mauritiana: Larval structure and growth[END_REF] Apostichopus

Fig. 1 .

 1 Fig. 1. Sampling location (zoom with a white box) and habitat of the sea-cucumber H. forskali (black arrows) in West Brittany, France (site located with a red box).

Fig. 3 .

 3 Fig. 3. Posture and spawning behavior of H. forskali. Males and females erected their body (A) and ejected a stream of sperm (sp) or oocytes (oo) (arrows, B, C) by the gonopore. Scale bars = 3 cm.

Fig. 4 .

 4 Fig. 4. Embryonic development of H. forskali at 17°C. A: mature spawned oocyte / unfertilized egg, and a spermatozoon (spz). B: fertilized egg with a fertilization envelop (fe) and polar bodies (pb). C: First cleavage, 2 cells egg formed with 2 blastomeres (b). D: second cleavage, 4 cells egg. E: 8 cells eggs. F: Morula. G: Rotary blastula, blastocoele (blc), ciliated cove (ci), and ectoderme (ect). H: early gastrula with blastopore (blp). I: mid gastrula. Scale bars = 100 µm.

Fig. 5 .

 5 Fig. 5. Larval development of H. forskali at 17°C. A: newly hatched auricularia; intestine (int), buccal cavity (bc), esophagus (eso) and enterocoeles (ent). B: early auricularia; cloaca (clo). C: Mid auricularia; ossicle (oss), axohydrocoel (axo), ciliary band (cb). D-F. Late auricularia; hyaline sphere (hs), mouth (mou), cardiac sphincter (cs). G: Early doliolaria. H: Doliolaria. I: Pentactula stage; tentacles (ten). Scale bars = 200 µm.

Fig. 6 .

 6 Fig. 6. Juveniles of H. forskali. A-B: 120-day juvenile; A: dorsal view; B: lateral view; tentacle (ten), podia (pod), calcareous ring (cr), ossicle (oss), papillae (pap), anus (an). C-D: 160-day juvenile; digestive tract (dt), ramified tentacle (rten). E-H: 250-day juveniles. G: zoom on juvenile dermis. I: 300-day juveniles.

Fig. 7 .

 7 Fig. 7. Average growth of the juveniles H. forskali. Grey bars represented minimum and maximum sizes, the black square the mean size and the blue line the regression line with the equation and R².

Fig. 8 .

 8 Fig. 8.Rearing time from embryonic to the benthic larval stage pentactula in different species of sea-cucumbers. In blue: embryonic stages; in light grey: auricularia stage; in green: doliolaria stage; in dark grey: pentactula (sources in table 2).

Fig. 5 .

 5 Fig. 5. Larval development of H. forskali at 17°C. A: newly hatched auricularia; intestine (int), buccal cavity (bc), esophagus (eso) and enterocoeles (ent). B: early auricularia; cloaca (clo). C: Mid auricularia; ossicle (oss), axohydrocoel (axo), ciliary band (cb). D-F. Late auricularia; hyaline sphere (hs), mouth (mou), cardiac sphincter (cs). G: Early doliolaria. H: Doliolaria. I: Pentactula stage; tentacles (ten). Scale bars = 200 µm.

Fig. 6 .

 6 Fig. 6. Juveniles of H. forskali. A-B: 120-day juvenile; A: dorsal view; B: lateral view; tentacle (ten), podia (pod), calcareous ring (cr), ossicle (oss), papillae (pap), anus (an). C-D: 160-day juvenile; digestive tract (dt), ramified tentacle (rten). E-H: 250-day juveniles. G: zoom on juvenile dermis. I: 300-day juveniles.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Average growth of the juveniles H. forskali. Grey bars represented minimum and maximum sizes, the black square the mean size and the blue line the regression line with equation and R².

Table 2 .

 2 Comparison of embryonic and larval development: rearing time, size and temperature of cultivated seacucumbers in aquaculture.

Table 2 .

 2 Comparison of embryonic and larval development for cultivated sea-cucumbers in aquaculture: rearing time, size and breeding temperature and distribution area.

	1
			Newly-
	Species	Fertilized egg size (µm)	hatched time AF (days auricularia -in	Auricularia (µm) size range	Pentactula AF (days) stage in time	Breeding (°C) temperature	Juvenile size AF (mm) after 90-120 days	Origine	References
			or hours)
	Holothuria forskali	189,07 ± 3,63	4 d	300-1450	42	17-18	8	NE Atlantic, Mediterranean -temperate waters	This study
	Holothuria tubulosa	172,1 ± 1,7	3 d	430-800	26	24	8,7	NE Atlantic, Mediterranean -temperate waters	Rakaj et al., 2018
	Holothuria mammata	167,17 ± 8,29	2 d	350-630	19-20	25-38	mortality	NE & E Atlantic, E mediterranean -temperate waters	Dominguez-Godino & Gonzalez-Wanguemert, 2018
	Holothuria polii	247,46 ± 15,48	3 d	300-750	9-10	24	21,0 ± 8,2	Mediterranean -temperate waters	Rakaj et al., 2019
	Holothuria arguinensis	171,66 ± 5,21	45 hr	500-920	13-14	27-28	-	NE & E Atlantic, E mediterranean -temperate waters	Dominguez-Godino et al., 2015
	Stichopus horrens	120,60 ± 11,53	45 hr	350-1200	19-27	25-27	18-40	Indo-Pacific, Red Sea, Australia -warm waters	Hu et al., 2010; Hu et al., 2013
	Holothuria leucospilota	140-145	36 hr	300-1500	20-23	29-33	40-70	Indo-Pacific -warm waters	Huang et al., 2018
	Holothuria scabra	154,4 ± 2,9	40-70 hr	300-890	12-15	25-27	30	Indo-Pacific -warm waters


Dr Xi Deng
email: deng.xi98@gmail.com

Yuya Shimizu

Bin Xie

Feng Xiao

Constructing higher order discontinuity-capturing schemes with upwind-biased interpolations and boundary variation diminishing algorithm

Keywords: shock capturing, very high order schemes, boundary variation diminishing

Introduction

Designing shock-capturing schemes for high speed compressible flows involving complex flow structures of wide range scales still remains an unresolved issue and attracts a lot of attention in the computational fluid dynamics community. The coexistence of discontinuities and small-scale flow features poses a big challenge to existing numerical schemes. High order and low dissipative schemes are demanded to resolve high-frequency waves and vortices featured in turbulent flows. However, high order schemes may cause spurious numerical oscillations when solving discontinuities such as shock waves, contacts and material interfaces. With the order increased, a scheme may become less robust or too oscillatory to use. To suppress numerical oscillations associated with high order interpolation, a common practice is to introduce certain amount of numerical dissipations by projecting high order interpolation polynomials to lower order or smoother ones, which is also known as limiting projection. An ideal limiting projection is expected to have as small as possible numerical dissipation, while still to effectively suppress numerical oscillation around critical regions. Excessive numerical dissipation introduced in limiting processes tends to undermine the accuracy even though high order interpolation is employed.

Over the decades, a great deal of efforts has been made to construct high order shock capturing schemes in finite volume or finite difference framework. Some representative methodologies to design shock capturing schemes are TVD (Total Variation Diminishing) [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF], ENO (Essentially Non-oscillatory) [START_REF] Harten | Uniformly high-order accurate non-oscillatory schemes[END_REF][START_REF] Harten | Uniformly high order accurate essentially non-oscillatory schemes III[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock capturing schemes[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock capturing schemes,II[END_REF], WENO (Weighted Essentially Non-oscillatory) [START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF][START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF] and MOOD (Multi-dimensional Optimal Order Detection) method [START_REF] Clain | Loubàere, A high-order finite volume method for systems of conservation laws-multi-dimensional optimal order detection (MOOD)[END_REF][START_REF] Diot | Loubàere, Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Diot | The MOOD method in the three-dimensional case: very-high-order finite volume method for hyperbolic systems[END_REF]. Among these schemes, the WENO family which is built on a weighted average of approximations from all candidate stencils is supposed to be one of widely-used methods. Lots of contributions have been made so far to improve the performance of the original WENO scheme. For example, it has been recognized that the classical WENO scheme generates excessive numerical dissipation that tends to smear out contact discontinuities or jumps in variables across material interface. Thus, a series of new smoothness indicators have been proposed in [START_REF] Henrick | Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points[END_REF][START_REF] Borges | An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[END_REF][START_REF] Ha | An improved weighted essentially non-oscillatory scheme with a new smoothness indicator[END_REF][START_REF] Fan | A new smoothness indicator for improving the weighted essentially non-oscillatory scheme[END_REF][START_REF] Acker | An improved WENO-Z scheme[END_REF] where contributions of the less smooth candidate stencils are optimized to reduce numerical dissipation. Moreover, it is observed that the spectral property of the WENO schemes is inferior to that of low-dissipation linear scheme with highorder polynomials. Some recent efforts have been made to improve the spectral property of the WENO schemes. In [START_REF] Van Lith | Embedded WENO: A design strategy to improve existing WENO schemes[END_REF], a new design strategy named embedded WENO is devised to utilise all adjacent smooth substencils to construct a desirable interpolation thus to improve the behavior of the spectral property. In spite of these works, there are still at least three remaining issues in existing high order shock-capturing schemes which need further investigations.

(1) In the vicinity of discontinuities the intrinsic numerical dissipation still persists and pollutes the numerical solution. The resolution of physical features such as contact surfaces, shear waves and material interfaces will evolve from bad to worse. To obtain a sharp discontinuous solution without numerical oscillations is challenging for any limiting processes.

(2) The limiting processes may regard high gradients and small-scale flow structures as discontinuities. Excessive use of limiting processes will undermine the accuracy and make the spectral property of designed schemes inferior to that of low-dissipation linear scheme. Especially, small-scale flow features in turbulence will be smeared out by falsely triggered limiting processes.

(3) Extension schemes to very high order in finite volume or finite difference method may not be efficient and economical. As stated in [START_REF] Balsara | Monotonicity preserving WENO schemes with increasingly high-order of accuracy[END_REF][START_REF] Gerolymos | Very-high-order WENO schemes[END_REF][START_REF] Balsara | An efficient class of WENO schemes with adaptive order[END_REF], the WENO schemes involve more sub-stencils with increasing order.

To evaluate the smoothness indicators is computationally expensive. Moreover, schemes may be less robust as order is increased. For example, there is a possibility of nonexistence of a smooth stencil when using higherorder schemes since the stencils used by the higher-order WENO schemes are quite large. Thus, in [START_REF] Gerolymos | Very-high-order WENO schemes[END_REF][START_REF] Titarev | Finite-volume WENO schemes for three-dimensional conservation laws[END_REF] a recursive-order-reduction procedure is further adopted to suppress numerical oscillations when using higherorder WENO schemes.

We have recently developed a new class of shock-capturing schemes [START_REF] Sun | Boundary Variation Diminishing (BVD) reconstruction: A new approach to improve Godunov schemes[END_REF][START_REF] Deng | High fidelity discontinuity-resolving reconstruction for compressible multiphase flows with moving interfaces[END_REF][START_REF] Deng | New Accurate and Efficient Method for Stiff Detonation Capturing[END_REF][START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF][START_REF] Deng | A fifth-order shock capturing scheme with two-stage boundary variation diminishing algorithm[END_REF] which hybrid polynomial based nonlinear reconstructions and the THINC (Tangent of Hyperbola for INterface Capturing) function as the reconstruction candidates. Realizing that polynomial-based reconstruction may not be suitable for discontinuities, in [START_REF] Sun | Boundary Variation Diminishing (BVD) reconstruction: A new approach to improve Godunov schemes[END_REF] a non-polynomial jump-like THINC function is employed to represent the discontinuity. The switch between WENO and THINC is guided by the BVD (Boundary Variation Diminishing) algorithms, which select reconstruction functions by minimizing the jumps of the reconstructed values at the cell boundaries. The spectral analysis and numerical experiments [START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF] reveal that THINC function with properly chosen steepness can outperform most existing MUSCL schemes for spatial reconstruction. A new scheme, the adaptive THINC-BVD scheme, was then designed under the BVD principle by employing THINC functions with different sharpnesses to solve both smooth and discontinuous solutions [START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF]. Since THINC is a monotonic and bounded function, the adaptive THINC-BVD scheme is able to eliminate numerical oscillation without any limiting projection. For better use of the numerical property of the high order polynomial function, a fifth-order scheme is devised by employing a linear fifth-order upwind scheme and THINC function with the BVD algorithm in our most recent work [START_REF] Deng | A fifth-order shock capturing scheme with two-stage boundary variation diminishing algorithm[END_REF]. The resulting scheme can retrieve linear upwind scheme in smooth region and capture sharp discontinuities without numerical oscillations. Our practice so far shows that the BVD principle can be used as a general and effective paradigm to design new numerical schemes.

In this paper, we propose a new general framework to construct high order shock capture scheme. These schemes, so-called P n T m -BVD (polynomial of n-degree and THINC function of m-level reconstruction based on BVD algorithm), are designed by employing n + 1 order upwind schemes with high order polynomials and THINC (Tangent of Hyperbola for INterface Capturing) functions with adaptive steepness as reconstruction candidates. The final reconstruction function in each cell is determined with multi-stage BVD (Boundary Variation Diminishing) algorithm so as to effectively control numerical oscillation and dissipation. As shown later in the paper, the proposed P n T m -BVD schemes have following desirable properties: 1) they solve sharp discontinuities with effectively suppressed numerical oscillations; 2) they retrieve the underlying linear high-order schemes for smooth solution over all wave numbers, and thus substantially reduce numerical dissipation errors for small-scale flow structures; 3) they can be extended to high order in an efficient way by simply increasing the order of the underlying upwind scheme.

The remainder of this paper is organized as follows. In Section 2, after a brief review of the finite volume method, the details of the new scheme for spatial reconstruction are presented. The spectral property of the new scheme will also be presented. In Section 3, the performance of the new scheme will be examined through benchmark tests in 1D and 2D. Some concluding remarks and prospective are given in Section 4.

Numerical methods

Finite volume method

In this paper, the 1D scalar conservation law in following form is used to introduce the new scheme

∂q ∂t + ∂ f (q) ∂x = 0, (1)
where q(x, t) is the solution function and f (q) is the flux function. We divide the computational domain into N nonoverlapping cell elements,

I i : x ∈ [x i-1/2 , x i+1/2 ], i = 1, 2, . . . , N, with a uniform grid spacing h = ∆x = x i+1/2 -x i-1/2 .
For a standard finite volume method, the volume-integrated average value qi (t) in cell I i is defined as

qi (t) ≈ 1 ∆x ˆxi+1/2 x i-1/2 q(x, t) dx. (2)
The semi-discrete version of Eq. (1) in the finite volume form can be expressed as an ordinary differential equation

(ODE) ∂ q(t) ∂t = - 1 ∆x ( fi+1/2 -fi-1/2 ), (3)
where the numerical fluxes f at cell boundaries can be computed by a Riemann solver

fi+1/2 = f Riemann i+1/2 (q L i+1/2 , q R i+1/2 ). (4)
The remaining main task is how to calculate q L i+1/2 and q R i+1/2 through the reconstruction process.

Reconstruction process

In this subsection, we give the details of how to calculate q L i+1/2 and q R i+1/2 using the BVD principle. The proposed P n T m -BVD schemes are designed by employing linear weight polynomial of n-degree and THINC function of mlevel as the candidate interpolants. The final reconstruction function in each cell is selected from these candidate interpolants with the BVD algorithm. Next, we introduce the candidate interpolants before the description of the BVD algorithm.

Candidate interpolant P n : linear upwind scheme of n-degree polynomial

A finite volume scheme of (n + 1)th order can be constructed from a spatial approximation for the solution in the target cell I i with a polynomial qPn i (x) of degree n. The n+1 unknown coefficients of the polynomial are determined by requiring that qPn i (x) has the same cell average on each cell over an appropriately selected stencil S = {i-n -, . . . , i+n + } with n -+ n + = n, which is expressed as

1 ∆x ˆxj+1/2 x j-1/2 qPn i (x)dx = q j , j = i -n -, i -n -+ 1, . . . , i + n + . (5)
To construct 2r -1 order upwind-biased (UW) finite volume schemes as detailed in [START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF][START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF][START_REF] Balsara | Monotonicity prserving WENO schemes with increasingly high-order of accuracy[END_REF], the stencil is defined with n -= n + = r -1. The unknown coefficients of polynomial of 2r -2 degree can be then calculated from [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock capturing schemes[END_REF]. With the polynomial qPn i (x), high order approximation for reconstructed values at the cell boundaries can be obtained by

q L,Pn i+ 1 2 = qPn i (x i+ 1
2 ) and q R,Pn i-

1 2 = qPn i (x i-1 2 ). ( 6
)
The analysis of [START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF][START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF][START_REF] Balsara | Monotonicity prserving WENO schemes with increasingly high-order of accuracy[END_REF] shows that in smooth region, the approximation with polynomial ( 6) can achieve 2r -1 order accuracy. In this work, we extend the proposed scheme from fifth order (r = 3) to eleventh order (r = 6) by using polynomials of 4th, 6th, 8th and 10th degree as underlying scheme for smooth solution. In order to facilitate the implementation of the proposed method, we give the explicit formulas of q L,Pn i+1/2 and q R,Pn i-1/2 for (n+1)th-order scheme as follows,

• 5th-order scheme

q L,P4 i+1/2 = 1 30 qi-2 - 13 60 qi-1 + 47 60 qi + 9 20 qi+1 - 1 20 qi+2 , q R,P4 i-1/2 = 1 30 qi+2 - 13 60 qi+1 + 47 60 qi + 9 20 qi-1 - 1 20 qi-2 . (7)
• 7th-order scheme

q L,P6 i+1/2 = - 1 140 qi-3 + 5 84 qi-2 - 101 420 qi-1 + 319 420 qi + 107 210 qi+1 - 19 210 qi+2 + 1 105 qi+3 , q R,P6 i-1/2 = - 1 140 qi+3 + 5 84 qi+2 - 101 420 qi+1 + 319 420 qi + 107 210 qi-1 - 19 210 qi-2 + 1 105 qi-3 . (8)
• 9th-order scheme

q L,P8 i+1/2 = 1 630 qi-4 - 41 2520 qi-3 + 199 2520 qi-2 - 641 2520 qi-1 + 1879 2520 qi + 275 504 qi+1 - 61 504 qi+2 + 11 501 qi+3 - 1 504 qi+4 , q R,P8 i-1/2 = 1 630 qi+4 - 41 2520 qi+3 + 199 2520 qi+2 - 641 2520 qi+1 + 1879 2520 qi + 275 504 qi-1 - 61 504 qi-2 + 11 501 qi-3 - 1 504 qi-4 . (9)
• 11th-order scheme q L,P10 i+1/2 = - Another candidate interpolation function in our scheme makes use of the THINC interpolation which is a differentiable and monotone Sigmoid function [START_REF] Xiao | Revisit to the THINC scheme: a simple algebraic VOF algorithm[END_REF][START_REF] Xiao | A simple algebraic interface capturing scheme using hyperbolic tangent function[END_REF]. The piecewise THINC reconstruction function is written as

qT i (x) = qmin + qmax 2 1 + θ tanh β x -x i-1/2 x i+1/2 -x i-1/2 -xi , (11)
where qmin = min( qi-1 , qi+1 ), qmax = max( qi-1 , qi+1 )qmin and θ = sgn( qi+1qi-1 ). The jump thickness is controlled by the parameter β, i.e. a small value of β leads to a smooth profile while a large one leads to a sharp jump-like distribution. The unknown xi , which represents the location of the jump center, is computed from constraint condition qi = 1 ∆x ˆxi+1/2

x i-1/2 qT i (x)dx. More details about the property of THINC function can be referred in our previous work [START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF].

Following [START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF], in this work we use THINC functions with β of m different values to represent different steepness, to realize non-oscillatory and less-dissipative reconstructions adaptively for various flow structures. A THINC reconstruction function qTk i (x) with β k gives the reconstructed values q L,T k i+1/2 and q R,T k i-1/2 , (k = 1, 2, . . . , m). We will use m up to three in present study.

The BVD algorithm

The underlying high order upwind schemes ( 7)-( 10) can achieve the targeted optimal order for smooth region.

However, numerical oscillations will appear for discontinuous solutions, such as the shock waves in high speed compressible flow. As aforementioned, in the conventional high resolution schemes, nonlinear limiting projections are designed to suppress such numerical oscillations in presence of discontinuous solutions. Unfortunately, these limiting processes usually undermine the accuracy and can hardly retrieve the original high-order linear schemes for smooth solutions of relatively small scales. An ideal limiting projection should maintain as much as possible the numerical properties of the scheme that uses the original polynomial of linear weights. In this work, we present a novel reconstruction scheme based on the BVD principle, where the linear high order schemes are directly used for smooth solutions.

In [START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF], a variant BVD algorithm was devised to minimize the total boundary variation (TBV), which implies that the reconstruction function which fits better with the flow field distribution will give a smaller boundary variation value. Thus using monotonic interpolations results in smaller boundary variation values in presence of discontinuous solution, whereas high-order interpolations are preferred to minimize the boundary variations for smooth region. In the P n T m -BVD schemes presented in this paper, reconstruction function is determined from the candidate interpolants with k-stage (k = m) BVD algorithm so as to minimize the TBV of the target cell. We denote the reconstruction function in the target cell I i after the k-th stage BVD as q<k>

T BV <k> i = q L,<k> i-1/2 -q R,<k> i-1/2 + q L,<k> i+1/2 -q R,<k> i+1/2 (12)
and from the THINC function qTk i (x) with a steepness β k as

T BV T k i = q L,T k i-1/2 -q R,T k i-1/2 + q L,T k i+1/2 -q R,T k i+1/2 . (13)
(II-III) Modify the reconstruction function for cells i -1, i and i + 1 according to the following BVD algorithm

q<k> j (x) = qTk j (x), j = i -1, i, i + 1; if T BV T k i < T BV <k> i . (14)
(III) The final BVD stage (k = m):

(III-I) Given the reconstruction functions q<m-1> i (x) from above stage, compute the TBV using the reconstructed cell boundary values from previous stage by

T BV <m-1> i = q L,<m-1> i-1/2 -q R,<m-1> i-1/2 + q L,<m-1> i+1/2 -q R,<m-1> i+1/2 , ( 15
)
and the TBV for THINC function of β m by

T BV T m i = q L,T m i-1/2 -q R,T m i-1/2 + q L,T m i+1/2 -q R,T m i+1/2 . (16)
(III-II) Determine the final reconstruction function for cell I i using the BVD algorithm as

q<m> i (x) =            qTm i ; if T BV T m i < T BV <m-1> i , q<m-1> i ; otherwise . ( 17
)
(III-III) Compute the reconstructed values on the left-side of x i+ 1 2 and the right-side of x i-1 2 respectively by

q L i+ 1 2 = q<m> i (x i+ 1 2 ) and q R i-1 2 = q<m> i (x i-1 2 ). ( 18
)
In this study, we propose and test fifth order scheme with P 4 T 2 -BVD, seventh order scheme with P 6 T 3 -BVD, ninth order scheme with P 8 T 3 -BVD and eleventh order scheme with P 10 T 3 -BVD. According to previous study in [START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF], in all tests of the present study we use Remark 2. This work provides a new framework to construct high order shock capturing schemes. As shown above, extending the scheme from 7th to 11th order is straightforward by simply applying the high order polynomials of the targeted order as given in ( 8)- [START_REF] Diot | Loubàere, Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF]. Our numerical experiments shows that higher-order schemes beyond 11th order can be also designed by adding more levels of THINC function and BVD algorithm in the same spirit.

β 1 = 1.

Spectral property

We study the spectral property of the proposed scheme by using approximate dispersion relation (ADR) analysis described in [START_REF] Pirozzoli | On the spectral properties of shock-capturing schemes[END_REF]. The numerical dissipation of a scheme can be quantified through the imaginary parts of the modified wavenumber while the numerical dispersion can be quantified with real parts. The spectral properties of different order schemes are shown in Fig. 1 in which the spectral properties of proposed P n T m -BVD schemes at different wave numbers are marked by the circles, and the solid lines in the same color represent the spectral properties of the corresponding high order-linear schemes. It can be seen that the proposed schemes have almost the same spectral property as their underlying linear upwind schemes even at high wavenumber band. However, for WENO-type high order schemes, as shown in [START_REF] Pirozzoli | On the spectral properties of shock-capturing schemes[END_REF][START_REF] Fu | A family of high-order targeted ENO schemes for compressible-fluid simulations[END_REF] accuracy is usually undermined at high wavenumber regime although they can recover to their underlying linear scheme at low wavenumber. The reason is that the WENO smoothness indicators tends to mis-interpret high frequency waves as discontinuities. Contrarily, the proposed P n T m -BVD schemes are designed to reduce numerical dissipations, and thus preserve the high-order upwind schemes even in high wavenumber regime, which is hardly realized by any existing high resolution scheme based on limiting projections using nonlinear weighting. It is noteworthy that the ADR analysis method proposed in [START_REF] Pirozzoli | On the spectral properties of shock-capturing schemes[END_REF] assumes that the time integration is exact and ignores the influence of time integration scheme. However, according to the works [START_REF] Bhumkar | A dispersion relation preserving optimized upwind compact difference scheme for high accuracy flow simulations[END_REF][START_REF] Bose | Analysis and design of a new dispersion relation preserving alternate direction bidiagonal compact scheme[END_REF][START_REF] Sengupta | High Accuracy Computing Methods: Fluid Flows and Wave Phenomenon[END_REF][START_REF] Vanharen | Revisiting the spectral analysis for high-order spectral discontinuous methods[END_REF] the time integration scheme may have strong effects on the overall scheme capability and may change the numerical dissipation and dispersion property of a scheme. Thus the methods proposed in [START_REF] Bhumkar | A dispersion relation preserving optimized upwind compact difference scheme for high accuracy flow simulations[END_REF][START_REF] Vanharen | Revisiting the spectral analysis for high-order spectral discontinuous methods[END_REF] should be implemented when considering the influence of time schemes.

Numerical results

In this section, some numerical experiments are performed as a demonstration of the proposed schemes. Linear advection equations and Euler equation systems will be numerically solved. In all examples, the ratio of specific heats is set to γ = 7/5 and CFL is set to 0.4. Time discretization schemes of the same order as the space-discretization are applied. In the present work, linear strong-stability-preserving Runge-Kutta algorithms developed in [START_REF] Gottlieb | Strong stability preserving properties of Runge-Kutta time discretization methods for linear constant coefficient operators[END_REF] are used.

To solve Euler equations, HLLC (Harten, Lax and van Leer with Contact) Riemann solver [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF] is used to calculate the numerical fluxes. The reconstruction process is carried out in terms of characteristic variables. To solve multidimensional problems, we follow [START_REF] Gerolymos | Very-high-order WENO schemes[END_REF] and use simple dimension-wise implementation for 2D Euler equations for the sake of simplicity. As reported in [START_REF] Zhang | On the order of accuracy and numerical performance of two classes of finite volume WENO schemes[END_REF], this simple multi-dimensional implementation is able to give adequate ] should be adopted.

Accuracy test for advection of one-dimensional sine wave

In order to evaluate the convergence rate of the proposed P n T m -BVD schemes, an advection test of smooth profile was conducted on gradually refined grids. The initial smooth distribution was given by

q (x) = sin (2πx) , x ∈ [-1, 1] . (19)
We ran the computation for one period (at t = 2.0) and summarized the numerical errors and the convergence rates for P n T m -BVD schemes in Table 1. In order to compare with the corresponding high order upwind schemes, we also summarized numerical errors calculated by high order upwind schemes in Table 2. As expected, the proposed P n T m -BVD schemes achieved n + 1 order convergence rates when grid elements were gradually refined. Importantly, we observed that the L 1 and L ∞ errors from the P n T m -BVD schemes were exactly same as those calculated by their corresponding high order linear upwind schemes, which was in line with the conclusion from the spectral property analysis in the previous section. These results confirm that the P n T m -BVD schemes are able to recover their underlying high order polynomial interpolations for smooth solutions.

In order to evaluate the computational cost of the proposed schemes, we measure the computational time for the reconstruction process to get interpolated values q L,R with 10 7 cells using one AMD CPU (AMD Ryzen 2970WX 4.20GHz). We also compare the cost with WENOM []. The measurement is summarized in Table 3. It can be seen that the proposed P n T m -BVD schemes indeed cost much more than WENO schemes. However, with the order increased, the cost is almost the same for P n T m -BVD schemes whereas the cost increases significantly for WENO schemes. To increase the order from seventh-to eleventh-, the P n T m -BVD schemes only change the underlying high order interpolations without additional complexity as WENO does.

Accuracy test for advection of a smooth profile containing critical points

We conducted the accuracy test which was more challenging for numerical schemes to distinguish smooth and non-smooth profiles because the initial distribution contains critical points. It has been reported in [START_REF] Henrick | Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points[END_REF] that WENO type schemes do not reach their formal order of accuracy at critical points where the high order derivative does not simultaneously vanish. Following [START_REF] Henrick | Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points[END_REF][START_REF] Gerolymos | Very-high-order WENO schemes[END_REF], the initial condition was given by

q (x) = sin(πx - sin(πx) π ), x ∈ [-1, 1] . ( 20
)
The computation was conducted for ten period (t = 20). We summarize the numerical errors L 1 and L ∞ of the proposed P n T m -BVD in Table 4 and those of the corresponding upwind schemes in Table 5. It can be seen that the P n T m -BVD schemes are capable of achieving their highest possible order of accuracy in the grid refinement tests.

Compared with their corresponding upwind schemes, the P n T m -BVD schemes almost realize the same L 1 and L ∞ errors, which demonstrates that P n T m -BVD schemes restore upwind schemes even at critical points. However, as shown in [START_REF] Henrick | Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points[END_REF][START_REF] Gerolymos | Very-high-order WENO schemes[END_REF], it is generally difficult for WENO type schemes to recover their underlying high order upwind schemes particularly around the critical points.

Advection of complex waves

In order to examine the performance of the proposed scheme in solving profiles of different smoothness, we further simulated the propagation of a complex wave [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF]. The initial profile contains both discontinuities and smooth regions with different smoothness, which is given by

q (x) =                                      1 6 G (x, β, z -δ) + G (x, β, z + δ) + 4G (x, β, z) if |x + 0.7| ≤ 0.1, 1 if |x + 0.3| ≤ 0.1, 1 -|10 (x -0.1) | if |x -0.1| ≤ 0.1, 1 6 [F (x, α, a -δ) + F (x, α, a + δ) + 4F (x, α, a))] if |x -0.5| ≤ 0.1, 0 otherwise, , (21)
where functions F and G are defined as

G (x, β, z) = exp -β (x -z) 2 , F (x, α, a) = max 1 -α 2 (x -a) 2 , 0 , (22)
a = 0.5, z = 0.7, δ = 0.005, α = 10.0, β = log 2 36δ 2 . ( 23
)
The computation was carried out for one period at t = 2.0 with a 200-cell mesh. The results calculated by the P n T m -BVD schemes were presented in Fig. 2. It can be seen that all of schemes are free of visible numerical oscillations and capture sharper discontinuities than conventional high order FVM schemes. With polynomial degree increased, the extreme points of the initial profile are better resolved. It is also noted that P n T m -BVD schemes produce almost same results for the discontinuity which is resolved by only four cells. Because the BVD algorithm can properly choose THINC reconstruction function across discontinuities, increasing order is only effective to the resolution of smooth regions.

Sod's problem

As one of widely used benchmark tests for shock-capturing schemes, the Sod's problem was employed to test the performance of present schemes in capturing the shock front, contact discontinuity, as well as the expansion wave.

The initial distribution on computational domain [0, 1] was specified as [START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF] (ρ 0 , u 0 ,

p 0 ) =           
(1, 0, 1) 0 ≤ x ≤ 0.5 (0.125, 0, 0.1) otherwise .

Two interacting blast waves

Interactive blast waves involving multiple interactions of strong shocks and rarefaction waves has been introduced

in [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF]. The initial distribution is given by

(ρ 0 , u 0 , p 0 ) =                    (1, 0, 1000) , if 0 ≤ x < 0.1,
(1, 0, 0.01) , if 0.1 ≤ x < 0.9,

(1, 0, 100) , if 0.9 ≤ x < 1.

(

) 27
Reflective boundary conditions are imposed at the two ends of computational domain. Two blast waves are generated by the initial jumps and then evolve, associated with violent interactions among different flow structures.

We used 400 mesh cells as used in most literature for this test problem, and depict the numerical density at time t = 0.038 in Fig. 6 against a reference solution obtained with WENO scheme using a very fine mesh. The solution produced by very high order WENO type schemes on 400-cell mesh can be found in the Fig. 14 [START_REF] Gerolymos | Very-high-order WENO schemes[END_REF] where solutions are smeared significantly. On the contrary, the proposed P n T m -BVD schemes show overall better resolution and can capture the left-most contact discontinuity with only three points.

2D Riemann problems

Two-dimensional Riemann problems which have been proposed and extensively studied in [START_REF] Schulz-Rinne | Classification of the Riemann problem for two-dimensional gas dynamics[END_REF][START_REF] Kurganov | Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers[END_REF] are employed to show that the extension of the new scheme to two dimensions can produce accurate, robust and non-oscillatory solutions for compressible flows of complex structures. With different initial configurations made of 2 × 2 constant states, the 2D Riemann problems are usually employed as benchmark tests to examine the behavior of numerical schemes. It has been noted that Kelvin-Helmholtz instability tends to develop along the interfaces between those states. As reported in [START_REF] Dumbser | A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws[END_REF][START_REF] Gerolymos | Very-high-order WENO schemes[END_REF][START_REF] Ha | An improved weighted essentially non-oscillatory scheme with a new smoothness indicator[END_REF][START_REF] Buchmüller | Improved accuracy of high-order WENO finite volume methods on Cartesian grids[END_REF][START_REF] Abedian | A high-order symmetrical weighted hybrid ENO-flux limiter scheme for hyperbolic conservation laws[END_REF][START_REF] Jung | Fine structures for the solutions of the two-dimensional Riemann problems by high-order WENO schemes[END_REF], capturing the discontinuities with adequate sharpness between different states is very crucial for the development of Kelvin-Helmholtz instability and small-scale structures which can only be generated by fine enough grids or very high order schemes.

We used a computational domain of [-0.5, 0.5] × [-0.5, 0.5], and set the initial conditions as follows,

(ρ 0 , u 0 , v 0 , p 0 ) =                              (
A uniform grid made of 800 × 800 cells was employed in all calculations. The numerical results solved by the proposed schemes were presented in Fig. 7. Compared with the reference solution which can be found in [START_REF] Jung | Fine structures for the solutions of the two-dimensional Riemann problems by high-order WENO schemes[END_REF], the new schemes resolve correctly the major flow features. Moreover, the new schemes obtain small-scale structures reconstruction, the scheme captures the extreme points more accurately with less numerical dissipation. From this case, it is obviously that very high order schemes are preferred to simulate problems involving small-scale structures which are common in shock-turbulence problems.

Shock vortex interaction

In this subsection, a shock-vortex interaction test is simulated. This problem is originally proposed in [START_REF] Rault | Shock-vortex interactions at high Mach numbers[END_REF] and reported in [START_REF] Dumbser | A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws[END_REF]. As stated in [START_REF] Dumbser | A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws[END_REF], this test is considered as a more realistic benchmark test for evaluating high order

where r = (xx c ) 2 + (yy c ) 2 . The strength of the vortex is described with M V = u m /c 0 in which c 0 is the sound speed in upstream of the shock. The specific values of these parameters are M V = 0.9, a = 0.075 and b = 0.175.

The results of density gradient field are presented in the Fig. 11. All schemes produce essentially non-oscillatory results. With order increased, the vortex after interacting with the shock begins to split, which agrees well with the numerical simulation in Fig. 9 of [START_REF] Dumbser | Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinear hyperbolic systems[END_REF]. From the observation in [START_REF] Dumbser | Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinear hyperbolic systems[END_REF], the splitting of the vortex can only be resolved with less dissipative scheme. It verifies hence that the high order schemes presented in this paper are able to obtain high resolution results for problems involving shock waves and smooth vortical features.

Concluding remarks

In this study, a new type of very high order reconstruction schemes, so-called P n T m -BVD, is proposed in the finite volume framework. The new scheme hybridizes polynomial of degree n and THINC function with m-level steepness. The effective reconstruction function in each cell is determined by the BVD algorithm which minimizes the numerical dissipation. High order accuracy can be realized in an efficient and straightforward way by directly increasing the degree of polynomial. We extend the new scheme to eleventh order in present study. The spectral analysis and numerical tests show that the P n T m -BVD schemes can retrieve the underlying low-dissipation linear schemes over all wave numbers for smooth solutions. The benchmark tests verify that the proposed schemes can capture sharp discontinuities with numerical oscillations effectively suppressed. Moreover, as the order is increased, small-scale flow features can be resolved with much less numerical dissipation. Thus this work provides an innovative and practical alternative approach for spatial reconstructions of very high order for finite volume method to solve hyperbolic conservative systems that contain discontinuous and smooth solutions of various scales. In future works, we will apply the proposed schemes to complex systems including acoustic waves, turbulence and moving interfaces.

For example, one of our recent work [START_REF] Deng | Implicit large eddy simulation of compressible turbulence flow with PnTm-BVD scheme[END_REF] shows the P n T m -BVD is able to achieve promising results for compressible turbulence flow.

i

  (x).The k-stage BVD algorithm is formulated as follows. (I): Initial stage (k = 0): (I-I) As the first step, use the linear high-order upwind scheme as the base reconstruction scheme and initialize the reconstructed function as q<0> i (x) = qPn(x) i . (II) The intermediate BVD stage (k = 1, . . . , m -1): (II-I) Set q<k> i (x) = q<k-1> i (x) (II-II) Calculate the TBV values for target cell I i from the reconstruction of q<k> i (x)

  1 and β 2 = 1.8 for P 4 T 2 -BVD, and β 1 = 1.2, β 2 = 1.1 and β 3 = 1.8 for P n T 3 -BVD(n = 6, 8, 10) schemes. Remark 1. The multi-stage BVD algorithm enables to reinforce the desired numerical properties at different stages. In the intermediate stage (k = 1, 2, . . . , m -1), the oscillation-free property is realized. The final stage (k = m) is devised to reduce numerical dissipation to capture sharp discontinuities. The parameter β in the intermediate stage is fixed to realize the non-oscillatory property, while in the final stage the parameter β can be chosen from 1.6 to 2.2. A larger value will result in a sharper discontinuity.

Figure 1 :

 1 Figure 1: Approximate dispersion and dissipation properties for different order schemes. Real parts of modified wavenumber are shown in the left panel, while imaginary parts are shown in the right. The circles represent the spectral property of proposed P n T m -BVD schemes at different wave numbers. The solid lines represent the spectral property of corresponding high order upwind schemes in the same color.

3 Figure 2 : 3 Figure 3 : 3 Figure 4 : 3 Figure 5 :

 32333435 Figure 2: Numerical results for advection of complex waves. The numerical solutions at t = 2.0 with 200 mesh cells are presented.

  1.0, -0.6259, 0.1, 1.0) x ≤ 0.0, y ≥ 0.0 (0.8, 0.1, 0.1, 1.0) x < 0.0, y < 0.0 (0.5197, 0.1, 0.1, 0.4) x > 0.0, y > 0.0 (1.0, 0.1, -0.6259, 1.0) x > 0.0, y < 0.0 .

Figure 6 :Figure 7 : 22 (Figure 8 : 2 (Figure 9 :

 6722829 Figure 6: Numerical results of two interacting blast waves problem for density field at t = 0.038 with 400 cells.

Figure 10 :

 10 Figure 10: Reference solution for the 2D shock-entropy wave interaction problem (left panel) and the density distribution along y=0 cross section by different order schemes (right panel).

Figure 11 :

 11 Figure 11: Numerical Schlieren images obtained by the proposed P n T m -BVD schemes of different orders. The computation was conducted with a mesh spacing of h = 1/100.

Table 1 :

 1 Numerical errors and convergence rate for linear advection test. Results are computed by the proposed P n T m -BVD schemes.

	Schemes Mesh L 1 errors	L 1 order L ∞ errors	L ∞ order
		10	2.493 × 10 -1		3.852 × 10 -1
		20	1.174 × 10 -2	4.41	1.815 × 10 -2	4.41
	P 4 T 2	40	3.986 × 10 -4	4.88	6.309 × 10 -4	4.85
		80	1.274 × 10 -5	4.97	2.002 × 10 -5	4.98
		10	8.518 × 10 -2		1.316 × 10 -1
		20	9.673 × 10 -4	6.46	1.495 × 10 -3	6.46
	P 6 T 3	40	8.350 × 10 -6	6.86	1.319 × 10 -5	6.82
		80	6.686 × 10 -8	6.96	1.052 × 10 -7	6.97
		10	2.733 × 10 -2		4.223 × 10 -2
		20	8.216 × 10 -5	8.38	1.269 × 10 -4	8.38
	P 8 T 3	40	1.816 × 10 -7	8.82	2.870 × 10 -7	8.79
		80	3.659 × 10 -10 8.96	5.756 × 10 -10 8.96
		10	8.716 × 10 -3		1.347 × 10 -2
		20	7.132 × 10 -6	10.26	1.102 × 10 -5	10.26
	P 10 T 3	40	4.041 × 10 -9	10.79	6.388 × 10 -9	10.75
		80	2.051 × 10 -12 10.94	3.227 × 10

-12

10.95 accuracy for problems involving shock waves. It should be noted here that in order to realize truly high-order schemes for multi-dimensional Euler equations, algorithms introduced in

[START_REF] Titarev | Finite-volume WENO schemes for three-dimensional conservation laws[END_REF][START_REF] Zhang | On the order of accuracy and numerical performance of two classes of finite volume WENO schemes[END_REF][START_REF] Buchmüller | Improved accuracy of high-order WENO finite volume methods on Cartesian grids[END_REF]

Table 2 :

 2 Same as table 1, but by the high-order linear upwind schemes.

	Schemes	Mesh L 1 errors	L 1 order L ∞ errors	L ∞ order
		10	2.493 × 10 -1		3.852 × 10 -1
		20	1.174 × 10 -2	4.41	1.815 × 10 -2	4.41
	5th order UW	40	3.986 × 10 -4	4.88	6.309 × 10 -4	4.85
		80	1.274 × 10 -5	4.97	2.002 × 10 -5	4.98
		10	8.518 × 10 -2		1.316 × 10 -1
		20	9.673 × 10 -4	6.46	1.495 × 10 -3	6.46
	7th order UW	40	8.350 × 10 -6	6.86	1.319 × 10 -5	6.82
		80	6.686 × 10 -8	6.96	1.052 × 10 -7	6.97
		10	2.733 × 10 -2		4.223 × 10 -2
		20	8.216 × 10 -5	8.38	1.269 × 10 -4	8.38
	9th order UW	40	1.816 × 10 -7	8.82	2.870 × 10 -7	8.79
		80	3.659 × 10 -10 8.96	5.756 × 10 -10 8.96
		10	8.716 × 10 -3		1.347 × 10 -2
		20	7.132 × 10 -6	10.26	1.102 × 10 -5	10.26
	11th order UW	40	4.041 × 10 -9	10.79	6.388 × 10 -9	10.75
		80	2.051 × 10 -12 10.94	3.227 × 10 -12 10.95

Table 3 :

 3 Numerical errors and convergence rate for linear advection test. Results are computed by the proposed P n T m -BVD schemes.

		7th order 9th order 11th order
	P n T m -BVD 5.25 s	5.33 s	5.41 s
	WENOM []	2.48 s	3.25 s	4.12 s

Table 4 :

 4 Numerical errors and convergence rate for advection of a smooth profile containing critical points. Results are calculated by P n T m -BVD schemes.

	Schemes Mesh L 1 errors	L 1 order L ∞ errors	L ∞ order
		10	1.068 × 10 -1		2.263 × 10 -1
		20	1.497 × 10 -2	2.83	3.772 × 10 -2	2.58
	P 4 T 2	40	7.138 × 10 -4	5.13	1.960 × 10 -3	4.27
		80	2.327 × 10 -5	4.94	6.582 × 10 -5	4.90
		160	7.334 × 10 -7	4.99	2.092 × 10 -6	4.98
		10	7.969 × 10 -2		1.306 × 10 -1
		20	2.539 × 10 -3	4.97	6.565 × 10 -3	4.31
	P 6 T 3	40	2.610 × 10 -5	6.60	8.020 × 10 -5	6.36
		80	2.135 × 10 -7	6.93	6.654 × 10 -7	6.91
		160	1.691 × 10 -9	6.98	5.269 × 10 -9	6.98
		10	3.702 × 10 -2		6.049 × 10 -2
		20	5.008 × 10 -4	6.21	1.287 × 10 -3	5.55
	P 8 T 3	40	1.413 × 10 -6	8.47	4.492 × 10 -6	8.16
		80	2.975 × 10 -9	8.89	9.779 × 10 -9	8.84
		160	5.942 × 10 -12 8.97	1.960 × 10 -11 8.96
		10	1.449 × 10 -2		2.837 × 10 -2
		20	1.222 × 10 -4	6.89	3.180 × 10 -4	6.48
	P 10 T 3	40	9.755 × 10 -8	10.29	3.357 × 10 -7	9.89
		80	5.472 × 10 -11 10.80	1.919 × 10 -10 10.77
		160	4.032 × 10 -14 10.40	1.472 × 10 -13 10.35

Table 5 :

 5 Same as table 4, but by high-order linear upwind schemes.

	Schemes	Mesh L 1 errors	L 1 order L ∞ errors	L ∞ order
		10	1.067 × 10 -1		2.265 × 10 -1
		20	1.793 × 10 -2	2.57	3.874 × 10 -2	2.55
	5th order UW	40	7.138 × 10 -4	4.65	1.960 × 10 -3	4.30
		80	2.327 × 10 -5	4.94	6.582 × 10 -5	4.90
		160	7.334 × 10 -7	4.99	2.092 × 10 -6	4.98
		10	7.969 × 10 -2		1.306 × 10 -1
		20	2.539 × 10 -3	4.97	6.565 × 10 -3	4.31
	7th order UW	40	2.610 × 10 -5	6.60	8.020 × 10 -5	6.36
		80	2.135 × 10 -7	6.93	6.654 × 10 -7	6.91
		160	1.691 × 10 -9	6.98	5.269 × 10 -9	6.98
		10	3.702 × 10 -2		6.048 × 10 -2
		20	5.008 × 10 -4	6.21	1.287 × 10 -3	5.55
	9th order UW	40	1.413 × 10 -6	8.47	4.492 × 10 -6	8.16
		80	2.975 × 10 -9	8.89	9.779 × 10 -9	8.84
		160	5.942 × 10 -12 8.97	1.960 × 10 -11 8.96
		10	1.449 × 10 -2		2.837 × 10 -2
		20	1.222 × 10 -4	6.89	3.180 × 10 -4	6.48
	11th order UW	40	9.755 × 10 -8	10.29	3.357 × 10 -7	9.89
		80	5.472 × 10 -11 10.80	1.919 × 10 -10 10.77
		160	4.032 × 10 -14 10.40	1.472 × 10 -13 10.35
	and the coefficients are

Acknowledgment

This work was supported in part by the fund from JSPS (Japan Society for the Promotion of Science) under Grant Nos. 15H03916, 17K18838 and 18H01366. The career of Xi Deng at M2P2 was made possible by Labex MEC (ANR-10-LABX-0092) and the A*MIDEX project (ANR-11-IDEX-0001-02). Especially, Xi Deng would like to express gratitude to Dr. Pierre Boivin who supported his life and research in France.

The computation was carried out on a mesh of 100 uniform cells up to t = 0.25. The numerical results calculated from the proposed scheme were shown in Fig. 3 for density fields. From the results, we observe that P n T m -BVD schemes can solve the contact discontinuity without obvious numerical oscillations. Compared with the results produced by other high order shock-capturing schemes [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF][START_REF] Henrick | Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points[END_REF][START_REF] Borges | An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[END_REF][START_REF] Ha | An improved weighted essentially non-oscillatory scheme with a new smoothness indicator[END_REF][START_REF] Fu | A family of high-order targeted ENO schemes for compressible-fluid simulations[END_REF][START_REF] Van Lith | Embedded WENO: A design strategy to improve existing WENO schemes[END_REF], the results of the proposed schemes are among the best with the jump contact resolved within only two cells. We also observe that P n T m -BVD schemes of different orders produce similar results across the contact where the THINC function is selected by the BVD algorithm. We also mark the cells where THINC function with β = 1.8 is selected in the Fig. 3. It can be seen the BVD algorithm is able to select the jump-like THINC function across discontinuities.

Lax's problem

Shock schemes should be able to deal with the possibility of existing strong discontinuities such as shock waves and detonation waves [START_REF] Teng | On the transition pattern of the oblique detonation structure[END_REF][START_REF] Teng | Numerical study on unstable surfaces of oblique detonations[END_REF]. In this subsection, we solved the Lax problem [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock capturing schemes[END_REF] in this subsection. The initial condition is given by

With the same number of cells as in the previous test case, we got the numerical results at t = 0.16. The density field is plotted presented in Fig. 4. Obviously, the proposed P n T m -BVD schemes obtain accurate solutions without numerical oscillations. Again, the BVD algorithm chooses the THINC function across the discontinuity, which produces similar results among P n T m -BVD schemes of different order. Compared with numerical results using very high order WENO schemes (see the Fig. 8 in [START_REF] Gerolymos | Very-high-order WENO schemes[END_REF]), the present results are one of best.

Shock density wave interaction problem

Here, we simulate the case proposed in [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock capturing schemes,II[END_REF] containing shocks and complex smooth flow features. The initial condition is set as

The numerical solutions at t = 0.18 computed on 200 are shown in Fig. 5, where the reference solution plotted by the solid line is computed by the classical 5th-order WENO scheme with 2000 mesh cells. With the order increased, the peak of the waves is resolved more accurately.

caused by Kelvin-Helmholtz instability, which is a challenging for other high order schemes as reported in [START_REF] Jung | Fine structures for the solutions of the two-dimensional Riemann problems by high-order WENO schemes[END_REF]. It can be also observed that the vortical structures become richer as the order of schemes increases, which agrees with the observations reported in [START_REF] Gerolymos | Very-high-order WENO schemes[END_REF] that the number of vortices depends on the order of the scheme and mesh resolution.

Double Mach Reflection

As a 2D problems involving strong shocks and vortical structures, the double Mach reflection generated by a Mach 10 hypersonic propagating planar shock reflected by a 30 • ramp [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF] is simulated here to evaluate the performance of the new schemes. This test has been adopted widely to assess the ability of a numerical scheme to capture accurately both the strong reflected shocks as well as the small-scale structures in the re-circulation zone. Without exact solution, the numerical dissipation of a numerical scheme can be visually assessed by the richness of vortex structures resulted from Kelvin-Helmholtz instabilities along the slip line in the re-circulation zone. The numerical schemes with large numerical viscosity tend to smear out these structures. In our test, the computational domain was [0, 3.2] × [0, 1].

A right-moving Mach 10 shock was imposed with 60 • angle relative to x-axis. The solution was computed up to time t = 0.2 with a grid resolution of ∆x = ∆y = 1 200 . The numerical results calculated by the proposed schemes of different order are shown in Fig. 8. Despite of the existence of strong shock waves, all P n T m -BVD schemes successfully reproduce essentially non-oscillatory numerical results. The solution quality of the present numerical results is competitive to other advanced high-resolution schemes, such as [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF][START_REF] Fu | A family of high-order targeted ENO schemes for compressible-fluid simulations[END_REF]. Of particular importance, the proposed schemes obtain well developed vortex structures along the slip line, which demonstrates the robust and low-dissipation property of the schemes.

2D shock-entropy wave interaction problem

In this subsection, we simulated the 2D shock entropy wave interaction problem which is usually used as one of simplest models for flow interactions involved in compressible turbulence. Similar to that specified in [START_REF] Acker | An improved WENO-Z scheme[END_REF], the initial condition is given by

(1 + 0.2 sin (14x cos θ + 14y sin θ) , 0, 0, 1) , otherwise.

with θ = π/6 which is the angle of the initial sine wave front with respect to the x axis. The computational domain is [-5, 5] × [-1, 1] with periodic conditions in y direction. The mesh size of h = 1/40 was used for simulation. We conducted the computation until t = 1.8. The numerical results calculated by P n T m -BVD schemes of different orders are presented in Fig. 9. We also calculated the reference solution by the original WENO scheme with a finer mesh in Fig. 10(a). With the order increased, it can be seen that the high-frequency density field is better resolved. We plot the density distribution along y = 0.0 cross section in Fig. 10(b). Similar to the 1D case, with a higher order