K Izumi

A.-M Lézine

Pollen-based biome reconstructions over the past 18,000 years and atmospheric CO 2 1 impacts on vegetation in equatorial mountains of Africa 2

This paper presents a quantitative vegetation reconstruction, based on a biomization 11 procedure, of two mountain sites in the west (Bambili; 5°56′ N, 10°14′ E, 2273 m) and east 12 (Rusaka; 3°26′ S, 29°37′ E, 2070 m) Congo basin in equatorial Africa during the last 18,000 13 years. These reconstructions clarify the response of vegetation to changes in climate, 14 atmospheric pressure, and CO 2 concentrations. Two major events characterize the biome changes 15 at both sites: the post-glacial development of all forest biomes ca. 14,500 years ago and their 16 rapid collapse during the last millennium. The rates of forest development between the biomes 17 are different; a progressive expansion of lowland biomes and an abrupt expansion of montane 18 biomes. The trends of pollen diagrams and biome affinity scores are not always consistent in 19 some periods such as the Younger Dryas interval and end of the Holocene Humid Period, 20 because the biomization method is not a simple summarization of the pollen data, but also takes 21 biodiversity into consideration. 22

Our sensitivity experiment and inverse-vegetation modeling approach show that changes 23 in atmospheric CO 2 concentration unequally influence vegetation in different local 24 environments. The study also suggests that the biome changes prior to the Holocene result from 25 both changes in the atmospheric CO 2 concentration and climate. The development of warm-26 mixed forest from xerophtic vegetation results from increases in atmospheric CO 2 concentration 27 and near-surface air temperature. Difference in local dryness results in the different biome 28 distributions, with more forest-type biomes at Bambili and more grass/shrub-type biomes at 29 Rusaka. 30 * Pollen nomenclature follows Vincens et al. (2007). 133 134

) are distributed along a decreasing rainfall gradient 135 from the Guineo-Congolian forests, which ranges from tropical rain forest (TRFO) and tropical 136 seasonal forest (TSFO) (near the Equator), to the Sudanian (to the North)/Zambezian (to the East

Introduction

Equatorial forests have been considered to be the most stable ecosystems on Earth. From 32 the viewpoint of their great floral and faunal diversity, they must have been in existence over the 33 last several million years even though their extent could have fluctuated with climate and 34 atmospheric composition. However, some previous palaeoecological studies (e.g., Maley and 35 Brenac, 1998;Vincens et al., 1999;[START_REF] Runge | Holocene landscape history and palaeohydrology evidenced by stable carbon 689 isotope (d13C) analysis of alluvial sediments in the Mbari valley (5°N/23°E), Central 690 African Republic[END_REF]Lézine et al., 2013a, b;Desjardin et al., 2013) 36 show that the equatorial forest ecosystems in Africa have undergone drastic modifications 37 (floristic, structural and palaeogeographic) in response to climate changes since the last glacial 38 maximum (LGM, ca. 21,000 years ago, 21 ka). These modifications include the possible 39 fragmentation of equatorial forested communities [START_REF] Maley | The African rain forest: main characteristics of changes in vegetation and 654 climate from the upper Cretaceous to the Quaternary[END_REF], the expansion of species 40 restricted to high elevations today [START_REF] Dupont | Vegetation change in equatorial West Africa: 576 time-slices for the last 150 ka[END_REF] at the LGM and the collapse of the 41 forests around 4-2 ka (Vincens et al., 1999;[START_REF] Lézine | Postglacial Pollen Records of Africa[END_REF][START_REF] Marchant | Rapid environmental change in African and South 661 American tropics around 4000 years before present: a review[END_REF]42 Lézine et al., 2013a,b). Changes in plant distribution and abundance from the last glacial 43 onwards were also observed in the lowlands all over North Africa, from the Equator to the 44 Northern Tropic (e.g., [START_REF] Lézine | Are modern pollen data 649 representative of west African vegetation?[END_REF][START_REF] Watrin | Plant migration and ecosystems at the 726 time of the "green Sahara[END_REF][START_REF] Hély | Holocene changes in African vegetation: 602 tradeoff between climate and water availability[END_REF]. 45

Moreover, 13 C measurements on leaf waxes implied the replacement of tropical montane forest 46 by scrub vegetation, the downward migration of alpine treelines and the marked shift towards 47 C4-plant dominance in the tropics during the last glacial period (e.g., Street-Perrott et al., 1997;48 Huang et al., 1999). 49

Climate factors, such as moisture and heat are commonly invoked to explain the changes 50 in ecosystem composition and structure (e.g., [START_REF] Lézine | Sahara and Sahel 643 vulnerability to climate changes, lessons from paleohydrological data[END_REF][START_REF] Anadón | Climate change will increase savannas at the 536 expense of forests and treeless vegetation in tropical and subtropical Americas[END_REF]. Other 51 factors, such as atmospheric pressure and atmospheric CO 2 could potentially also have some 52 impacts on the vegetation through physiological processes. Reduced partial pressure of CO 2 and 53 O 2 associated with an increase in altitude related to lower sea level during glacial times could 54 influence photosynthesis (e.g. [START_REF] Friend | Evolutionary and ecophysiological responses of mountain 586 plants to the growing season environment[END_REF][START_REF] Terashima | Is photosynthesis suppressed at higher 701 elevation due to low CO 2 pressure[END_REF]Sakata and 55 Yokoi, 2002). A decrease in atmospheric CO 2 generally results in a reduction of in the 56 abundance of plants with C3-photosynthesis pathway because of the required increased rates of 57 photorespiration, and an expansion of C4-plants due to their adaptation anatomically and 58 physiologically to low atmospheric CO 2 concentrations (e.g., [START_REF] Ehleringer | C4 photosynthesis, atmospheric CO 2 and 579 climate[END_REF]Cowling 59 and Sykes, 1999). Simulations with the BIOME3 equilibrium vegetation model also indicated 60 that low atmospheric CO 2 alone could result in the observed replacement of tropical montane 61 forest by scrub vegetation at the LGM [START_REF] Jolly | Effect of low glacial atmospheric CO 2 on tropical African 617 montane vegetation[END_REF]. To understand vegetation 62 changes at tropical high-elevation sites in paleoecological context, we therefore need to 63 understand the effect of changes in the atmospheric CO 2 concentration and air pressure to the 64 vegetation as well as purely climatic effects. 65

In this study we focus on changes in biomes at roughly similar altitude (2000~2300 m) in 66 tropical Africa: to the West in the Cameroon volcanic line (Bambili) and to the East in the Kivu 67 montane range (Rusaka). Detailed pollen studies at both sites (Lézine et al., 2013a;Bonnefille et 68 al., 1995) have shown that vegetation composition varied considerably through time from 18 ka 69 to the present. The goals of our paper are 1) to discuss biome changes in these mountain areas 70 over the past 18 ka and 2) to investigate the impacts of changes in atmospheric CO 2 71 concentration on vegetation at the early in the last deglaciation (i.e., 18 ka) using an inverse 72 vegetation-modeling approach. The comparison between two distant sites will enable us to 73 identify the more prominent climate-change events that have affected the Afromontane forests. where the Afromontane forest is preserved today [START_REF] Letouzey | Etude phytogéographique du Cameroun[END_REF][START_REF] Letouzey | Notice de la carte phytogéographique du Cameroun au 1 : 500 000[END_REF]Momo Solefak, 2009). 83

Regional precipitation, with the rainy season from March to October, is due to the West African 84 monsoon. The temperature is lower relative to lowlands over the tropical Africa due to the 85 altitude of the site. 86 87

Figure 1. Location of the two pollen sites: Bambili (Lézine et al., 2013a) and Rusaka (Bonnefille 88 et al., 1995); Guineo-Congolian forests (blue): rain forest (dark blue), deciduous forests (medium 89 blue) and mosaic of forest and savanna (light blue). Sudanian, Zambezian savannas, Sahelian and 90

Somalia-Masai steppes and deserts (yellow) [START_REF] White | The vegetation of Africa[END_REF]. Arrows indicate the direction of main 91 wind flow at 925 hPa in winter (left) and summer (right). 92

The Burundi highlands are a part of the Albertine Rift Mountains that enclose the western 93 branch of the East African Rift, following a roughly North-South direction. The mountain ranges 94 include high mountains, such as the Virunga Mountains (4507 m) and the Rwenzori Mountains 95 (5109 m). These altitudes are not reached in Burundi where the highest peak reaches 2684 m 96 only. Rusaka is a swamp lying at 3°26' S, 29°37' E and 2070 m in altitude (Figure 1). The regional climate of Rusaka is related to the South African monsoon, the rainy season is from 98 November to May. 99

The African biomes 100

Fossil pollen data are generally expressed in the form of abundances of individual plant 101 taxa, and detailed pollen descriptions at Bombili and Rusaka have been provided elsewhere 102 [START_REF] Bonnefille | 546 Glacial/Interglacial record from Intertropical Africa, high resolution pollen and carbon data 547 at Rusaka, Burundi[END_REF]Lezine et al., 2013). The pollen sequence is continuous and extends from 103 0 to 18,071 years ago at Bambili and extends from 750 to 18,061 years ago at Rusaka. Here we 104 use biomes, which are geographically and climatically broadly distributed physiognomic 105 vegetation types, for representing vegetation changes in this study. The biomes are represented 106 by assemblages of plant functional types (PFTs) that are defined on the basis of plant traits (e.g., 107 the life form, leaf form, phenology, and bioclimatic tolerances) that reflect their preferable 108 environments, in which the species maximize productivity and minimize environmental stress 109 (Table 1). The use of biomes and PFTs helps to solve the problem of classifying paleoecological 110 records by reducing the number of entities considered and by proving an ecological basis for 111 treating plants from different regions in a comparable way. A method for converting pollen taxa 112 to biomes (i.e. biomization) is described in section 3.1. 115

The Afromontane vegetation of Africa is discontinuous, with patches separated from one 116 other by lowlands, and thus is referred to as the "Afromontane archipelago" [START_REF] White | The vegetation of Africa[END_REF]. 117

Despite the geographic discontinuity, they share numerous plant species that are distinct from the 118 surrounding lowland regions. Three main biomes (Table 2), common to all Afromontane regions, 119 are distinguished. They correspond to an elevation gradient from roughly 1600 m to the top of 120 the highlands: 121 --Warm mixed forest (WAMF) occurs in a lower ombrophilous areas and in the lowland 122 Guineo-Congolian forests; 123 --Afro-alpine forest (AAF) occurs in an areas at higher elevation typically above 2000 m, 124

with upper limit of the forest is typified by the presence of abundant Ericaceae; and 125 --Afro-alpine grassland (AAG) is cool afro-alpine grasslands that consist of C3 grasses 126 which are found neR the top of the mountains above 2800 m. 127

Regional differences between the eastern and western mountain ranges are observed, however, 128 with e.g., Hagenia, Clifforita, Afrocrania and Junioerus absent from the western sector 129 (Cameroon), as well as Artemisia that has never been collected here. and South) tropical dry forests (TDFO), to a mixture of woodlands and grasslands (i.e. savanna 138 (SAVA)), to the Sahelian (to the North)/Somalia Masai (to the East) steppes (STEP), and finally 139 to desert (DESE). 140

In this study, we focus on nine biomes (three biomes for highlands and six biomes for 141 lowlands) on the basis of our own field expertise and local descriptive botanical literature 142 [START_REF] Troupin | Flore des plantes ligneuses du Rwanda[END_REF][START_REF] Letouzey | Etude phytogéographique du Cameroun[END_REF][START_REF] Letouzey | Notice de la carte phytogéographique du Cameroun au 1 : 500 000[END_REF]Momo Solefak, 2009). 143

Methods 144

Biomization procedure 145

Biomization is a quantitative procedure that reconstructs biomes on the basis of the 146 characteristic signature in the pollen record of different PFTs [START_REF] Prentice | BIOME 6000: reconstructing global mid-Holocene vegetation 687 patterns from palaeoecological records[END_REF]. There 147 are five steps in the biomization method used here: (1) assignment of taxa represented in the 148 pollen assemblages to PFTs (i.e., defining a taxon × PFT matrix), (2) definition of biomes as 149 combinations of PFTs (defining a PFT × biome matrix), (3) combination of the above two 150 matrices (defining a taxon × biome matrix by simple matrix multiplication of the matrices from 151 steps 1 and 2), (4) for a particular pollen assemblage, calculation of affinity scores for each 152 biome, and ( 5) selection of a biome with the highest affinity score as the dominant biome 153 represented by that pollen assemblage. The detailed formula for step 4 is described in Prentice 154 and Webb (1998). The affinity score for each biomes can be thought of as a measure of the likely 155 presence at a site of that biome given the particular pollen spectrum, such that the lower affinity 156 the score, the less likely the biome is to be present. The score is not equivalent to the proportion 157 of the area covered by an individual biome at the site, but simply describes presence or absence.

The biomization procedure has been successfully used worldwide to reconstruct modern 159 and past biomes for selected time periods, typically the last glacial period and the Holocene (e.g. 160 Prentice et al., 2000;Elenga et al., 2000;Jolly et al., 1998a). Within tropical Africa, the method 161 has been applied to modern pollen samples at the continental scale [START_REF] Jolly | Simulated climate and biomes of 614 Africa during the late Quaternary: comparison with pollen and lake status data[END_REF] or more 162 regionally in East [START_REF] Vincens | Modern pollen-based 712 biome reconstructions in East Africa expanded to southern Tanzania[END_REF] and West Africa [START_REF] Lézine | Are modern pollen data 649 representative of west African vegetation?[END_REF]. Past biome 163 reconstructions have also been performed for Plio-Pleistocene [START_REF] Bonnefille | High-resolution vegetation and 543 climate change associated with Pliocene Australopithecus afarensis[END_REF]Novello 164 et al., 2015) and more recent (from the last glacial to the present) sections [START_REF] Lebamba | Pollen, vegetation and climate at Lake Barombi Mbo 632 (Cameroon) during the last ca. 33 000 cal yr BP: a numerical approach[END_REF]165 Amaral et al., 2013). 166

We have used the complete matrices, pollen taxa-PFTs-biomes for both Bambili and 167 Rusaka defined by [START_REF] Vincens | Modern pollen-based 712 biome reconstructions in East Africa expanded to southern Tanzania[END_REF] and [START_REF] Lézine | Are modern pollen data 649 representative of west African vegetation?[END_REF], respectively, and the PFTs-168 biomes matrix (Table 1, 2) is the same for both sites. The lakeshore aquatic plant taxa and ferns 169 taxa were removed as these respond to local hydrological conditions, rather than being reflective 170 of broader scale climate controls. While a 0.5 % threshold for all taxa has generally been used in 171 calculating affinity scores to reduce the incidence of misassignment among relatively species-172 poor assemblages [START_REF] Prentice | BIOME 6000: reconstructing global mid-Holocene vegetation 687 patterns from palaeoecological records[END_REF], we selected 0.2 % threshold for this study after an 173 examination of several biomization practices. We also used a threshold (40%) for Poaceae in 174 order to minimize the over-representation of this taxon in individual pollen spectra. 175

A sensitivity experiment for a simple CO 2 effect on vegetation 176

In order to explore the impact of change in atmospheric CO 2 alone on vegetation at 177

Bambili and Rusaka, we first performed a sensitivity experiment. The experimental design is 178 similar to the one in [START_REF] Jolly | Effect of low glacial atmospheric CO 2 on tropical African 617 montane vegetation[END_REF]; the observed seasonal climate data is held 179 constant and only atmospheric CO 2 concentration is varied (from 400 ppm to 180 ppm, 10 ppm [START_REF] Izumi | Application of paleoenvironmental data for testing climate models and 608 understanding past and future climate variations[END_REF] and BIOME4 (Kaplan et al., 2003). These models simulate common equilibrium 182 vegetation and bioclimatic variables, but they have different vegetation responses to change in 183 atmospheric CO 2 concentration (e.g., BIOME5-beta has lower carbon use efficiency, the ratio of 184 net primary production to gross primary production). In running these vegetation models, we use 185 monthly climate from CRU CL 2.0 [START_REF] New | A high-resolution data set of surface climate 666 over global land areas[END_REF] and soil (FAO, 1995) at each site. 186

Altitude-adjusted air pressure is used at each site (7.7 × 10 4 Pa at Bambili and 7.9 × 10 4 Pa at 187 Rusaka). 188

Inverse-modeling through an iterative forward modeling approach 189

In order to examine the potential impacts of altitude and changes in atmospheric CO 2 190 concentration on vegetation from the viewpoint of pollen-based climate reconstruction, we used 191 an "inverse modeling through iterative forward modeling" (IMIFM) approach [START_REF] Izumi | Application of paleoenvironmental data for testing climate models and 608 understanding past and future climate variations[END_REF]192 Izumi and Bartlein, in revision), which can be compared with the forward modeling approach 193 that uses inputs of climate and CO 2 concentrations to mechanistically simulate vegetation. The 194 IMIFM approach (also called "inverse vegetation modeling" approach for climate reconstruction 195 in [START_REF] Guiot | Inverse vegetation 595 modeling by Monte Carlo sampling to reconstruct palaeoclimates under changed 596 precipitation seasonality and CO2 conditions: application to glacial climate in 597 Mediterranean region[END_REF] and Wu et al. (2007a)), was developed to overcome some disadvantages of 196 conventional statistical reconstruction approaches, such as modern-analogue, regression, and 197 response-surface techniques. These conventional approaches generally require several ecological 198 assumptions for climate reconstruction from pollen data, in particular that climate is the ultimate 199 cause of change in vegetation, and the modern data contain all the information necessary to 200 interpret the paleodata [START_REF] Guiot | A few prospective ideas on climate 599 reconstruction: from a statistical single proxy approach towards a multi-proxy and 600 dynamical approach[END_REF]. However, plant-climate interactions are very sensitive 201 to atmospheric CO 2 concentration (e.g., [START_REF] Cowling | Physiological significance of low atmospheric CO 2 for 568 plant-climate interactions[END_REF][START_REF] Prentice | Ecosystem effects of CO2 concentration: evidence from 679 past climates[END_REF] and thus modern pollen samples influenced by higher CO 2 concentrations of the past century are 203 not necessarily good analogs for climates under lower CO 2 concentrations. 204

The basic assumption of the IMIFM approach is that it should be possible to reconstruct 205 the climate data that gave rise to a "target" paleovegetation sample by searching for the set of 206 The IMIFM approach has the potential to provide more accurate quantitative climate 214 estimates from pollen records than statistical approaches because it allows the mechanistic 215 effects of non-climatic variables, such as the atmospheric CO 2 concentration and atmospheric 216 pressure, to be explicitly considered in the reconstruction. However, the approach is strongly 217 dependent on the quality of the forward vegetation model. Therefore, to reduce the dependency 218 of our results on a single model, we use two equilibrium vegetation models, BIOME5-beta and 219 BIOME4. The detailed methodology and vegetation models consulted for individual papers. 220

In the IMIFM approach for climate reconstruction, we need to compare pollen-based 221 observed biomes with the simulated biomes in the vegetation models. We define the simulated 222 biomes as follows: the afro-alpine forest (AAF) is composed of simulated scrubland and mixed 223 forest, and the afro-alpine grassland (AAG) is the simulated temperate grassland, but the climate spaces for both AAF and AAG are based on the climatic requirements of temperate 225 microphyllous shrub vegetation (i.e., the minimum mean temperature of the coldest month 226 (MTCO) is 5 °C, and the maximum of the MTCO is 10.5 °C) in Jolly et al. (1998b). The climate 227 spaces for simulated warm mixed forest (WAMF) are based on BIOME4's description of 228 WAMF, but the minimum of MTCO is increased to 10.5 °C. Below that threshold in MTCO the 229 simulated biome is regarded as AAF. In searching the climate scenarios over the deep tropical 230 areas, we allow for differences between mean temperature of the warmest month (MTWA) and 231

MTCO of less than 7 °C. 232

In order to investigate the effects of changes in atmospheric CO 2 concentration to climate 233 spaces of a target biome (i.e., AAF at Bambili and AAG at Rusaka) at the early in the last 234 deglaciation (i.e., 18 ka), we set two experiments using the IMIFM approach at each site; exp. 1) 235 paleo atmospheric CO 2 concentration (194 ppm for ca. 18 ka; [START_REF] Bazin | Carbon dioxide composite data on AICC2012 chronology[END_REF] and exp. 2) 236 modern atmospheric CO 2 concentration (331 ppm). The difference between exp.1 and exp. 2 237

shows the effects of atmospheric CO 2 difference to the target biome. In running these vegetation 238 models, we use monthly climate from CRU CL 2.0 [START_REF] New | A high-resolution data set of surface climate 666 over global land areas[END_REF] and soil (FAO, 1995) at 239 each site. Altitude-adjusted air pressure is used at each site (7.7 × 10 4 Pa at Bambili and 7.9 × 240 10 4 Pa at Rusaka). 241

Results 242

Vegetation changes 243

First, we present summary pollen diagrams showing percentages of the three main pollen 244 groups: trees, herbs and undifferentiated (which correspond to pollen grains determined at a too 245 low taxonomic level and/or pollen grains corresponding to plants with a variety of life forms) in order to illustrate the main physiognomical changes (i.e., ratio between trees and herbs) over the 247 last 18 ka at Bambili and Rusaka (Figure 2a-b). Both sites have varying degrees of vegetation 248 change, but the pollen percentages show some common trends between the two sites over the last 249 18 ka; a progressive expansion of forests over the last glacial-interglacial transition period, in 250 particular during the Bølling/Allerød warm period (ca. 13.8 ka) and the Holocene, and 251 degradation of the forests during the Younger Dryas interval (ca. 12.9 to 11.7 ka), and at ca. biome scores are a measure of the likely presence of the biome; both highland (afro-alpine forest 258 (AAF) and warm mixed forest (WAMF)) and lowland (tropical rain forest (TRFO), tropical 259 seasonal forest (TSFO) and tropical dry forest (TDFO)) forest biomes show relatively high 260 affinity scores through the entire period at Bambili (Figure 2g and2i). From 18 ka, they 261 increased, through fluctuating, to an optimum that they reach at ca. 14.5 (for the AAF and 262 WAMF) and 10.3 ka (for the TRFO, TSFO and TDFO) during the last deglaciation. The highest 263 affinity scores of all these forest biomes are reached during the phase of developing forests (ca. 264 11.5 ka) according to the pollen diagram (Figure 2a). During the Holocene, in contrast to the 265 relative stability of the montane forest biomes, affinity scores of lowland forests ones decreased 266 from ca. 10.3 ka onward. While they remained relatively stable in spite of a period of slight 267 decline from ca. 4.7 to ca. 1.15 ka, the 3.3ka event (the significant fall of tree pollen percentage) 268 is not recorded in our forest biome reconstructions. On the contrary, all the forest biomes 269 abruptly collapsed at ca. 1.15 ka as shown by the dramatic fall of their cumulative scores, of 50% 270 within ca. 70 years only. It was only after ca. 0.5 ka that forest biomes recovered and then 271 reached their Holocene affinity-score values (Figure 2g

and 2i). 272

The open grass/shrubland biomes (afro-alpine grassland (AAG), steppe (STEP) and 273 desert (DESE)) progressively diminished from ca. 18 ka to ca. 8.2 ka, and then gradually 274 increased until the present at Bambili (Figure 2c and2k). The higher affinity scores of some 275 xerophytic biomes (STEP and DESE) at ca. 3 and ca. 1.15 ka are also consistent with the herb 276 pollen percentage (Figure 2a and2c). The trends of these biome scores are similar to the trend of 277 herb pollen percentages. On the other hand, the affinity score of savanna (SAVA) was relatively 278 high throughout the whole period (Figure 2e).

The time interval between ca. 18 and ca. 15 ka had unfavorable conditions for tropical 280 lowland forests (TRFO, TRSO, and TDFO) and montane forests (AAF and WAMF) at Rusaka 281 (Figure 2h and2j). The tropical lowland forests started to develop at ca. 14.5ka and then 282 progressively increased until ca. 8.5 ka. After the forest expansion was abruptly interrupted 283 between ca. 8.5 and ca. 8 ka, the lowland forests remained stable until ca. 1.5 ka. On the other 284 hand, the montane forests abruptly increased at ca. 15 ka, then followed different trends: in 285 contrast to the AAF which remained stable, the WAMF continuously expanded until ca. 4 ka and 286 then remained relatively stable until ca. 1 ka. Among all the phases of forest decrease which 287 punctuated the African Humid Period (ca. 9 ka to ca. 6 ka), only that corresponding to the 8.2 ka 288 event is clearly reflected by the sharp decrease of he scores of all the lowland forest biomes. The 289 tropical lowland forests and WAMF recorded a dramatic decline during the last few hundred 290

years. 291

The affinity scores of xerophytic biomes (AAG, STEP, DESE and SAVA) are higher 292 than those for forest biomes during the period between ca. 18 ka and ca. 15 ka, which is 293 consistent with the pollen percentages at Rusaka (Figure 2b,2d,2f and 2l). The affinity scores of 294 the xerophytic biomes did not largely change through the entire period, but there is an opposite 295 trend between highlands and lowlands: the AAG progressively decreases and the other biomes 296 progressively increase from the glacial period to the present. The Holocene evolution of SAVA, 297 STEP and DESE was punctuated by phases of slight reduction at ca. 12, 10.5, 3.7-3.2, and 2.1 298

ka. 299

Figure 3 shows the first and second dominant biomes for representing the main biome 300 changes at Bambili and Rusaka sites. We chose a biome with the highest affinity score as the 301 first dominant biome and a biome with the next highest score as the second dominant one. The score express the potential co-existence of these biomes. The first dominant biome over the 304 Holocene is the warm-mixed forest (WAMF), which is in particular the sole dominant biome 305 except ca. 1 ka and ca. 9 ka, and it has relatively high affinity scores at both sites, but there are 306 different biome distributions and have lower affinity scores (than the Holocene) between the two 307 sites prior to the Holocene. The WAMF was already the first dominant biome from ca. 15 ka, 308 and before the establishment of the forest, savanna (SAVA) was the first dominant biome at 309 Bambili. On the other hand, the WAMF was first dominant biome from ca. 13 ka, and before the 310 forest establishment, afro-alpine forest (AAF, ca. 15 ka to ca. 13 ka) and afro-alpine grass (AAG, 311

Effects of atmospheric CO 2 concentration to vegetation 318

Unlike the sensitivity experiment in [START_REF] Jolly | Effect of low glacial atmospheric CO 2 on tropical African 617 montane vegetation[END_REF], our sensitivity 319 experiment does not show a shift of the warm mixed forest (WAMF) to the target xeric biomes (afro-alpine forest (AAF) at Bambili and afro-alpine grassland (AAG) at Rusaka) by the 321 decreases in atmospheric CO 2 alone (not shown). But, under the lower CO 2 concentrations, the 322 simulated net primary production (NPP) largely decreases for most of the forest-related plant 323 functional types (PFTs) and slightly decreases for shrub-and grass-related PFTs (Fig. 4). The 324 sensitivity analysis thus shows that the xerophytic biomes have relatively higher NPP and thus 325 more easily expand than forest biomes under lower atmospheric CO 2 concentrations. To a greater 326 or lesser extent, these responses are consistent among the BIOME models (i.e., BIOME3, 327 BIOME4, and BIOME5-beta). To illustrate the possible climate spaces for target vegetation (i.e., AAF at Bambili and 334 AAG at Rusaka) under the different atmospheric CO 2 concentration, we estimated mean annual temperature (MAT, °C) and mean annual precipitation (MAP, mm/year) using the IMIFM 336 approach with both BIOME5-beta and BIOME4 vegetation models (Fig. 5). This target 337 vegetation is the observed first dominant biome, which has the highest affinity score, at the early 338 in the last deglaciation (i.e., 18 ka) (the rightest part in Fig. 3). According to observed climate 339 data, CRU CL 2.0 [START_REF] New | A high-resolution data set of surface climate 666 over global land areas[END_REF], although MAT is similar at both sites, MAP at Rusaka is 340 less than half at Bambili (Figure 5), and thus Rusaka is likely to experience drier conditions than 341

Bambili for vegetation at both the present and past. 342 The difference of atmospheric CO 2 concentration (i.e., 331 ppm at 0 ka vs. 194 ppm at 353 ca. 18 ka) also influences plant productivity and climate spaces for the target vegetation (i.e. 354 AAF at Bambili and AAG at Rusaka) (Figure 5 [1] and [2]). The change in atmospheric CO 2 355 level (i.e., from 331 ppm to 194 ppm) decreases the simulated GPP, AR, and NPP by 30-45% at vegetation models at Rusaka; BIOME4 shows a greater response than the BIOME5-beta related 358 to the different carbon-use efficiency simulation between models. At Rusaka, BIOME5-beta 359 slightly increases in AR and decreases in both GPP and NPP, and there is different responses 360 between grass-type biomes and woody-type biomes in the single model. Compared to the 361 modern CO 2 experiment (Figure 5 [2]), the paleo CO 2 experiment (Figure 5 [1]) with both 362 models shows that the MAT drops by around 4 °C and MAP decreases about 150 mm/year at 363 Bambili, and that MAT drops about 0.5-1 °C and MAP increases about 200-300 mm/year at 364

Rusaka. The response of the forest-type biome (i.e. AAF at Bambili) to the change in 365 atmospheric CO 2 concentration is much larger than one of the grass-type biome (i.e. AAG at 366 Rusaka). 367 5. Discussion 368

Vegetation changes 369

Inspection of the formula for calculating affinity scores in the biomization procedure, 370

shows that the scores for each biome do not simply reflect the sum of pollen percentages that 371 constitutes each biome. If the percentage of the pollen sum that constitutes each biome is fixed, 372 greater diversity among taxa for a given biome will automatically yield a higher affinity score 373 [START_REF] Prentice | BIOME 6000: reconstructing global mid-Holocene vegetation 687 patterns from palaeoecological records[END_REF]. Therefore, the trends of pollen abundance and affinity scores will not 374 necessarily be consistent even if we use the exact same pollen data (Figure 2). However, the 375 general trend of forest expansion, forest establishment, and forest degradation described using 376 both pollen percentages and biome scores at Bambili and Rusaka is mutually consistent here. migration of the tree line in mountains (Figure 6a), and TRFO/STEP (i.e., the most humid biome 379 versus the one of the most arid biomes) for implying the lowland changes in biomes mainly 380 linked to moisture changes (Figure 6b). AAG/AAF displays a parallel trend at both sites with 381 high values at levels dated from 18 to 14.5ka. AAG/AAG values abruptly fall at 14.5 (Ruzaka) 382 or 13.7ka (Bambili) then remain stable during the Holocene until the last millennium where a 383 sharp increase is observed. Conversely TRFO/STEP values increase from 18 ka to 14.7 384 (Rusaka)-14ka (Bambili) with values substantially higher at Bambili compared to Rusaka. In this 385 latter site, they remain stable during the Holocene, whereas the dramatically increase during the 386 mid-and late Holocene at Bambili. 387 388 Figure 6. AAG/AAF ratio illustrating the descent of the tree line at Rusaka and Bambili (a) and 389 TRFO/STEP ratio illustrating the expansion of the lowland forests in relation with increased 390 humidity (b). 30 °N JJAS insolation from Berger and Loutre, 1991. 391 One of the remarkable results of our analysis is the characterization of the period prior to 392 the Holocene forest establishment. The first dominant biome over the Holocene is the warm 393 mixed forest (WAMF) at both Bambili and Rusaka, but the biome distributions were different 394 between the two sites prior to the Holocene (i.e., WAMF and savanna (SAVA) at Bambili and 395 AAF and afro-alpine grass (AAG) at Rusaka) (Figure 3). While the extreme reduction of all 396 kinds of forests is coeval with both the expansion of the lowland grasslands and shrublands and the downward shift of the Afroalpine grasslands (Figure 6a) during glacial times at Rusaka, such 398 a shrinkage was not observed at Bambili where all forest biomes were present at that time -even 399 with noticeable reduction of their population size -and where the downward shift of the 400 Afroapline belt was much less pronounced (Figure 6a) 401

Our result proposes that the eastern and western sectors of the Guineo-Congolian forest 402 domain had the following distinct environmental conditions: northwest Cameroon benefited from 403 conditions favorable to forest persistence at both low and high elevation in spite of dryness 404 related to Heinrich event 1 (H1, ca. 17 ka) [START_REF] Stager | Cooling cycles, Heinrich event 1, and the 695 desiccation of Lake Victoria[END_REF]. The forest domain was probably 405 not as fragmented as previously thought [START_REF] Maley | The African rain forest: main characteristics of changes in vegetation and 654 climate from the upper Cretaceous to the Quaternary[END_REF]. Lower atmospheric CO 2 concentration 406 and regional dryness allowed for xerophytic biomes to expand in the lowlands. For some 407 vegetation, the impact of lower CO 2 concentration is equivalent to that of increased dryness 408 [START_REF] Loehle | Predicting Pleistocene climate from vegetation in North America[END_REF]. By comparison, the environment was too dry in the Burundi highlands for the 409 persistence of any kind of forests. 410

Atmospheric CO 2 concentrations over the Holocene were still lower than modern, but the 411 impact on the vegetation is generally considered to be negligible for practical purposes. Here, the 412 expansion of xerophytic biomes mainly results from increased dryness. As shown by the 413 TRFO/STEP ratio (Figure 6b), the expansion of moist forest biomes in the lowlands at both sites 414 closely matches the insolation trend (Braconnot et al., 2007a(Braconnot et al., , 2007b) ) and related fluctuations in 415 monsoon rainfall since ca. 15 ka (e.g., deMenocal et al., 2000;[START_REF] Weldeab | Holocene African droughts relate to 729 eastern equatorial Atlantic cooling[END_REF]Gasse et al., 416 2000Gasse et al., 416 , 2008;;[START_REF] Lézine | Sahara and Sahel 643 vulnerability to climate changes, lessons from paleohydrological data[END_REF][START_REF] Tierney | Northern hemisphere controls on tropical southeast African climate during the past 704 60,000 years[END_REF][START_REF] Tierney | Abrupt shifts in Horn of Africa hydroclimate since the last 706 glacial maximum[END_REF]. Unlike 417 Rusaka where the driest biomes continued into the Holocene, Bambili recorded a drastic decline 418 of STEP and DESE, which corresponded to the onset of the "African Humid Period" i.e the 419 widespread expansion of tropical forests and woodlands in northern Africa (deMenocal et al., cannot however exclude the role of human populations who widely expanded at that time 444 (Lézine et al., 2013b) in the forest degradation. 445

Effects of atmospheric CO 2 concentration on vegetation 446

Our study shows that changes in atmospheric CO 2 concentration impact on plant 447 productivity and vegetation distribution and thus possibly influence climate spaces for the 448 vegetation in equatorial mountains of African. These results are consistent with previous studies. 449

However, although [START_REF] Jolly | Effect of low glacial atmospheric CO 2 on tropical African 617 montane vegetation[END_REF] already proposed that low atmospheric CO 2 alone 450 could result in the observed replacement of tropical montane forest by scrub vegetation, this 451 dramatic change apparently depends on the local climate and vegetation model used. As a result, 452 whether the vegetation prior to the Holocene is due to lower atmospheric CO 2 concentration 453 relative to present, changes in the climate, or their combination is an unresolved question (Wu et 454 al., 2007b). Sensitivity experiments can be used to assess the impacts of a change in atmospheric 455 CO 2 concentration on vegetation, but it is difficult to consider the impacts of climate changes on 456 vegetation at the same time because similar vegetation can exist under a wide climate space. 457 Thus, the inverse-vegetation modeling approach (e.g., Wu et al., 2007a, b;Izumi and Bartlein, in 458 revision) can be one option for investigating impacts of changes in the atmospheric CO 2 459 concentration and climate on vegetation. 460

The changes in atmospheric CO 2 concentration potentially result in large impacts on 461 vegetation and climate reconstruction based on the vegetation in our sites. Our MAT estimation 462 at ca. 18 ka shows the similar values at Bambili and Rusaka because the target biomes are similar 463 afromontane biomes at the two sites. When we take into account modern MAT at each site, we 464 can therefore suggest similar temperature decreases at ca. 18ka in equatorial mountains of Africa. vegetation changes from xerophytic biomes (afro-alpine forest (AAF) and afro-alpine grass 467 (AAG)) to warm mixed forest (WAMF) prior to the Holocene result from both the increase in 468 atmospheric CO 2 concentration and surface air temperature. 469

Our mean annual precipitation (MAP) estimation does not vary accompanying changes in 470 atmospheric CO 2 concentration in particular at Bambili, partly because of the impact of lower 471 CO 2 concentration on water-use efficiency (WUE; the ratio of rate of photosynthesis to the rate 472 of transpiration) in vegetation. If climate did not vary between the present and the past (at ca. 18 473 ka), the paleo vegetation in particular xerophytic biomes would still appear drier than present, 474 simply because of decreased WUE. WUE is sensitive to changes in atmospheric CO 2 through 475 effects on stomatal conductance, and consequently a decrease in WUE results from low CO 2 476 concentration (Lochle, 2007;[START_REF] Prentice | Ecosystem effects of CO2 concentration: evidence from 679 past climates[END_REF]. Moreover, in spite of the similar 477 afromontane biomes at the two site, the MAP reconstructions are opposite between the two sites 478 because of the different WUE responses between forest-type and grass-type biomes to the lower 479 CO 2 concentration. This MAP increase at Rusaka also implies that climate reconstruction at ca. 480 18 ka would appear too dry if the decrease in atmospheric CO 2 concentrations was not 481 considered. There are different responses in plant productivity between AAG (at Rusaka) and 482 AAF (at Bambili) to a change in atmospheric CO 2 concentration (e.g., AAF has larger 483 responses), and this has a large impact on climate reconstruction but one that depends on the 484 vegetation models and vegetation types. The different types of biomes at Bambili and Rusaka 485 possibly result from both the regionally different dryness and the lower atmospheric CO 2 486 concentration.

significant issues is that the results depend on the vegetation model used [START_REF] Guiot | A few prospective ideas on climate 599 reconstruction: from a statistical single proxy approach towards a multi-proxy and 600 dynamical approach[END_REF]. In 489 this study, we used the two equilibrium vegetation models, and there were some opposing 490 responses between the models, and thus we need to further evaluate the models if we wish to 491 produce robust results. Moreover, the inverse-vegetation modeling approach in Wu et al. (2007a, 492 b) assumed that last glacial period was drier than the present and thus their potentially prescribed 493 precipitation range was lopsided (i.e., -90% to 50% of modern values). But, if we do not have 494 any evidence for the potential climate space (e.g. from lake-status data), we should search over a 495 wider climate space (e.g., -90% to 90% of modern values) for the target vegetation. 496

Conclusion

497

Comparing two mountain sites located in different environments of equatorial Africa 498 sheds some new light on the response of plant formation (i.e. biomes) to climate changes: 499 (i) two major events characterize the changes biome distribution over the last 18 ka. Taking into 500 account uncertainties of the age models, their timing is remarkably similar east and west of the 501 Congo basin: the post-glacial development of all the forest biomes ca. 14.5 ka and their collapse 502 during the last millennium; (ii) contrary to the lowlands where forests biomes expanded 503 progressively, the montane forest development was abrupt, occurring within centuries. Mountain 504 biomes then remained remarkably stable throughout the Holocene contrasting with the repeated 505 fluctuations in the arboreal forest cover revealed by the tree pollen percentages. In particular, the 506 forest decline during the Younger Dryas and at the end of the Holocene Humid Period at ca. 3.3 507 ka is not or only slightly reflected in the forest biome scores; (iii) the amplitude of the collapse of 508 all the forest biome during the last millennium is remarkable and points to the major impact of 509 the Medieval warm period in the deep tropics.

Our sensitivity experiment and inverse-vegetation modeling approach show that 511 atmopsheric CO 2 concentration unequally impact on vegetation due to different local 512 environments such as climates at each site. This study also suggests that the biome changes prior 513 to the Holocene resulted from both changes in the atmopsheric CO 2 concentration and climate; 514 the development and establishment of warm mixed forest from the xerophytic vegetation results 515 from increases in atmopsheric CO 2 concentration and near-surface air temperature. The 516 difference of local dryness also influences the biome distribution between the two sites, more 517 forest-type biome at Bambili and more grass/shrub-type biome at Rusaka. Finally, our climate 518 reconstruction proposes that the post-glacial climate in equatorial Africa may have been more 519 mesic than previous studies suggest. 520

  extends from the Gulf of Guinea to the Southwest to the Adamawa plateau to the Northeast. 78 Mean altitude decreases from West (around 2000 m in the "Grassfield" region) to East (around 79 1000 m in the Adamawa plateaus), punctuated by high mountains, such as Mount Cameroon on 80 the coastline (4095 m) and Mount Oku in the Western plateaus (3011 m). Bambili is a crater lake 81 located in the Western plateaus (05°56′ N, 10°14′ E, 2273 m; Figure 1) close to Mount Oku 82

  climate scenarios, which input to a forward vegetation model yields a simulated vegetation that 207 resembles the vegetation represented by a target fossil-pollen sample. The application of the 208 IMIFM approach involves the generation of many thousands of candidate sets of individual 209 climate-variable values that are individually discarded or retained depending on their ability to 210 correctly generate the observed vegetation using a specific forward model. The retained climate-211 variable values, which allow a correct simulation of the target vegetation, are then statistically 212 summarized to provide the reconstructed or estimated values of the climate variables. 213

  pollen diagrams (trees, herbs and undifferentiated; a-b) and affinity scores of 255 each biome (c-l) at Bambili and Rusaka.

  ca. 18 ka to ca. 15 ka) was the first dominant biome at Rusaka. 312 313 Figure 3. The first and second dominant biomes and their affinity scores at Rusaka and Bambili. 314 The first dominant biome affinity score is shown with the absolute value (top). The affinity 315 scores of the second dominant biomes are shown with the ratio of the second biome affinity 316 score to the first biome one (bottom). 317

328 329Figure 4 .

 3284 Change in net primary production (NPP) in the sensitivity experiment. In running the 330 vegetation models (i.e., BIOME4 and BIOME5-beta), we use the altitude-adjusted air pressure 331 (7.7 × 10 4 Pa at Bambili, 7.9 × 10 4 Pa at Rusaka). The climate input comes from CRU CL 2.0 332 data at each site. 333

343

  Figure 5. Potential climate (mean annual temperature (MAT) and mean annual precipitation 344 (MAP)) spaces for target vegetation (i.e., afro-alpine forest (AAF) at Bambili and afro-alpne 345 grass (AAG) at Rusaka) at ca. 18 ka (18,071 years ago at Bambili and 18,061 years ago at 346 Rusaka) using the IMIFM approach with two vegetation models through two experiments: [1] 347 paleo CO 2 concentration (194. ppm) and [2] modern CO 2 concentration (331. ppm). Dots 348 indicate median values, bold vertical lines indicate interquartile intervals (25 th to 75 th percentile), 349 and thin vertical lines indicate 90 % interval (5 th to 95 th percentile) for probability distribution of 350 each climatic variable from BIOME5-beta (blue) and BIOME4 (red) respectively. The gray lines 351 indicate observed modern values from the CRU CL2.0 (New et al., 2002) data at each site. 352

113Table 1 .

 1 Plant functional types proposed for the west and east African areas 114

	wte1	Warm-temperate broad-and needle-leaved evergreen (higher elevation)
	wte2	Warm-temperate broad-and needle-leaved evergreen (lower elevation)
	c3a	C3 herb from lowlands areas at the boundary between the Saharan and the Mediterranean zones
	c3b	C3 herb from the afroalpine grasslands (top of the mountains >2800m alt)
	c4	C4 herb (tropical grasslands and savannas)
	PFT
	te	Tropical evergreen
	tr1	Wet tropical raingreen
	tr2	Dry tropical raingreen
	tr3	Driest tropical raingreen
	sf	Steppe forb/shrub
	df	Desert forb/shrub

Table 2 .

 2 West and east African biomes and their characteristic plant functional types (PFTs), 131 main phytogeographical affinities and main vegetation types[START_REF] White | The vegetation of Africa[END_REF]) 132

	Biomes	PFTs	main phytogeographical affinities	main vegetation types (main taxa *)
	AAG	c3b	Afromontane archipelago-like centre of	Afromontane udifferentiated montane vegetation (C3
			endemisms (VIII)	grasses mainly Poaceae, Alchemilla, Asteraceae)
	AAF	wte1	Afromontane archipelago-like centre of	Afromontane udifferentiated montane vegetation
			endemisms (VIII)	(Podocarpus, Prunus, Syzygium, Nuxia, Rapanea,
				Maesa, Myrsine, Hagenia, Juniperus, Hypericum,
				Ericaceae)
	WAMF	wte2	Afromontane archipelago-like centre of	Afromontane udifferentiated montane vegetation
			endemisms (VIII)	(Podocarpus, Ilex, Schefflera, Syzygium, Olea,
				Strombosia, Ficus, Celtis, Rutaceae, Sapindaceae,
				Euphorbiaceae)
	TRFO	te	Guineo-Congolian regional center of	Lowland rain forest: wetter types (Caesalpiniaceae,
			endemism (I)	Mimosaceae, Moraceae, Meliaceae, Irvingiaceae)
	TSFO	tr1	Guineo-Congolian regional centre of	Lowland rain forest: drier types / mosaic of lowland
			endemism (I)/	rain forest and secondary grassland / Swamp forest
			Guineo-Congolian/Sudanian/Zambezian	(Ulmaceae, Sterculiaceae)
			regional transition zone (XI)
	TDFO	tr2	Sudanian (III)/Zambezian (II) regional	Woodland / dry forest / Miombo (Brachystegia,
			centre of endemism	Julbernardia, Isoberlinia, Combretaceae, Lannea,
				Prosopis, Hallea, Monotes, Protea, Euphorbiaceae.)
	SAVA	tr3+c	Southern Sahel regional transition zone	Wooded grassland/ deciduous bushland/ mosaic of
		4	(XVI)/Guineo-	lowland rain forest and secondary grassland (C4 tall
			Congolian/Sudanian regional transition	(up to 3m high) grasses, Terminalia, Lophira,
			zone (XI) -Zambezian transition zone	Mitragyna inermis, Borassus aethiopium)
			(X)/Lake Victoria mosaic (XII)
	STEP	sf+c4 Northern Sahel regional transition zone	Semi-desert grassland and shrubland (C4 short
			(XVI)/Somalia-Masai regional centre of	grasses, Acacia, Commiphora, Balanites,
			endemism (IV)	Capparidaceae)
	DESE	df+c4	Sahara regional transition zone (XVII)	Desert (C4 short grasses,
		+c3a		Chenopodiaceae/Amaranthaceae, Resedaceae,
				Brassicaceae)

Equatorial mountains of Africa

2.1 Geographical features

interglacial variations in the treeline elevation in tropical Africa. Proceedings of the

National Academy of Sciences of the United States of America 104, 23, 9720-9724.

Acknowledgements

521

This research was funded by the National Research Funding Agency in France (ANR-09-PEXT-522 001 C3A) and the Belgian Federal Science Policy Office (BR/132/A1/AFRIFORD). Rusaka 523 pollen data were extracted from the African Pollen Database (http://fpd.sedoo.fr/fpd/) and 524

Bambili pollen data are available from the AML. We thank H. Chevillotte and S. Janicot for their 525 contribution. We are also grateful Patrick Bartlein and two anonymous reviewers for their 526 insightful comments on the manuscript. KI is supported by Labex L-IPSL and AML is supported 527 by CNRS. 528


Richard Pasquetti
email: richard.pasquetti@unice.fr

Comparison of some isoparametric mappings for curved triangular spectral elements

Keywords: spectral element method, isoparametric mappings, curved triangular elements, Fekete points

Using the spectral element method (SEM), or more generally hp-finite elements (hp-FEM), it is possible to solve with high accuracy various kinds of problems governed by partial differential equations (PDEs), see e.g. [START_REF] Komatitsch | Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles[END_REF][START_REF] Mercerat | Triangular spectral element simulation of two-dimensional elastic wave propagation using unstructured triangular grids[END_REF]. However, as soon as the physical domain is not polygonal, the accuracy quickly deteriorates if curved elements are not implemented. This is the reason why various methods have been developed during the last decades, starting from the celebrated transfinite interpolation proposed for quadrangular elements in [START_REF] Gordon | Construction of curvilinear co-ordinate systems and applications to mesh generation[END_REF]. In this note we revisit this problem for triangular elements, based on the use of Fekete points for interpolations and of Gauss points for quadratures, i.e. when using the so-called Fekete-Gauss approximation. As detailed in [START_REF] Pasquetti | Spectral element methods on unstructured meshes: comparisons and recent advances[END_REF], such an approach shows the so-called spectral accuracy. However, differently to the quadrangles based SEM, it does not involve diagonal mass matrices, see e.g. [START_REF] Cohen | Higher order triangular finite elements with mass lumping for the wave equation[END_REF][START_REF] Giraldo | A diagonal mass matrix triangular spectral element method based on cubature points[END_REF][START_REF] Mulder | New triangular mass-lumped finite elements of degree six for wave propagation[END_REF] and references herein for works trying to preserve this nice property that is especially useful when addressing evolution problems with an explicit time marching. In the frame of the Fekete-Gauss T SEM (T , for triangle), the present study clearly points out the importance of a good choice of the bending procedure by comparing different isoparametric mappings for the Poisson and Grad-Shafranov PDEs.

Let Ω be a two-dimensional domain and consider an affine finite element mesh T h of Ω composed of simplices. As well known, if Ω is not polygonal then the mesh T h is no-longer satisfactory if using a high order P N approximation, where N is the total degree of the polynomial approximation. Usually, one substitutes an isoparametric P N mapping to the affine P 1 one for the triangles that should approximate the boundary Γ of Ω. If T ∈ T h is such a triangle, we assume it to be in the simple case where one of the vertices, say A 1 , is in Ω whereas the two others, A 2 and A 3 , are on Γ, and that the segment A 2 A 3 is curved whereas A 1 A 2 and A 1 A 3 remain straight, see Fig. 1. This may not be sufficient to address some more complex situations where the whole mesh should be deformed and not only the boundary cells, see e.g. [START_REF] Ding | An efficient dynamic mesh generation method for complex multi-block structured grids[END_REF][START_REF] Løvgren | Global C 1 maps on general domains[END_REF][START_REF] Persson | Curved Mesh Generation and Mesh Refinement using Lagrangian Solid Mechanics[END_REF]. Here we keep in mind the situation where the required mapping is local and not global, so that its computational cost remains negligible.

In the frame of a nodal T SEM, everything is given once knowing a set of interpolation points and a set of quadrature points in the reference element T = {(r, s) : r ∈ (-1, 1), s ∈ (-1, -r)}, and the images of the interpolation points in each element T of the mesh. As first described in [START_REF] Pasquetti | Spectral element methods on unstructured meshes: comparisons and recent advances[END_REF], for the Fekete-Gauss approach matrices are set up to compute at the Gauss points of T the derivatives, with respect to r and s, of the polynomial interpolant u N of any scalar function u defined by its values at the Fekete points. Such matrices, say D r and D s , are rectangular and of size m p × n p , where m p ≥ n p is the number of quadrature points, generally chosen for exact integration of polynomials of degree m = 2N , and n p = (N + 1)(N + 2)/2 is the number of interpolation points. Then, since knowing the images of the interpolation points, one can compute the Jacobian matrix and determinant at the quadrature points without expressing explicitly the isoparametric mapping from T to T . Thus, the problem resumes to localize at best the interpolation points in the element T .

Four approaches are investigated hereafter:

-The bending procedure that we introduced in [START_REF] Lazar | Fekete-Gauss spectral elements for incompressible Navier-Stokes flows: The twodimensional case[END_REF]; -The transfinite interpolation discussed for the triangles in [START_REF] Perronnet | Triangle, tetrahedron, pentahedron transfinite interpolations. Application to the generation of C 0 or G 1continuous algebraic meshes[END_REF]; -The harmonic extension; -The linear elasticity approach, see e.g. [START_REF] Persson | Curved Mesh Generation and Mesh Refinement using Lagrangian Solid Mechanics[END_REF].

Before going into the details of these different approaches, we assume to have at hand a parametric description of the boundary of Ω, i.e. , Γ is defined by a vector function, say x Γ (t) = (x Γ (t), y Γ (t)), where t ∈ [t min , t max ]. Then, for all the four approaches it is assumed that the boundary nodes, i.e. the interpolation points of Γ, are at the intersection of the lines joining the inner vertex A 1 and the interpolation points of the straight line A 2 A 3 , i.e. the edge of the element, say T , provided by the piecewise linear P 1 mapping. This requires one to solve, in general using a simple numerical procedure, N -1 (number of nodes of an edge without the end points) generally non-linear equations. We prefer using such a geometrical approach rather than the parametrization itself, simply because the same curve may be obtained in very different manners. It remains to define n ′ p inner interpolation points, with

n ′ p = n p -3N = (N -1)(N -2)/2.
A The bending procedure: Assume that F is the interpolation point obtained with the affine P 1 mapping. Let G and G be the points at the intersections of A 1 F with the line A 2 A 3 and the boundary Γ, respectively. Then, the bending procedure [START_REF] Lazar | Fekete-Gauss spectral elements for incompressible Navier-Stokes flows: The twodimensional case[END_REF] consists of stating that F is homothetic to F by the homothety of center A 1 and of ratio A 1 G/A 1 G, so that:

A 1 F = A 1 G A 1 G A 1 F and d ≡ F F = ( A 1 G A 1 G -1)A 1 F
where d stands for displacement. One notices that for this bending procedure, only the point G of Γ influences the location of the interpolation point F . In order to define the inner interpolation points, one has to solve n ′ p equations.

The transfinite interpolation: For the triangle, the method makes use of the barycentric coordinates, say (λ 1 , λ 2 , λ 3 ). Recall that λ i = 0 stands for the edge opposite to the vertex i, λ i = Constant to a line parallel to this edge and that vertex i belongs to λ i = 1. General transfinite interpolation formula are given in [START_REF] Perronnet | Triangle, tetrahedron, pentahedron transfinite interpolations. Application to the generation of C 0 or G 1continuous algebraic meshes[END_REF] (different approaches are however possible, see e.g. [START_REF] Demkowicz | Computing with hp-adaptative finite elements[END_REF]). For the triangle, if the displacement d vanishes at the two edges A 1 A 2 and A 1 A 3 , again assuming that A 2 A 3 is the curved edge, one has:

d(λ 1 , λ 2 , λ 3 ) = λ 2 d(0, 1 -λ 3 , λ 3 ) + λ 3 d(0, λ 2 , 1 -λ 2 ) .
This means that two points of the boundary Γ influence the location of an inner interpolation point F . They are those at the intersections of Γ with the lines parallel to A 1 A 2 and A 1 A 3 and passing by the point F . One notices that the method requires solving 2n ′ p equations.

Harmonic extension:

The goal is here to define the inner interpolation points by solving, for each curved element T , the weak form of the Laplace Dirichlet problem: ∆d = 0 in T , d| ∂ T = g with g given on A 2 A 3 and g = 0 on the two other sides. To resolve this vector Laplace problem, that in fact yields two uncoupled scalar problems, one has to set up the differentiation matrices, say D x and D ỹ that allows one to compute derivatives in T . To this end, one makes use of the differentiation matrices D r and D s and applies the chain rule. Since the mapping from T to T is affine, the Jacobian matrix and determinant Jacobian are constant. For the same reason, one may solve directly for the interpolation point coordinates. Finally, one should invert a matrix of size n ′ p × n ′ p , in order to compute the coordinates of the inner interpolation points. The approach is thus rather simple and in practice not costly. In [START_REF] Løvgren | Global C 1 maps on general domains[END_REF], where global mappings are considered, it is however outlined that the harmonic extension may fail to define a mapping.

Linear elasticity: Here the curved domain is viewed as the deformation of triangles, that is governed by the equation of linear elasticity. Introducing the Lame coefficients λ and µ, the displacement field is governed by the Navier-Cauchy equation:

µ∆d + (λ + µ)∇(∇ • d) = 0 in T , d| ∂ T = g .
There is now a coupling between the components of the displacement field. Using the ingredients previously described for the harmonic extension, the weak form of this elasticity problem yields a matrix system of size 2n ′ p . The solution of this system provides the displacements of the inner interpolation points. Again, one can also compute directly their coordinates. Because the forcing term is zero, the mapping is here parametrized by the ratio µ/(λ + µ) = 1 -2η, where η ∈ (-1, 0.5) is the so called Poisson ratio.

Fig. 1 shows the distribution of the interpolation points in a triangle, as obtained with the different strategies for a polynomial approximation degree N = 6. Of course, only the inner nodes differ. In these examples, the curved boundary is defined parametrically by : x = t, t ∈ (-1, 1), and y = 0.3 sin(πt) (left), y = 0.3 cos(πt/2) (center) and y = -0.3 cos(πt/2) (right). In the rest of the paper the different isoparametric mappings are compared for the Poisson equation, -∆u = f , with exact solution u ex = cos(10x) cos(10y) so that f = 200 u ex , and, from an example given in [START_REF] Ltjens | The CHEASE code for toroidal MHD equilibria Computer[END_REF], for the Grad-Shafranov PDE

-∇ • ( 1 µx ∇u) = f , with u ex = - κ 2ρ 3 q 1 4 (x 2 -ρ 2 ) + 1 κ 2 x 2 y 2 -a 2 ρ 2 , f = 1 + κ 2 µκρ 3 q x .
In both cases the boundary Γ of the computational domain is defined by Γ = {(x, y) :

x 2 = ρ 2 + 2aρ g(t) , y = κa ρ x sin(t) , 0 ≤ t < 2π} , g(t) = cos(t) + 0.75 exp(-2(t -π) 2 )
and Dirichlet boundary conditions are used. Moreover, as in [START_REF] Heumann | A high order Galerkin method for computing geometric coefficients of axisymmetric magnetohydrodynamic equilibria[END_REF][START_REF] Lee | ECOM: A fast and accurate solver for toroidal axisymmetric MHD equilibria[END_REF] we take ρ = 1, a = 0.32, κ = 1.7 and q = 1, but with g(t) = cos(t) to obtain a non convex domain that includes different types of curved elements.

For the Grad-Shafranov problem, the 1/(µx) value at the Gauss points is computed by interpolation of (µx) from the Fekete to the Gauss points, so that it is exact for non deformed triangles. For the linear elasticity extension, as in [START_REF] Persson | Curved Mesh Generation and Mesh Refinement using Lagrangian Solid Mechanics[END_REF] we use η = 0.4, yielding µ/(λ + µ) = 0.2. Computations have also been carried out without curving the boundary triangles, then (i) assigning to the boundary points of the straight edge the corresponding values of u ex on the curved boundary or (ii) simply using u ex , just like if the computational domain was a polygon. Computations have been carried out for the three P 1 meshes shown in Fig. 2 and for the polynomial approximation degrees N ∈ {1, 2, 3, 6, 9}. With u N for the numerical solution and {x k } for the interpolation points, log-log plots of the l ∞ error, ǫ = max k |u N -u ex |(x k )/ max k |u ex |(x k ), with respect to the polynomial degree N are shown in Fig. 3. Of course, if the boundary triangles are not curved then the results are bad for the real domain and on the contrary excellent for the polygonal one. If curved elements are involved, the convergence curves look more algebraic than exponential, with however high convergence rates. One notices that for N ≤ 2 all methods show the same accuracy, since then there is no inner interpolation point. With respect to the polygonal case, the loss of accuracy may be attributed to a failure of the Gauss quadrature formula: Indeed, since the Jacobian matrix and determinant are no longer constant in each element, polynomials of degree greater than 2N are numerically integrated. We have checked that this was negligible, since using values m > 2N yields no valuable gain in accuracy, see Table 1. For curved elements the differentiation matrices are however no-longer a simple combination of the matrices D r and D s defined for the reference element, which implies a sensitive loss of accuracy. With respect to the standard SEM, i.e. based on quadrangular elements, the fact that the interpolation and quadrature points do not coincide certainly constitutes another source of inaccuracy, see e.g. the SEM results in [START_REF] Heumann | A high order Galerkin method for computing geometric coefficients of axisymmetric magnetohydrodynamic equilibria[END_REF]. From the results of this Note, the bending and the transfinite interpolation methods appear to be less satisfactory for elements of high degree, especially in the Grad-Shafranov case, since a slow-down in the decrease of the error is clearly observed. Our explanation is that isoparametric mappings based on PDEs are more consistent than those based on interpolations, in the sense that the locations of the inner nodes do not only result from some particular points of the actual boundary but from its polynomial interpolant, which indeed defines the boundary of the computational domain. One may conjecture that such a conclusion could extend to the tetrahedron and to T SEM approximations different from the Fekete-Gauss one. Table 1: Max norm of the error for the Laplace and Grad-Shafranov problems for m ≥ 2N , using the transfinite interpolation, the medium mesh and for N = 6.

Figure 1 :

 1 Figure 1: Curved triangles and location of the interpolation points for the different isoparametric mappings (N = 6).

Figure 2 :

 2 Figure 2: Computational domain and meshes: coarse (69 elts), medium (151 elts) and fine (293 elts).

Figure 3 :

 3 Figure 3: Max norm of the error vs the polynomial degree N for the Laplace (at left) and Grad-Shafranov (at right) problems, using the coarse (top), medium (middle) and fine (bottom) meshes.

  10 -4 4.0246 10 -4 4.0065 10 -4 4.0078 10 -4 Gr. Shaf. 2.4653 10 -5 2.6421 10 -5 2.6306 10 -5 2.6313 10 -5

Acknowledgments:

The P 1 FEM meshes have been generated with the free software "Triangle". Special thanks to my Colleague Holger Heumann and to Lucas Drescher for providing me with the opportunity to carry out the present study and for stimulating discussions.