Thomas Lecocq
email: thomas.lecocq@umons.ac.be

Maxence Gérard

Kevin Maebe

Nicolas Brasero

Lauren Dehon

Guy Smagghe

Irena Valterová

Thibaut De Meulemeester

Pierre Rasmont

Denis Michez

Chemical reproductive traits of diploid Bombus terrestris males: Consequences on bumblebee conservation

Keywords: Bee decline, Bumblebees, Conservation, Diploid males, Pre-mating recognition, Reproductive traits

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

Introduction

About 90% of the world's plant species are pollinated by animals, and mainly bees [START_REF] Ollerton | How many flowering plants are pollinated by animals?[END_REF]. This makes bees a critically important functional group for ecosystem services in wild environments as well as in crops [START_REF] Free | Insect Pollination of Crops[END_REF][START_REF] Klein | Importance of pollinators in changing landscapes for world crops[END_REF]. Among bees, females of the social genus Bombus (i.e. bumblebees) are among the most important pollinators of the northern temperate zone [START_REF] Williams | A preliminary cladistic investigation of relationships among the bumble bees (Hymenoptera, Apidae)[END_REF]. However, some bumblebee species, as well as other wild bee species, are obviously experiencing a strong decline in their populations [START_REF] Carvalheiro | Species richness declines and biotic homogenisation have slowed down for NW-European pollinators and plants[END_REF][START_REF] Rasmont | Climatic risk and distribution atlas of European bumblebees[END_REF]. Several factors were proposed to explain this decline: landscape modifications/fragmentations, intensive use of agrochemicals, pathogen infections, competition with alien species, climate change, and interactions between all of these factors (review in Winfree 2010; [START_REF] Potts | Developing European conservation and mitigation tools for pollination services: Approaches of the STEP (Status and Trends of European Pollinators) project[END_REF]. Ultimate consequences are a decrease of population size and an increase of population fragmentation that leads to a loss of genetic diversity through inbreeding and genetic drift (Zayed, 2009;[START_REF] Frankham | Introduction to Conservation Genetics[END_REF].

Loss of genetic diversity is particularly problematic in haplo-diploid species such as bees (Zayed & Packer, 2005;Zayed, 2009). Like most Hymenoptera, bees have a single-locus complementary sex determiner (sl-CSD): males are hemizygote while females are heterozygote for this locus [START_REF] Crozier | Evolutionary genetics of the Hymenoptera[END_REF][START_REF] Heimpel | Sex determination in the Hymenoptera[END_REF].

Inbreeding in haplo-diploid species with sl-CSD brings a heavy load through the production of homozygotes which show male phenotype [START_REF] Duchateau | Diploid males in the bumble bee Bombus terrestris. Sex determination, sex alleles and viability[END_REF][START_REF] Duchateau | Sexual biology of haploid and diploid males in the bumble bee Bombus terrestris[END_REF]. Diploid male larvae can be recognized and killed by workers in social bee species like Apis mellifera, [START_REF] Woyke | Drone larvae from fertilized eggs of the honeybee[END_REF] but they can develop into viable adults in bumblebee species (e.g. Bombus terrestris (L.); [START_REF] Duchateau | Diploid males in the bumble bee Bombus terrestris. Sex determination, sex alleles and viability[END_REF]. The production of diploid bumblebee males can be substantial: half of the diploid eggs, which normally develop into females, result in males after brother/sister mating [START_REF] Duchateau | Diploid males in the bumble bee Bombus terrestris. Sex determination, sex alleles and viability[END_REF][START_REF] Duchateau | Sexual biology of haploid and diploid males in the bumble bee Bombus terrestris[END_REF]. These diploid males decrease the fitness of the local populations [START_REF] Whitehorn | Impacts of inbreeding on bumblebee colony fitness under field conditions[END_REF]. First, the production of diploid males in place of workers implies that the number of the workers decreases up to 50% resulting in a small colony size and a low production of new queens [START_REF] Duchateau | Diploid males in the bumble bee Bombus terrestris. Sex determination, sex alleles and viability[END_REF][START_REF] Duchateau | Sexual biology of haploid and diploid males in the bumble bee Bombus terrestris[END_REF]. Second, as diploid males produce diploid sperm [START_REF] Ayabe | Cytological evidence for triploid males and females in the bumblebee, Bombus terrestris[END_REF][START_REF] Darvill | Triploid bumblebees indicate a direct cost of inbreeding in fragmented populations[END_REF], the mating between monandrous virgin queens and diploid males produces 100% of triploid (i.e. sterile) offspring [START_REF] Payne | The frequency of multiple paternity in bumble bee (Bombus) colonies based on microsatellite DNA at the B10 locus[END_REF][START_REF] Takahashi | Diploid male production in a rare and locally distributed bumblebee, Bombus florilegus (Hymenoptera, Apidae)[END_REF].

Therefore, this process leads to a vortex of extinction: populations decreasing based on extrinsic factors (e.g. habitat fragmentation) show higher inbreeding, lower allelic diversity and production of diploid males, this production leading to an additional decrease of the population (Zayed & Packer, 2005). Consequently populations can be driven to quick extinction (Zayed & Packer, 2005) as currently observed in small and fragmented bumblebee populations [START_REF] Darvill | Population structure and inbreeding in a rare and declining bumblebee, Bombus muscorum (Hymenoptera: Apidae)[END_REF][START_REF] Takahashi | Diploid male production in a rare and locally distributed bumblebee, Bombus florilegus (Hymenoptera, Apidae)[END_REF]. This vortex of extinction could be reduced if females are able to recognize and avoid mating with diploid males.

In monandrous species, like most bumblebee species, the female choice requires pre-mating indicators of male genetic equipment [START_REF] Herzner | Brothers smell similar: variation in the sex pheromone of male European Beewolves Philanthus triangulum F. (Hymenoptera: Crabronidae) and its implications for inbreeding avoidance[END_REF]. Therefore, pre-mating mechanisms of avoidance of diploid males could be selected in monandrous bumblebee queens in view of the deleterious consequences of mating with diploid males. First bioassays show that Bombus diploid males have the same mating success as haploid males in experimental conditions [START_REF] Duchateau | Sexual biology of haploid and diploid males in the bumble bee Bombus terrestris[END_REF]. However, such experiments could lead to unnatural mating because queens have no choice between different sexual partners. For example, sibling individuals mate in experimental condition even if B. terrestris queens seem to be less willing to mate with their siblings than with an unrelated partner when they have the choice [START_REF] Whitehorn | Impacts of inbreeding on bumblebee colony fitness under field conditions[END_REF]. Therefore, we expect that premating cues produced by males could allow female pre-mating discrimination to avoid diploid males in natural conditions. As diploid males show differences in body size and physiology [START_REF] Duchateau | Sexual biology of haploid and diploid males in the bumble bee Bombus terrestris[END_REF][START_REF] Gerloff | Effects of inbreeding on immune response and body size in a social insect, Bombus terrestris[END_REF], we hypothesize that they could display original chemical reproductive traits.

The present study investigates the cephalic labial gland secretions (CLGS) of haploid versus diploid males of B.

terrestris. Courtship signals of bumblebees include both behavioral and chemical features (see Baer 2003 for a review). Among the features, the CLGS act as a key trait for the pre-mating recognition [START_REF] Bergström | Complexity of bumble bee marking pheromones: biochemical, ecological and systematical interpretations[END_REF][START_REF] Baer | Bumblebees as model organisms to study male sexual selection in social insects[END_REF][START_REF] Ayasse | Chemical ecology of bumble bees[END_REF]: Most bumblebee males patrol along paths where they scent-mark objects with their CLGS that attract conspecific virgin females (see experimental demonstration that CLGS of bumblebee males attract virgin queens in [START_REF] Coppée | Age-dependent attractivity of males' sexual pheromones in Bombus terrestris (L.) [Hymenoptera, Apidae][END_REF]Lecocq et al., 2015b). This chemical reproductive trait consists of a species-specific complex mixture of (mainly aliphatic) compounds synthesized de novo [START_REF] Calam | Species and sex-specific compounds from the heads of male bumblebees (Bombus spp.)[END_REF][START_REF] Lecocq | Molecular and chemical characters to evaluate species status of two cuckoo bumblebees: Bombus barbutellus and Bombus maxillosus (Hymenoptera, Apidae, Bombini)[END_REF][START_REF] Lecocq | Scent of a break-up: phylogeography and reproductive trait divergences in the redtailed bumblebee (Bombus lapidarius)[END_REF] with a interindividual variation [START_REF] Coppée | Intraspecific variation of the cephalic labial gland secretions in Bombus terrestris (L.) (Hymenoptera: Apidae)[END_REF][START_REF] Žáček | Comparison of age-dependent quantitative changes in the male labial gland secretion of bombus terrestris and Bombus lucorum[END_REF][START_REF] Lecocq | Scent of a break-up: phylogeography and reproductive trait divergences in the redtailed bumblebee (Bombus lapidarius)[END_REF] which leads to differential female preferences for different kinds of males (e.g. young versus old males, sibling versus unrelated males; [START_REF] Coppée | Bombus terrestris (L. 1758) : A complex species or a species complex ? -Intraspecific pheromonal and genetic variations of Bombus terrestris (L.), Impacts on the speciation[END_REF][START_REF] Coppée | Age-dependent attractivity of males' sexual pheromones in Bombus terrestris (L.) [Hymenoptera, Apidae][END_REF]Lecocq et al., 2015c).

Therefore variations in this trait could allow the virgin queen to detect diploids and avoid mating with them.

Material and Methods

Collection and Rearing of Specimens

In order to produce diploid males (diploid male samples were also used in [START_REF] Gérard | Discrimination of haploid and diploid males of Bombus terrestris (Hymenoptera; Apidae) based on wing shape[END_REF], we obtained three colonies of Bombus terrestris terrestris (L.) from bee breeder Biobest N.V. We were reared these colonies at the laboratory of Zoology in Mons (UMons, Belgium) in a controlled climate room (see rearing protocol in [START_REF] Lhomme | Born in an alien nest: how do social parasite male offspring escape from host aggression?[END_REF][START_REF] Lhomme | A method for yearround rearing of cuckoo bumblebees (Hymenoptera: Apoidea: Bombus subgenus Psithyrus )[END_REF][START_REF] Gérard | Discrimination of haploid and diploid males of Bombus terrestris (Hymenoptera; Apidae) based on wing shape[END_REF]: we obtained queens (n=11) and males (n=21) among the first generation from only one colony. We forced brother-sister mating by placing the newly emerged sibling queens and males in a flight cage exposed to natural light (see mating protocol in [START_REF] Lhomme | A method for yearround rearing of cuckoo bumblebees (Hymenoptera: Apoidea: Bombus subgenus Psithyrus )[END_REF][START_REF] Gérard | Discrimination of haploid and diploid males of Bombus terrestris (Hymenoptera; Apidae) based on wing shape[END_REF]. We used these mated queens to produce a second generation of inbred colonies. These queens overwintered during two months prior to a activation by CO2 narcosis (see protocol in [START_REF] Lhomme | A method for yearround rearing of cuckoo bumblebees (Hymenoptera: Apoidea: Bombus subgenus Psithyrus )[END_REF]. Then, we placed the queens in small boxes with two workers and a cocoon to (i.e. to maximize the probability of colony initiation see protocol in [START_REF] Lhomme | Born in an alien nest: how do social parasite male offspring escape from host aggression?[END_REF][START_REF] Lhomme | A method for yearround rearing of cuckoo bumblebees (Hymenoptera: Apoidea: Bombus subgenus Psithyrus )[END_REF]. All 11 mated queens of the first generation initiated a colony, producing a total of 62 workers and 11 males of second generation (Table S1).

In order to produce haploid males, we reared eight queenless micro-colonies with five workers from one colony from bee breeder Biobest N.V. Each micro-colony was maintained in a plastic box in a controlled climate room (same rearing conditions than for colonies). The micro-colonies produced 41 haploid males (Table S1).

We increased the diversity in our dataset by sampling 38 wild males of B. terrestris terrestris in Belgium in natura with a net (Table S1). As wild males of B. terrestris could be confused with the closely related species B. cryptarum, the species identification was checked by CLGS analyses: the two species have species-specific CLGS composition used as diagnostic traits (e.g. Bertsch et al 2005; see CLGS analyses protocol below).

All wild (n=38) and reared (n=52) males (n total=90) were killed by freezing at -20°C. For the reared males, we used five to ten days old males when the B. terrestris males are the most attractive to females [START_REF] Coppée | Age-dependent attractivity of males' sexual pheromones in Bombus terrestris (L.) [Hymenoptera, Apidae][END_REF]. For the wild males, we examined the wing condition (Tkalců, 1969) in order to avoid the sampling of old specimens that have a CLGS modification composition and are less attractive to females [START_REF] Žáček | Comparison of age-dependent quantitative changes in the male labial gland secretion of bombus terrestris and Bombus lucorum[END_REF][START_REF] Coppée | Age-dependent attractivity of males' sexual pheromones in Bombus terrestris (L.) [Hymenoptera, Apidae][END_REF].Determination of the Degree of Ploidy

We assessed the degree of ploidy by genotyped males with 13 microsatellite markers developed for B. terrestris: B11, B100, B121, B126 and B132 [START_REF] Estoup | Characterization of (GT)n and (CT)n microsatellites in two insect species Apis mellifera and Bombus terrestris[END_REF] and 0007_47n22, 0078_59o4, 0222_63d21, 0583_22I4, 0636_34m4, 0646_83e8, 0930_40o1, 0939_33h17 [START_REF] Stolle | Novel microsatellite DNA loci for Bombus terrestris (Linnaeus, 1758)[END_REF]. We extracted sample DNA following the method described by [START_REF] Maebe | Microsatellite analysis in museum samples reveals inbreeding before the regression of Bombus veteranus[END_REF] with some small modifications: 200 μl of 5% Chelex (InstaGeneTM matrix, BioRad) and 10 μl of proteinase K (20 mg/ml) were added to one middle leg and incubated at 56°C for 2h and 97°C for 15 min.The microsatellite loci: B11, B100, B121, B126 and B132 were amplified in multiplex by PCR as described in [START_REF] Maebe | Microsatellite analysis in museum samples reveals inbreeding before the regression of Bombus veteranus[END_REF]. For allele detection of the primers originated from Stolle et al. (2011), we used a tailed-primer approach [START_REF] Schuelke | An economic method for the fluorescent labeling of PCR fragments[END_REF]. We used the PCR protocol described in [START_REF] Maebe | Quantitative trait loci for light sensitivity, body weight, body size, and morphological eye parameters in the bumblebee, Bombus terrestris[END_REF]: simplex PCR reactions were performed, containing 1.5 μl template DNA, 1 μl of 10x PCR buffer (Qiagen), 0.2 μl of 10 mM dNTP's (Qiagen), 0.1 μl of 10 μM forward primer, 0.4 μl of 10 μM reverse primer, 0.4 μl of 10 μM labeled M13-primer and 0.05 μl of 2.5 units/reaction Hotstar Taq DNA Polymerase (Qiagen). Initial denaturing at 95°C for 15 min was followed by 30 cycles of denaturing at 94°C for 30 s, annealing at 48, 52 or 58°C for 30 s, and extension at 72°C for 30 s, and ended at 72°C for 10 min with a final extension step.

Pooled samples were visualized on a ABI-3730xl sequencer (Applied Biosystems) and fragments were scored manually using Peak Scanner Software v 1.0 (Applied Biosystems).

Comparative Chemical Analyses of Cephalic Labial Gland Secretions

We extracted the CLGS of 90 males following the entire cut head protocol described in De [START_REF] De Meulemeester | A simplified protocol for bumble bee species identification by cephalic secretion analysis[END_REF]. Vials containing the solvent and sample were (i) kept for 24 h at room temperature (20°C) to fulfill the extraction and (ii) stored at -40°C for two months prior to the analyses. We (i) analyzed the CLGS samples and identified compounds using gas chromatography-mass spectrometry (GC/MS) and (ii) quantified them by gas chromatography with flame ionization detector (GC/FID) (see GC/MS and GC/FID specifications and temperature program in Lecocq et al., 2015). We identified the CLGS in Xcalibur TM with their mass spectra compared to those in the National Institute of Standards and Technology Library (NIST, U.S.A) with NIST MS Search 2.0. We determined the double bond positions from (i) mass spectra of dimethyl disulphide adducts of unsaturated components [START_REF] Francis | Alkylthiolation for the determination of double-bond position in unsaturated fatty acid esters[END_REF] (reaction time: four hours) and (ii) by chemical ionization with acetonitrile as a reaction gas [START_REF] Oldham | Determination of the double bond position in functionalized monoenes by chemical ionization ion-trap mass spectrometry using acetonitrile as a reagent gas[END_REF]) (see protocol details in [START_REF] Lecocq | Scent of a break-up: phylogeography and reproductive trait divergences in the redtailed bumblebee (Bombus lapidarius)[END_REF]. We quantified the peak areas of CLGS compounds in GCsolution Postrun (Shimadzu Corporation) with automatic peak detection and noise measurement. We elaborated the data matrix (Table S2) with the relative amounts of each compound for each individual (the relative amounts are calculated by dividing the peak areas of compounds by the total area of compounds in each sample). We based our data matrix on the alignment of each compound between all samples performed with GCAligner 1.0 (Dellicour & Lecocq, 2013a,b).

We compared the CLGS between haploid and diploid males with non-metric multidimensional scaling (nMDS) ordination using a Bray-Curtis similarity matrix, three dimensions and 50 runs (R-package ecodist, [START_REF] Goslee | The ecodist package for dissimilarity-based analysis of ecological data[END_REF]. We transformed (log (x+1)) the data matrix to reduce the great difference of abundance between compounds in high and low concentration [START_REF] De Meulemeester | A simplified protocol for bumble bee species identification by cephalic secretion analysis[END_REF]. An hypothetic CLGS differentiation between haploids and diploids was assessed by performing permutation-based version of the multivariate analysis of variance (perMANOVA) using the Bray-Curtis similarity matrix and 10,000 permutations (Rpackage vegan, [START_REF] Oksanen | Vegan: Community Ecology Package[END_REF] (perMANOVA assumption; i.e. similar variance homogeneity between groups was checked by performing a distance-based test for multivariate homogeneity of group dispersions for a one-way ANOVA design; [START_REF] Anderson | Distance-based tests for homogeneity of multivariate dispersions[END_REF]. We determine CLGS compounds that are potentially specific to and regular to diploid or haploid males with the indicator value (IndVal) method [START_REF] Dufrêne | Species assemblages and indicator species: The need for a flexible asymmetrical approach[END_REF]; see details in Lecocq et al., 2015b). All statistical analyses were performed in R (R Development Core Team, 2013).

Results

Genetic analyses showed that reared queen-right colonies produced nine diploid males (Table S1). One diploid male was also found in wild specimens (Table S1). Other males were unambiguously haploid (Table S1).

Fifty six compounds were detected in the CLGS of analyzed specimens (Table S2). The general pattern obtained for CLGS of diploid males corresponded to those of haploid males. No qualitative difference was detected between haploid and diploid males (including no differences in double bound positions). The main compounds (i.e. compounds that have the highest relative amount in one individual at least) were farnesol, 2,3-Dihydrofarnesol, geranylcitronellol, and 2,3-dihydrofarnesyl dodecanoate. The nMDS (stress value = 0.16; R2=0.83) did not show any obvious CLGS differentiations between haploid and diploid males except for a slight trend observed for the reared diploid males in the nMDS (i.e. these specimens are more or less grouped together at the edge of the nMDS) (Fig. 1). The perMANOVA (DF = 1, F = -9.17, P =0.79) analyses confirmed the lack of significant CLGS differentiation between the two groups. The IndVal method did not reveal indicator compounds with strong significance (IndVal>0.70) in haploid or diploid males (Table S2).

Discussion

Diagnostic Features of Diploid Bumblebee Males

In view of the dramatic consequences of diploid male production on population (Zayed & Packer, 2005), we expected that B. terrestris would display pre-mating mechanisms avoiding virgin queen mating with diploid males as already observed in some Hymenoptera [START_REF] Herzner | Brothers smell similar: variation in the sex pheromone of male European Beewolves Philanthus triangulum F. (Hymenoptera: Crabronidae) and its implications for inbreeding avoidance[END_REF]. Such pre-mating mechanisms are most likely based on differences between reproductive traits used in pre-mating recognition between diploid and haploid males. Therefore, we investigated one reproductive trait, the CLGS that acts as a key trait for the bumblebee pre-mating recognition. However, our comparative statistical analyses based on CLGS composition did not detect obvious differences between diploid and haploid males in their CLGS, a main key trait for mate attraction (Fig. 1). The slight trend observed for the reared diploid males in the nMDS (i.e. these specimens are more or less grouped together at the edge of the nMDS; Fig. 1) could reflect a "family-effect" since all queens that produced these males are sister mated with their brothers as already observed in other study [START_REF] Lecocq | Molecular and chemical characters to evaluate species status of two cuckoo bumblebees: Bombus barbutellus and Bombus maxillosus (Hymenoptera, Apidae, Bombini)[END_REF]. Nevertheless, our statistical analyses did not show significant differentiation between CLGS of diploid and haploid males (see perMANOVA results). Even if the components of the male CLGS that constitute the pheromones are not known (i.e. those components that release a female response or convey information for females), this does not invalidate the basic result that there is no significant differentiation in diploid male CLGS since there is no specific and regular chemical compounds to diploid or haploid males (see IndVal results in Table S2). Although we are aware of the limited power of our statistical analyses due to our small sample size (10 diploid males), our data clearly demonstrates that there is no strong difference in the CLGS composition between haploid and diploid males.

According to our results we speculate that the queens are not able to discriminate the diploid males through their CLGS. However, since the bumblebee male pre-mating signal [START_REF] Baer | Bumblebees as model organisms to study male sexual selection in social insects[END_REF][START_REF] Ayasse | Chemical ecology of bumble bees[END_REF][START_REF] Streinzer | Functional morphology of the visual system and mating strategies in bumblebees (Hymenoptera, Apidae, Bombus )[END_REF] includes behavioral and other chemical features such as tendon glands [START_REF] Jarau | Leg tendon glands in male bumblebees (Bombus terrestris): structure, secretion chemistry, and possible functions[END_REF] that are mainly poorly known, we cannot rule out the fact that the whole pre-mating signal of diploid males are evaluated by females as being strictly similar to those of haploid males. Indeed other traits involved in the species communication (e.g. cuticular hydrocarbons) or in the reproduction could be different in diploid B. terrestris males such as observed for cuticular compounds and testis size of diploid drone in the close relative Apis mellifera L. [START_REF] Herrmann | Characters that differ between diploid and haploid honey bee (Apis mellifera) drones[END_REF]. Moreover, diploid B. terrestris males are distinct in other characters such as wing shape, body size, and physiology [START_REF] Duchateau | Sexual biology of haploid and diploid males in the bumble bee Bombus terrestris[END_REF][START_REF] Gerloff | Effects of inbreeding on immune response and body size in a social insect, Bombus terrestris[END_REF][START_REF] Gérard | Discrimination of haploid and diploid males of Bombus terrestris (Hymenoptera; Apidae) based on wing shape[END_REF]. These features could allow the female to diagnose diploid males and/or decrease the fitness of diploid males (e.g. fitness of a bumblebee male is correlated with body size, Cueva del [START_REF] Del Castillo | Macroevolutionary patterns of bumblebee body size: detecting the interplay between natural and sexual selection[END_REF].

Further studies on other traits and ethological test allowing choices of virgin females could help to define whether the B. terrestris females are not able to recognize diploid males. Nevertheless, current evidence suggests that queens are not able to discriminate the diploid males during the pre-mating recognition because (i) our results show that there are no diploid specific cephalic labial gland secretions that presumably contain the main long-distance signal of the bumblebee pre-mating recognition (e.g. [START_REF] Ayasse | Chemical ecology of bumble bees[END_REF], and (ii) diploid males have the same mating success as haploid males in experimental conditions where short distance pre-mating recognition signals should act [START_REF] Duchateau | Sexual biology of haploid and diploid males in the bumble bee Bombus terrestris[END_REF].

Diploid Males in the Wild

Diploid males have been previously detected, mainly in bumblebee species with small and fragmented distribution like B. florilegus in Japan and B. muscorum in Great Britain, with a considerable proportion of such males of up to 50% [START_REF] Darvill | Population structure and inbreeding in a rare and declining bumblebee, Bombus muscorum (Hymenoptera: Apidae)[END_REF][START_REF] Takahashi | Diploid male production in a rare and locally distributed bumblebee, Bombus florilegus (Hymenoptera, Apidae)[END_REF]. However, diploid males were detected in Table Table S1 Sampling table. The table gives the origin of samples and the ploidy determination by genotyping each specimen with 13 microsatellite loci with indication of the fluorescent dye used in the PCR amplification.

Table S2

Data matrix of cephalic labial gland secretions (relative amounts of each compound) and list of the identified compounds in the sampling. MW is the molecular weight of the compound. IndVal results are indicator value obtained in IndVal method.

Figure captions

Fig. 1

 1 Fig. 1 Statistical analyses of cephalic labial gland secretion. Non-metric multidimensional scaling (nMDS) ordination analysis of chemical reproductive trait of diploid and haploid Bombus terrestris males (n= 90). Three first axes of the nMDS ordination plot based on Bray-Curtis distances calculated on a male cephalic labial gland secretion matrix of B terrestris.

Acknowledgments

This research was supported by the Belgian Science Policy (project BR/132/A1/BELBEES) and the Institute of Organic Chemistry and Biochemistry of the Academy of Sciences of the Czech Republic (RVO61388963). TL was grant holders of the Fonds pour la Recherche dans l'Industrie et l'Agriculture. Authors acknowledge David Notton (the Natural History Museum, London, UK) for correcting the English.

populations of widespread species like B. occidentalis, B. perplexus, and B. terricola in about 5% frequency [START_REF] Whidden | Frequencies of Diploid Males in Natural Populations of Three North American Bumble Bee (Bombus) Species (Hymenoptera: Apidae)[END_REF]. Our results show a ratio (≈ 2%) of wild diploid males in B. terrestris similar to those found in other widespread bumblebee species. This low ratio was expected since B. terrestris is one of the most common and widespread bumblebee species in Europe [START_REF] Rasmont | An overview of the Bombus terrestris (L. 1758) subspecies (Hymenoptera: Apidae)[END_REF]. Indeed, the large distribution and population density of this species [START_REF] Rasmont | An overview of the Bombus terrestris (L. 1758) subspecies (Hymenoptera: Apidae)[END_REF][START_REF] Rasmont | Climatic risk and distribution atlas of European bumblebees[END_REF] make unlikely a limited genetic flow (i.e. genetic drift) between populations and therefore a high diploid ratio. However, our current sampling is limited to a small geographic region and to a small sample. Further studies on the diploid male ratio on large B. terrestris sampling and other European and non-European species are needed to validate this hypothesis.

Consequences for Bumblebee Conservation

If future studies confirm current evidence (see our results and Duchateau & Marien 1995) that diploid males cannot be recognized through their reproductive traits by females, this rules out the potential role of pre-mating recognition in curbing an ongoing vortex of extinction of diploid males in bumblebee species. Zayed & Packer (2005) and [START_REF] Hein | What stops the "diploid male vortex"?-A simulation study for species with single locus complementary sex determination[END_REF] underline four crucial conditions that must be met to initialize the diploid male vortex: (i) fixed sex ratio, (ii) low reproductive rate, (iii) no mate choice (i.e. females cannot detect if potential mates are related mates), and (iv) extreme population fragmentation. Conservation plans should act on the two last factors since the two first conditions depend on intrinsic features of bumblebee species. The decrease of population fragmentation through the establishment of corridors [START_REF] Holzschuh | Grass strip corridors in agricultural landscapes enhance nest-site colonization by solitary wasps[END_REF] between populations could increase the genetic diversity of populations and thus prevent the initialization of diploid male vortex. Alternatively, individual translocations between populations [START_REF] Arrendal | Genetic evaluation of an otter translocation program[END_REF] could neutralize the fragmentation and enable females to choose exotic (unrelated) mates. Nevertheless, the beneficial effects of such translocations would appear to be restricted to areas in the immediate vicinity of the release sites [START_REF] Arrendal | Genetic evaluation of an otter translocation program[END_REF]. Moreover, potential translocations should consider the interpopulational differentiations to maximize the genetic distinctiveness of translocated individuals while preserving the species diversity (i.e. avoid genetic homogenization through translocation; Lecocq et al., 2015b).

Disclosure

The authors declare that they have no conflict of interest.


Kaidi Wan

Luc Vervisch

Jun Xia
email: jun.xia@brunel.ac.uk

Pascale Domingo

Zhihua Wang
email: wangzh@zju.edu.cn

Yingzu Liu

Kefa Cen

Alkali metal emissions in early stage of a pulverized-coal flame: DNS analysis of reacting layers and chemistry tabulation

Keywords: Pulverized-coal combustion, Direct numerical simulation, Sodium chemistry, Alkali metal, Chemistry tabulation

The intricate coupling between coal pyrolysis, gas phase combustion and the emissions of alkali metal, such as sodium, is studied in early stage of a temporally evolving three-dimensional planar turbulent jet carrying pulverized-coal particles. Complex chemistry is used to account for both the combustion of volatile hydrocarbons and the sodium containing species. The response of the sodium chemistry is analyzed in the mixture fraction space, along with the topology of the reactions zones. Combustion is found to start preferentially in partially premixed flames, which then evolve toward diffusion-like reactive layers and reach chemical equilibrium. From the direct numerical simulation (DNS) database, the possibility of modeling the dynamics of sodium species using one-dimensional premixed flamelet generated manifolds (FGM) is investigated. A chemical lookup table is constructed for the combustion of the partially premixed volatiles and an additional three-dimensional simulation is performed to compare the tabulated sodium species against their reference counterparts with complex chemistry. Quantitative analysis of the performance of the developed chemistry tabulation confirms the validity of the approach. Perspectives for the modeling of sodium emissions in pulverized-coal furnaces and boilers are finally drawn.

Introduction

Solid fuels, such as coal and biomass, are widely utilized to support the worldwide energy consumption due to their broad availability and the overall flexibility of solid fuel combustion systems [START_REF]Key world energy statistics[END_REF]. In practice, the presence of alkali metals such as sodium (Na) in coal and potassium (K) in biomass leads to massive ash depositions and corrosions in the solid-fuel-fired boilers [START_REF] Niu | [END_REF]. These alkali-induced issues severely limit the utilization of coal and biomass with a relatively high concentration of alkali, such as Zhundong coal in China [3]. Therefore, a better understanding is urgently required of the fundamental mechanisms driving the formation and transformation of alkali species and their interactions with complex multi-phase turbulent reacting flows during coal and biomass combustion.

Over the past few decades, experimental research on sodium release and reacting dynamics evolves from offline to online measuring techniques using advanced laser diagnostics, e.g., [4][5][6]. On the numerical side, a one-step Arrhenius model has been firstly proposed by van Eyk et al. [7] for sodium release during the combustion of a single char pellet. In recent works [5,6], a twostep kinetics model, which is able to describe sodium release in all stages of coal combustion, has been developed based on the simultaneous measurements of the sodium release, surface temperature and diameter of a burning coal pellet. These sodium release models are suitable to be applied in circulating fluidized bed (CFB) combustion, since the coal pellet size and burnout time employed are of the order of magnitude of those of a typical CFB boiler [6]. However, the sodium release model of pulverized-coal combustion (PCC) must still be developed.

Considering the homogeneous reactions of sodium species, well-established chemical reaction mechanisms [8] are available to model the final forms of alkali species in post-combustion gases. To the best of our knowledge, however, the dynamic of sodium release and the reaction characteristics in turbulent pulverized-coal combustion have not been investigated yet. Flow-resolved or quasi-direct numerical simulation of pulverized-coal combustion (PCC-DNS) has been attempted recently [9,10], as the sharply increasing computing capacity now allows for investigating the complex multi-physics PCC phenomenon. Since the turbulence-chemistry interaction is not modeled but resolved in these simulations, the results provide precious insights into the combustion characteristics of turbulent pulverized-coal flames and also serve as reference for subgrid scale modeling development.

Since its pioneer introduction in Reynolds-Averaged Navier-Stokes (RANS) modeling of turbulent premixed flames [11], the tabulated chemistry methods, e.g., flamelet generated manifolds (FGM) [12,13], have shown great potential on predicting major hydrocarbon and also minor species, such as NO x [14,15].

Within this context, the aim of the present study is twofold. First, the responses of sodium species in early stage of a three-dimensional (3D) turbulent pulverizedcoal flame are simulated and analyzed, with char reaction ignored. Second, the possibility of modeling the reaction dynamics of these chemical species using onedimensional (1D) premixed FGM is examined. In a previous study [16], a tabulation method of sodium species in PCC was proposed assuming the sodium reactions to take place only once the hydrocarbon combustion has reached its equilibrium state. Here, the coupling between sodium reactions and hydrocarbon combustion is directly accounted for in the three-dimensional simulations, in which the flow and the flame are fully resolved. The coal particles are modelled as Lagrangian point sources, with a two-way coupling with the flow (mass, momentum and energy). In the vast literature on multi-phase flow, these simulations are sometimes called 'quasi-DNS' or 'point-particle DNS'. To ease the reading, our simulation will be denoted 'DNS' in the following, still keeping in mind that the internal structure and the boundary layers of the coal particles are not resolved.

Pulverized-coal flame configuration and numerics

A three-dimensional temporally evolving planar pulverized-coal jet flame is studied (Fig. 1). The physical dimensions of the computational domain composed of 8.38M cells are L x = 25.6 mm, L y = 25.6 mm and L z = 12.8 mm. A uniform mesh of h = 100 µm is employed and the time step is 4.8 µs (CFL < 0.8) with subiteration in a stiff chemistry dedicated solver, which is further discussed below ( § 2.4) in terms of flame resolution. Initially, high-speed primary air (bulk velocity: 15 m/s, 300 K) ladened with coal particles is introduced for |y| < 2.5 mm. This flow is surrounded by a low-speed hot coflow of burnt gases (3 m/s, 1559 K). The locations of the coal particles follow a random uniform distribution, with a particle density of 1400 kg/m 3 . The initial number of particles is 21945 with a mono-disperse distribution of diameter of 25 µm. The coal properties are those of the Loy Yang brown coal, for which coal analysis data can be found in [16,17]. The precentage of sodium that is releasable during the coal pyrolysis stage Figure 1: Instantaneous pulverized-coal particle distribution and iso-Q criterion [19], Q = 10 5 s -2 , colored by gas temperature. The coal particles are colored by the burnout, which denotes the ratio of the released volatile mass to the initial mass of volatile in a pulverizedcoal particle.

is set to 19.1%, according to [17]. The coflowing stream carries the burnout products of the same coal volatiles at an equivalence ratio of 0.45. This stream promotes the ignition of the coal particles, as mixing with burnt products would do in a real burner. Periodic boundary conditions are employed in all directions. Turbulent fluctuations of 1.2 m/s (8% of the primary air velocity) are initially added within the shear layers between the air jet and the coflow, to facilitate the flow destabilisation. The parameters of the present PCC case, e.g., the coflow temperature and mass loading ratio of coal particles, are chosen according to jet flame studies [16,18], with adjustments to ease DNS but without change to the multiphase reacting flow dynamics.

Coal pyrolysis, including sodium release, gas phase volatile combustion and sodium reactions are simulated. Heterogeneous char reaction is not considered here, because the contribution of char reaction should be weak in such a small-scale PCC flame, as already discussed in [13,16,20]. In addition, char reaction is much slower than homogeneous gas phase combustion and, under the present conditions, sodium reactions are found to be completed within these shorter time scales.

Fig. 1 shows an instantaneous snapshot of the planar pulverized-coal flame. The turbulent flow structures are visualized via the iso-surface of Q criterion (second invariant of the velocity gradient tensor) [19] = 10 5 s -2 and colored by gas temperature. The instantaneous distribution and burnout of each particle are also illustrated.

Gas phase and pulverized-coal particles

The governing equations for the gas phase and the pulverized-coal particles are solved in the Eulerian and Lagrangian frameworks, respectively, with a low-Mach number in-house code specialised in the simulation of turbulent coal flames [16,18,21]. The conservation equations for mass, momentum, species and temperature of the gas phase are those of [18], in the present study without the subgrid scale modeling terms. The Lewis number (Le = 1.0) and the Prandtl number (Pr = 0.7) are assumed to be constant.

The momentum and the temperature equations of the Lagrangian coal particles are also those of [18]. Radiative heat transfer is solved with the Discrete Ordinates Method (DOM) [START_REF] Chandrasekhar | Radiative transfer[END_REF]. The particle emissivity ε p is set to 0.9 [START_REF] Wan | [END_REF]. The gas absorption coefficient is determined by the weighted-sum-of-the-gray-gases model (WSGGM) [24]. Finally, the mass loss rate of each coal particle (dm p /dt), due to the pyrolysis, is predicted by the single first-order reaction model (SFOM) proposed by Badzioch and Hawskley [25]. The kinetic parameters of the SFOM model, i.e., A v (9×10 5 s -1 ), E v (6×10 4 J/mol) and Q v (1.05), have been calibrated by the chemical percolation devolatilization (CPD) model, as in [21].

Gas phase chemistry

The volatiles react according to the previously validated DRM22 skeletal mechanism by Kazakov and Frenklach [START_REF] Kazakov | Reduced reaction sets based on grimech 1.2[END_REF], involving 104 elementary reactions and 22 chemical species. The compositions/mass fractions of the volatile fuels, including CH 4 /0.037, C 2 H 2 /0.309, CO/0.208, CO 2 /0.155, H 2 /0.030, H 2 O/0.261, are obtained from the CPD model [18,21]. The Tar species is replaced by C 2 H 2 [20]. The compositions predicted by the CPD model are slightly adjusted to fulfil the elemental mass conservation.

The detailed reaction mechanism of alkali metal species by Glarborg and Marshall [8] includes 105 elementary reactions over the elements Na, K, C, H, O, S and Cl. In the present study, only the subset with the elements Na, C, H, and O is considered, which includes five sodium species, i.e., Na, NaO, NaO 2 , NaOH and Na 2 O 2 H 2 , and 24 elementary reactions. Considering that the sodium vapor generated inside the porous structure of a coal particle will be transported outward by the volatile yielded during the pyrolysis stage and that the sodium release was found proportional to the burnout of a coal particle during the early combustion stage [6], the sodium release rate is assumed to be proportional to the volatile release rate [16]. As in our previous study [16], atomic sodium is the alkali species released from coal in the simulation.

Numerical schemes

The numerics is based on an approach previously used for both DNS and large eddy simulation (LES) [START_REF] Pierce | [END_REF]28]. The time advancement uses a second-order Crank-Nicolson scheme. A Quadratic Upstream Interpolation for Convective Kinematics (QUICK) scheme is used for the scalar advection terms in the species and temperature equations, while a second-order central difference scheme is applied to the scalar diffusion terms in the species and temperature equations and all terms in the momentum equation. An Alternating Direction Implicit (ADI) method has been employed, in which semiimplicit tridiagonal/pentadiagonal equations are solved separately for each direction. The particle equations are explicitly advanced using a second-order Runge-Kutta (RK2) scheme.

Verification of the grid resolution

Before proceeding with the analysis, the resolution level of the scalar signals within the reaction zones is analyzed. Figure 2 shows the comparison of O radical profiles from the DNS and from one-dimensional counterflow diffusion flames computed for various strain rates with CANTERA [START_REF] Goodwin | Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF], for the present volatile fuel and oxidizer conditions. The DNS profile is extracted from the location where the mixture fraction scalar dissipation rate, χ Z , reaches its maximum, which matches here the maximum burning rate. According to non-premixed flame theory, the thinner reaction zones should be found for such intense burning rates and high level of χ Z [START_REF] Peters | Turbulent combustion[END_REF]. It can be observed in Fig. 2 that the reaction zone is well resolved, as there are more than 10 grid points inside this thinner flame layer. In addition, the maximum scalar dissipation rate in the DNS (χ Z = 937 s -1 ) is found well below the quenching point (χ Z = 3124 s -1 ). Interestingly, the Y O DNS signal also fits within the expected ranges observed in the diffusion flamelet computations. A similar analysis has been carried out with the partially premixed flames introduced thereafter for chemistry tabulation. The preheating of the reactants combined with the overall low levels of equivalence ratios due to the pyrolyzed character of the fuel make the flame thick, and a similar conclusion is reached concerning the required resolution.

Analysis of flame and sodium emission dynamics

Flame regimes

The combustion characteristics of the pulverized-coal flame are studied using the flame index, F.I. [9,[START_REF] Domingo | [END_REF] fractions of the volatile fuel and the oxidizer as: F.I.

= ∇Y f • ∇Y O 2 , where Y f = Y CH 4 + Y CO + Y C 2 H 2 + Y H 2 .
Positive values of F.I. are likely to indicate premixed flame regimes, while negative ones probe diffusion flame regimes. Figure 3 shows the instantaneous distributions (plane of z = 0) of F.I. in the regions where the heat release rate exceeds 10 7 W/m 3 and therefore where the combustion actively happens. The maximum heat release rate in the 3D domain is 2.37×10 10 W/m 3 . First, the coal particles in the shear layers are heated up by the high-temperature coflow. Pyrolysis takes place and the volatile starts to escape from the particles to mix with the surrounding air. The first ignition is thus controlled by homogeneous combustion, with spots dominated by premixed regimes (t 3 ms). However, in the region around particles with the highest temperatures, the diffusive combustion mode also appears because of higher fluxes of volatile release, which allow for the fuel to accumulate before burning against air (t = 2.5 and 3 ms). Later on, more and more particles take part in pyrolysis due to turbulent mixing and the volatile mass fraction further increases, leading to a dominant diffusion combustion mode (t = 5 and 10 ms). A global observation that is now refined from a direct analysis of the scalar fields.

Sodium species dynamics

Figure 4 shows the scatter plots of instantaneous mass fractions of the five sodium species against the mixture fraction, Z, at t = 5 ms. The mixture fraction is defined as Z = 1.0 -Y N 2 /0.767 (Z = 1 in the pure volatile fuel, Z = 0 in the oxidizer and Z = 0.068 in the coflow). The data is obtained from the entire 3D domain. In agreement with our previous study [16], Na and NaOH are the two major sodium species. Y Na stays almost zero for Z < 0.05, to then increase almost linearly with the Interestingly, the sodium chemistry response differs from the one of usual hydrocarbon diffusion flames. For instance at high-strain rate (red line) in the 1D counterflow flame, fast mixing promotes low temperature sodium reactions, with the production of NaOH on the rich side and a decrease of Na concentration at Z = 0.4 (Fig. 4). Compared with the sodium profiles obtained from 1D counterflow diffusion flames, many points of the turbulent flame fall outside the region between the profiles corresponding to the lowest and highest strain rates (Fig. 4). Because the maximum scalar dissipation rate in the DNS stays below the quenching limit of steady flamelets (Fig. 2), these points are likely to be representative of unsteadiness and/or partial premixing of the reactants. Accordingly, the tabulation of sodium species will now be attempted using premixed flamelets rather than steady diffusion flamelets. Indeed, the fresh mixture gets partially premixed quite rapidly after be-ing ejected from the particles and the gradients of thermochemical quantities across iso-equivalence ratio surfaces stay moderate compared to gradient in other directions (as also previously reported in the literature at many places for for fuels released from a disperse phase, see for instance [32].)

Premixed flamelet generated manifold for sodium

Lookup table construction

The FGM approach was first applied in pulverizedcoal combustion by Knappstein et al. [13] recently. Here it is used to tabulate the sodium chemistry. A set of one-dimensional freely propagating and laminar premixed flames are computed with CANTERA [START_REF] Goodwin | Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF], as in the DNS with the DRM22 mechanism [START_REF] Kazakov | Reduced reaction sets based on grimech 1.2[END_REF] for the combustion of volatiles and the detailed sodium reaction mechanism [8]. Volatile species including atomic sodium are mixed with air at various equivalence ratios. The reaction progress variable is defined as

Y c = Y H 2 O +Y CO +Y CO 2 .
To build the lookup table, the flamelet computation is repeated over a range of equivalence ratios, mapped in the mixture fraction Z, and also varying total enthalpy levels, with the Lewis number chosen at unity for all the species. The responses in the physical space are then remapped (x = x(Y c )) into the progress variable space (Y i (x) = Y i (Y c )), for every mixture fraction and total enthalpy levels (H), leading to Y i = Y i (Z, Y c , H). The progress variable Y c is normalized by its bounds versus equivalence ratio before being stored in the table, to define a normalized progress variable C (Y c is solved with the flow not C). When Z is outside the flammability limits, linear interpolations between chemical equilibrium and the mixing line are performed, to obtain the required information during the construction of the table. Hence, the final lookup table is parameterized by three coordinates: the mixture fraction, the total enthalpy and the progress variable. The manifold is mapped on to 101 × 8 × 101 data points for Z × H × C. Figure 5 shows the responses of a major sodium species, Y NaOH , to Z and Y c . To verify the lookup table discretization and the multi-linear interpolations, as in [16], one-dimensional premixed flames were first solved with the progress variable and the lookup table, selecting conditions between tabulated ones, and the results compared with success against the detailed chemistry solution.

Comparison against DNS

The reference three-dimensional DNS including the two kinetics (the DRM22 mechanism for hydrocarbon and the transported Na species) is complemented by a simulation in which the DRM22 is kept, but the FGM lookup table replaces the transport by the flow of the Na species. The thermal effect induced by the coflowing stream of burnt gases is accounted for through the tabulation versus total enthalpy, while the dilution of Na by those inert burnt gases is calibrated with an additional passive scalar. The results closely follow those of the DNS for both major and minor sodium species (Fig. 6). The hydrocarbon species are almost unchanged, because of the very small concentrations of Na species, which do not affect much the hydrocarbon kinetics (not shown for brevity).

Figure 7 shows the scatter plots of instantaneous mass fractions of the five sodium species predicted by the FGM tabulation against those by DNS at t = 10 ms. For the two major sodium species, i.e., Y Na and Y NaOH , the conditional means of the FGM predictions agree well with the DNS results. Among the other minor sodium species, the FGM conditional means of Y Na 2 O 2 H 2 match the DNS results well with a slightly large uncertainty. For Y NaO and Y NaO 2 , the FGM conditional means and DNS results agree well, except that the FGM tends to slightly under predict them under a high concentration of Y NaO and Y NaO 2 . It can be found that accurate predic- tion of minor sodium species with FGM is more difficult, since the mass fractions of minor sodium species are much lower and they are more sensitive to the temperature and compositions of the gas phase. To better quantify the performance of the lookup table, the Pearson correlation coefficient (r ∈ [0, 1]) [33] is computed for the predictions of the five sodium species in DNS and FGM. As shown in Table 1, r > 0.9 is achieved for all the sodium species at various time instants, considering the grid points where Y sodium > 0 of the entire domain, thus confirming the potential of the proposed strategy.

Summary

A temporally evolving planar pulverized-coal jet flame is simulated with complex chemistry, to account for both the combustion of volatile hydrocarbons [START_REF] Kazakov | Reduced reaction sets based on grimech 1.2[END_REF] and the reactions of sodium species [8], which are released from the burning coal particles. From the flame index, combustion is found to start preferentially in partially premixed flames during the ignition stage, to then evolve toward diffusion-like reactive layers where the diffusive combustion mode dominates, with however partial premixing of the reactants. Among the five sodium species, Na and NaOH are found to be the two major ones. Y Na increases along with the mixture fraction Z; while Y NaOH feature a complex behaviour, as it reaches higher values at both fuel-lean and fuel-rich sides, to drop to a much lower level under the stoichiometric condition.

A sodium chemistry lookup table is constructed from laminar one-dimensional premixed flames at various mixture fractions and enthalpy levels, thus following an FGM tabulation strategy. Three-dimensional simulations are performed with the lookup table for sodium species, keeping the remaining parts of the modeling unchanged. Both major and minor sodium species predicted by FGM tabulation method closely follows the reference simulation, with a Pearson correlation coefficient larger than 0.9. Prof. Peter Glarborg of DTU, who provided us the detailed mechanism of alkali metal species. Computing resources were provided by CRIANN (http://www. criann.fr).

t = 2 Figure 3 :

 23 Figure 3: Instantaneous distributions of flame index at the plane z = 0, at locations where heat release rate exceeds 10 7 W/m 3 . Coal particles in the region of |z| 100 µm are shown and colored by particle temperature.

Figure 4 :

 4 Figure 4: Scatter plots of instantaneous mass fractions of the five sodium species against Z at t = 5 ms. Blue thick line and red thin line are the sodium profiles obtained from 1D counterflow diffusion flames, corresponding to the lowest (37 s -1 ) and highest strain rates (5246 s -1 ).

Figure 5 :

 5 Figure 5: FGM lookup table of Y NaOH vs mixture fraction and progress of reaction. Baseline enthalpy level, i.e., no heat loss.

Figure 6 :

 6 Figure 6: Comparison of Na, NaOH and NaO predicted by the DNS (top) and the FGM (bottom) sodium tabulation approach at t = 10 ms.

Figure 7 :

 7 Figure 7: Scatter plots of instantaneous mass fractions of the five sodium species predicted by the FGM tabulation method against those by DNS at t = 10 ms. Red line: conditional mean and fluctuation of the FGM predictions. Blue line: Perfect correlation.

  . It is calculated from the spatial gradients of the mass

		0.004
				max scalar dissipation rate (1/s):
					15
		0.003			29 58
					112
					220
	O (-) Y	0.002			428 833 1637
		0.001			3124 937 (DNS)
		0.000	-0.001	0.000	0.001	0.002
				Distance (m)
	Figure 2: Comparison of O radical profiles between DNS results and
	1D counterflow diffusion flames at various strain rates. The maximum
	scalar dissipation rate of each profile is denoted.

Table 1 :

 1 Pearson correlation coefficients of the five sodium species predicted by the DNS and the FGM sodium chemistry table at different time instants.

	Time	Na	NaOH Na 2 O 2 H 2 NaO	NaO 2
	3 ms	0.989 0.973	0.947	0.963 0.969
	5 ms	0.999 0.979	0.912	0.977 0.969
	10 ms 0.997 0.971	0.923	0.981 0.948

Acknowledgments

This work was jointly supported by the National Natural Science Foundation of China (51706200, 51390491, 51422605), the Project funded by China Postdoctoral Science Foundation, the Royal Society and the Engineering and Physical Sciences Research Council (EPSRC) of the UK. Special thanks are due to


Bastien Duboc

Guillaume Ribert
email: guillaume.ribert@coria.fr

Pascale Domingo

Evaluation of chemistry models on methane/air edge flame simulation

Keywords: Hybrid chemistry, Compressible solver, Edge flame, DNS

The integration of chemistry into a numerical fully compressible solver is carried out in this study using three models: detailed chemistry, fully tabulated chemistry (CTC) and a model coupling both approaches called HTTC, for hybrid transported-tabulated chemistry. With HTTC major species are transported while most minor species are tabulated. As minor species are no longer transported with the flow, the time step is close to the values usually encountered for non-reactive flows, far beyond what is found in detailed chemistry. The performance of HTTC for reproducing the dynamics of a methane/air edge flame featuring a very strong mixture fraction gradient is also investigated. The results agree favorably with the reference case simulated with detailed chemistry unlike the CTC model which is unable to predict the topology of the flame. Finally, the shape of the flame, the flame speed and the flame stabilization height are reasonably well captured with HTTC with a calculation cost divided by about 5 compared to the reference case.

• Corresponding author:

-Last name: Ribert -First name: Guillaume -Mailing address:

CORIA -CNRS and INSA de Rouen, Avenue de l'Université, 76800 Saint-Etiennedu-Rouvray, France -Phone number: +33(0)2 32 95 97 92 -Email: guillaume.ribert@coria.fr

• Colloquium 4: Laminar flames

• Abstract: 160 words

• Article, total length: 6000 = 6 × 900 + 2.2 × 270 mm with the two-column method (LaTeX)

• Color figures: not applicable

Introduction

Introducing fully detailed schemes in numerical simulations for combustion is still a scientific challenge. In practical situations, it cannot be achieved because the number of species and reactions are way too large [1]. Instead of detailed chemistry (acronym FTC), reduced kinetics or tabulated thermochemistry (CTC) may be used [2]. The first strategy may lead to inaccuracy if the reduced scheme is not optimized [3] because minor species and radicals are missing. If done carefully, the resulting reduced schemes compare well with the original detailed mechanism on the dedicated range of validity. The second approach is based on the tabulation of chemical responses of canonical combustion problems such as one-dimensional laminar premixed flames [4,5]. These structures are often projected into a progress variable and mixture fraction space to build a look-up table. Thus, only these two variables need to be transported with the flow, dramatically reducing the computational cost. However, such tables are cumbersome to create, lack of flexibility and may lead to very large database not suited to the context of high-performance computing. Hence, table downsizing methods have been discussed in the literature, using the self-similarity behavior of the radical species in laminar flamelets [6,7] or ignition phenomena [8,9]. This property has been further exploited by Ribert et al. [10] to develop a strategy combining the detailed-chemistry solving for the main species with the tabulation of the intermediate species, called Hybrid Transported-Tabulated Chemistry (HTTC). CTC and HTTC models are presently evaluated on the challenging configuration of a methane/air edge flame featuring a large gradient of mixture fraction, and compared with results coming from detailed chemistry. In such configurations, the reactants are partially premixed before burning, and an edge flame is present consisting of a premixed flame front divided into a lean and a rich zone, followed by a trailing diffusion flame, that burns the excess of the reactants downstream [11]. The whole range of equivalence ratio (φ) from pure fuel to pure oxidizer is then present making the simulation challenging for any combustion models.

Chemistry modeling and numerical solver

Dealing with a chemical system composed of

N S species (S = {1, • • • , N S }) reacting through N R reactions (R = {1, • • • , N R }
) requires to write one transport equation for each species k ∈ S when using a FTC solver:

∂ρY k ∂t + ∂ ∂x i ρ[u i + V k,i + V c i ]Y k = ωk , k ∈ S . (1)
x i , t and u i are the spatial coordinates, time and ithvelocity components, respectively. ρ is the density defined as ρ = S ρ k , with ρ k = ρY k . Y k is the mass fraction of species k with S Y k = 1, V k,i is the diffusion velocity of species k computed with the Hirschfelder and Curtiss approximation [START_REF] Hirschfelder | Molecular Theory of Gases and Liquids[END_REF], V c i is the correction velocity to ensure the mass conservation and ωk is the chemical source term of species k with S ωk = 0. At the extreme opposite, CTC methods [START_REF] Vicquelin | [END_REF]14] assume that chemical evolutions in the composition space are parameterized by a reduced set of N t variables such as the progress variable, Y c , the mixture fraction, Z, enthalpy, etc. that are transported with the flow. The N S transport equations used with the FTC solver are replaced in CTC by N t ( N S ) equations plus a look-up table containing all the expected flame structures. The balance equation for Y c is formally written:

∂ρY c ∂t + ∂ ∂x i (ρu i Y c ) = ∂ ∂x i Å ρD Y c ∂Y c ∂x i ã + ωc , (2)
with

Y c = k∈S α k Y k = Y CO 2 + Y CO in
D Y c = k∈S α k ∂Y k ∂Y c D k = ∂Y CO 2 ∂Y c D CO 2 + ∂Y CO ∂Y c D CO , (3)
to ensure a proper flame speed. ωc and D Y c are part of the look-up table and depend on N t variables. Z is transported with a unity Lewis number assumption as in [15]. With HTTC [10], the whole kinetic scheme is kept unaltered meaning that the knowledge of the N S species mass fractions is required, but the set of chemical species is splitted in a set of major

(M = {1, • • • , N M }) and minor (m = {1, • • • , N m }) species: N S = N M + N m .
Major species are transported with the flow and the mass fractions of minor species come from generic laws observed in the simulations of canonical problems. For minor species, the self-similar flame tabulation (S2FT) technique [6,7,8] is presently exploited. This reduced look-up table is accessible by Y c and Z (or φ) [16]. Eq. ( 1) is then replaced by

∂ρY k ∂t + ∂ ∂x i ρ[u i + V k,i + V c i ]Y k = ωk , k ∈ M (4)
plus a look-up table (S2FT),

Y k = f (Y c , Z), k ∈ m. (5)
In Eq. ( 4), the computation of ωk requires the evaluation of the rate-of-progress of all elementary reactions of the kinetic scheme. For any species k ∈ m, Y k is directly evaluated with Eq. ( 5) and not through a transport equation anymore. Finally, V c i in Eq. ( 4) is evaluated based on transported species only (k ∈ M):

V c i = - M V k,i Y k M Y k .
This strategy has been successfully applied to the computation of 1D laminar methane/air premixed flames with the REGATH numerical code [17,18] in which transport equations are solved for a constant pressure with a Newton algorithm. The three solvers FTC, CTC and HTTC are presently considered within the finite volume solver DNS/LES SiTCom-B [14,19,20] that explicitly solves the unsteady fully compressible and reactive Navier-Stokes equations on cartesian meshes. From a numerical point of view, solving Eqs. (4) in addition to the momentum and energy equations in SiTCom-B requires to evaluate the density at each time step based on the contribution of major and minor species. However, the mass fractions of minor species are only accessible through the knowledge of the table parameters (Y c , Z) which are unknown at each time step n. As a consequence,

ρ n = M ρ n k + m ρ n-1 k (6)
with ( m ρ k ) n-1 ( M ρ k ) n for all n [10]. The error on ρ n is estimated to be very small since the variations of the quantities from two consecutive time steps should be small with a compressible code. Once the density is known, the mass fraction of major species can be determined by Y n k = ρ n k /ρ n , k ∈ M as well as the table parameters: Y n c and Z n . Minor species are given by Eq. ( 5) and any ωk , k ∈ M can be properly computed. The implementation of the three solvers into SiTCom-B has been validated by simulating a stoichiometric 1D laminar premixed methane/air flame with a variable Lewis number and without NO x chemistry. The kinetic scheme of Lindstedt and co-workers is used [21] (

N S = 29: H, OH, O, HO 2 , H 2 , H 2 O, O 2 , CO, CO 2 , N 2 , CH, HCO, CH 2 (S), CH 2 , CH 2 O, CH 3 , CH 3 O, CH 2 OH, CH 4 , C 2 H, HCCO, C 2 H 2 , CH 2 CO, C 2 H 3 , C 2 H 4 , C 2 H 5 , C 2 H 6 ,
C, and C 2 ; N R = 141) as in [10]. A comparison with the solution given by REGATH is provided in Fig. 1. With HTTC, the following criterion is used to split the N S species in major and minor species: the reactants and products having a non-zero mass fraction,

∆P Th = ρ u S 2 L (1 -ρ u /ρ b ) = -1.
03 Pa with the flame speed S L = 37.42 cm/s and the unburnt and burnt density set to ρ u = 1.130 kg/m 3 and ρ b = 0.1498 kg/m 3 , respectively. To ensure a stable temporal integration, the maximum time step used by the three solvers are 9.9 × 10 -9 s for FTC and 3.1 × 10 -8 s for HTTC and CTC. For the last two, the computational cost is reduced by a factor 3 thanks to the increase of the global time step, the cost per time step being roughly 10% less with HTTC than with FTC using Lindstedt's mechanism for methane combustion. Indeed, tabulating the minor species with HTTC allows for a strong increase of the chemical time step, which is the bottleneck of fully explicit compressible FTC solvers. In the present simulation, the data mining for CTC and HTTC has no impact on the CPU time.

Simulation setup and HTTC generalization

Simulation setup

The simulation setup is a pure methane slot injector surrounded by a co-flow of air, so that a steady laminar edge flame can stabilize above the burner. The width of the slot is D f = 2 mm, and the thickness of the injector wall is 0.5 mm. The simulated area along with the boundary conditions is given in Fig. 2 and 18 millimeters high, beginning at the injector outlet, the goal being to simulate the tip of the edge flame.

T o = 300 K U y,o (x) = U 0 y,o = 0.4 m/s Methane injector NSCBC Inlet U y,f (x) = U 0 y,f " 1 ✓ 2x D f ◆ 2 # T f = 300 K U 0 y,f = 6.0 m/s NSCBC Outlet P = 1 bar U y = U 0 y,o
The domain is uniformly meshed with 50 micrometer cells. Simulations are performed with the three solvers: FTC, CTC and HTTC. The detailed kinetic mechanism of Lindstedt [21] without NOx has been used as in [10].

For FTC and HTTC, species are transported with variable Lewis numbers. The tabulated species mass fractions are stored using uniform Y c and Z meshes, with small discretization step ∆Y c = 5×10 -4 and ∆Z = 0.001, to ensure that the parts of the table with strong derivative of the tabulated species mass fraction with respect to Z and Y c are described with a sufficient accuracy. A change of table size would have a limited impact on CPU time because the code used to access the table has been carefully optimized. Unsteady simulations are run until a steady state of the laminar flow is reached, i.e. when the stabilisation height of the flame tip and the maximum temperature in the domain become constant. The flame tip (located at a height y 0 ) is defined as the intersection of the stoichiometric mixture fraction (Z st ) isoline and an isoline of a small value of the progress variable (Y c = 0.005).

Generalization of HTTC

The flame databases required in CTC and HTTC formulations have been initially generated for a range of equivalence ratio between the flammability limits of CH 4 /air premixed flames, i.e. for φ ∈ [0.6, 1.4]. For cases where chemical reactions take place on a wider range of φ, such as in triple or edge flames, look-up tables need to be completed [15,22,23] or prolonged [14,24]. The latter strategy is presently used for the CTC method. It consists of an extension of all the tabulated variables (temperature, transport coefficients, etc.) out of the flammability range by a linear interpolation between the fresh and equilibrium states. The chemistry is not extended ( ωc = 0) since the chemical activity or S L is low out of the flammability limits. However, such a method cannot be directly transposed for HTTC. Indeed, the tabulated minor species mass fractions are zero both in the fresh gases and in the burnt gases, at equilibrium. Then, for equivalence ratios out of the flammability limits, their interpolated values would be 0. For many minor species this would lead to a gross error since they may take quite large values at the flammability limit. In figure 3, the mass fraction of C and C 2 H 2 in methane/air flames, extracted from the 1D premixed database are plotted in a (Y c , φ) space. It appears that both these minor species are not equal to zero at the rich boundary (φ = 1.4). At the lean boundary (φ = 0.6), some radical species such as C 2 H 5 , do not fade to zero either (not shown). The proposed extension method consists in generating additional flames out of the flammability limits, using the 1D flame solver RE-GATH. Even if it has been shown in experiments that premixed flames cannot propagate beyond these limits, REGATH is still able to provide data to extend the table. Accordingly, the range of φ has been augmented from [0.6, 1.4] to [0.5, 5] even if the results given by the 1D flame solver cannot strictly be considered as premixed propagating flames for very lean or very rich mixtures. Indeed for φ > 2.5, the flame front propagates at velocities smaller than 1 cm/s and is very wide. The 10 meter-long domains used to generate the database are not even sufficient to reach equilibrium values for all minor species. The extended HTTC table is also displayed in Fig. 3. For most of the tabulated species, their mass fraction reaches zero for values of the progress variable smaller than the equilibrium value. Moreover, for φ > 2, these species have a negligible mass fraction, whatever the value of Y c , so the part of the flame database for high values of φ is not used. This is, for example, the case for C. Such an extension of the table will then be a satisfactory approach for such species. Nevertheless, a few radical species are still present at high equivalence ratios with a non-negligible mass fraction (e.g. C 2 H 2 ). For these high values of φ, they do not even reach equilibrium at the end of the computational domain, so a part of the table is missing for these species. Besides, their mass fractions at the rich equivalence ratio boundary of the table are still not zero. These species will then be transported instead of being tabulated. Consequently, when using the extended table, 7 species (C 2 H 2 , C 2 H 3 , C 2 H 4 , C 2 H 5 , C 2 H 6 , CH 2 CO, CH 3 ) are added to the set M of the transported species (N M = 20). However, this additional cost does not affect the maximum time step used in simulations, the limiting species being still tabulated, and does not change the code efficiency.

Edge flame simulations

Flame structure with FTC

The field of heat release rate (HRR) is displayed in Fig. 4(a). A lean and a rich premixed zone exist along with a trailing diffusion flame but the premixed wings are merged with the diffusion flame exhibiting a monobrachial structure, as in [11], because of the large values of ∇Z (60 m -1 < ∇Z < 120 m -1 on the isoline Z st ) and the low value of the radius of curvature of the flame, R c , compared to the thermal flame thickness, δ T . A minimum value of R c = 0.5 mm is found on the isoline HRR = 0.05 GW/m 3 , significantly smaller than δ T ≈ 1 mm. This is consistent with the observations made by Kim et al. [25], where the curvatures of edgeflames is found close to R c = 1 mm in flows featuring gradients of mixture fraction around 50 m -1 on the isoline Z st . The maximum of HRR is located nearby the stoichiometric line, at the triple point, where the three parts of the flame merge. The temperature of the flow increases along the stoichiometric line downstream of the premixed front (Fig. 4 approximately 2 mm above the burner rim, at a radial location where the velocity of the flow is small enough to allow for the flame stabilization.

The contours of HRR plotted in Fig. 4(a) show that the premixed flame front is inclined with respect to the axial direction of the flow. Kim et al. [25] have attributed this phenomenon to the effect of the velocity gradient, which is usually strong in jet flames. As a consequence, the propagation velocity of the flame tip cannot be equal to the axial velocity of the flow at the triple point. Following [25], the propagation velocity is assumed to be equal to the velocity U n normal to the premixed front, which is presently calculated from the velocity vector U at the flame tip:

U n = U • ∇Y c / |∇Y c |.
U n = 0.209 m/s while the 1D laminar flame speed is S L (φ = 1) = 0.367 m/s. This smaller propagation speed is due to the strong gradient of mixture fraction at the flame tip. Finally, the reactive points of the simulated domain that contribute to the HRR (HRR > 1% of max (HRR)) are plotted in a (Y c , Z) space in Fig. 5. The flammability limits are here defined as the range of φ where reactive points are found (as in [26]). Even if this definition of the flammability range is different of what is used in experimental studies, it is clear that the limits are extended in the premixed zone, for φ ranging from 0.15 to 2.73. However, the reactive points that contribute to the diffusion flame are mainly found between φ = 0.6 and φ = 1.4, and for large values of Y c , because they are located in the burnt gases downstream of the premixed flame front. The maximum of HRR is found in the premixed area of the flame, but HRR is still significant in the diffusion trailing flame. Finally, reactive points are found beyond the equilibrium line, which is a consequence, the profiles of temperature taken at two heights above the flame tip differ from the reference case (FTC), even if the maximum values of the temperature are close (Fig. 7). Enlarging the range of φ as for the HTTC procedure has been tested without any further improvement on the propagation speed. To alleviate this issue inherent to tabulation methods, Nguyen et al. [START_REF] Nguyen | Combust[END_REF] proposed a multi-dimensional tabulation approach that takes the fluxes in the Z direction into account during the generation of the manifold. However, a five dimensions look-up table is needed, complicating the CTC procedure and increasing its numerical cost. height y 0 + 1 mm or diffusion at height y 0 + 5 mm and y 0 + 10 mm. The differences between the two simulations come from the tabulated values for radicals in HTTC that are taken from premixed flamelets, which are slightly different from what is observed in FTC. Thus, some small discrepancies appear in the physical space as shown in Fig. 8 for C 2 H 5 species. Despite these imperfections, the mass fraction field of major species the temperature and the propagation velocity of the flame tip are well reproduced by HTTC. With HTTC, the chemical time step is still larger than the convective time step, and simulations are performed approximately 5 times faster than with FTC. Actually, simulations with HTTC were first realized before switching to FTC ones in order to speedup the flame establishment and convergence.

Even though the HTTC simulations lead to satisfactory results by making use of tables based on 1D premixed flamelets, an even better agreement could be reached if diffusion effects were included in S2FT tables.

Conclusion

A CH 4 /air edge flame with a strong mixture fraction gradient is simulated with three solvers featuring different models for detailed chemistry. The fully transported chemistry (FTC) is the reference case for the tabulated complex thermochemistry (CTC) and the hybrid chemistry (HTTC) approaches. Conventional techniques such as CTC modeling fail to predict such complex flame topology because only the structure of premixed flamelets are tabulated without strain-rate-induced effects. However, it has been shown that the HTTC solver is capable of better capturing the shape, the propagation speed and the height of flame stabilization, even if the tabulated profiles are extracted from laminar premixed flames. With HTTC, the time-step is increased compared with detailed chemistry, and reaches values usually encountered for non-reactive flows with compressible solvers. Thus, the edge flame configuration is simulated with a computational cost divided by about 5. Future works should focus on the table generation procedure that could integrate constrained chemical equilibrium, on a direct comparison with a reduced chemistry on a challenging case and finally on the application of HTTC in the framework of Large-Eddy Simulation.

Figure 1 :

 1 Figure 1: One-dimensional atmospheric CH 4 /air flame at φ = 1. Mesh resolution in the flame front: ∆x = 20 µm. Squares: REGATH. Lines: SiTCom-B (black solid: FTC, blue solid: HTTC, green solid: CTC).

Figure 2 :

 2 Figure 2: Dimensions and boundary conditions of the simulation domain.

Figure 3 :

 3 Figure 3: Species mass fraction of radical species C (left) and C 2 H 2 (right) in the original 1D methane flame database (φ ∈ [0.5, 5.0]), as a function of Y c and φ.

Figure 4 :

 4 Figure 4: Edge flame simulated with the FTC solver. In green: isolines of Z corresponding to φ = 0.6, 1.0, 1.4 and 5.0. In pink: isolines of heat release rate (0.05, 0.2, 0.5, 1.0, 1.5 and 2.0 GW/m 3 ).

Figure 5 :

 5 Figure 5: Scatter plot of HRR (only for HRR > 1% of max (HRR)) for FTC simulations. (a) premixed flame, |FI| = 1, (b) diffusion flame, |FI| = 0. The red line is the equilibrium extracted from the 1D database. FI: flame index [26].

4. 2 .Figure 6 :Figure 7 :

 267 Figure 6: Comparison between CTC/HTTC and FTC simulations. Isolines of Z st (green) and Y c = 0.005 (black) for FTC (solid lines) and CTC/HTTC (dashed lines). HRR in red.

Figure 8 :Figure 9 :

 89 Figure 8: Mass fraction of C 2 H 5 (zoom) computed Z and Y c = Y CO + Y CO 2 . Green: isolines of φ. Black: isoline of Y c = 0.005.

  . It consists in a two-dimensional domain, 15 millimeters wide

	Symmetry Axis	NSCBC Side
		Air coflow
		NSCBC Inlet

Acknowledgments

Bastien Duboc was financially supported by the Region Haute-Normandie. Computations were performed using HPC resources from CRIANN. The preprocessing code to build up the HTTC table is available upon request at https://www.coria-cfd. fr/index.php/HTTC.