Christophe Corona
email: christophe.corona@gmail.com

Edouard Jean-Louis
email: jean-louis.edouard@mmsh-univ.fr

Guibal Frédéric

Lambert Georges Noël

Py Vanessa

Guiot Joël

Thomas André

Last-millennium summer-temperature variations in Brianconnais (French Alps) based on a composite tree-ring larch chronology

Keywords: Briançonnais, dendrochronology, Larix decidua Mill., construction wood, composite multicentennial chronologies, climate reconstruction, Mots clés: Briançonnais, dendrochronologie, Larix decidua Mill., bois de construction, chronologies composites pluri-centenaires, reconstruction climatique

We present a tree-ring-based reconstruction of last millenium summer temperature from Briançonnais in the northern French Alps. The sample material is derived from multi-centennial living trees and also includes larch (Larix decidua Mill.) series from historical construction timbers. The Regional Curve Standardization method was applied to preserve interannual to multicentennial variations in this highelevation dataset. The proxies are inferred using a linear model and calibrated using the June to August mean temperatures from the last revised version of the HISTALP database spanning the 1760-2003 and adjusted to take into account warm bias before 1850. The transfer function has crossvalidated coefficients of correlation (r) of 0.48 (0.42), reductions of error (RE) of 0.21 (0.27) and coefficients of efficiency of 0.20 (0.16). The concordance of the Briançonnais temperature reconstruction with other Alpine records suggests that the transfer function has successfully reconstructed past summer temperature during the last two millennia. Evidence was found of a warm episode between and 900 and 1200AD, and a cold phase between 1500 and 1850AD. These events were possibly correlated to the so-called 'Mediaeval Warm Period' (MWP) and 'Little Ice Age' (LIA). model reconstructed a summer temperature decrease more than 2°C for the LIA compared with the temperature prevailing during the MWP.

Introduction

To date, numerous temperature reconstructions have been produced from local tree-ring chronologies and regional-scale network compilations. However, most of them span the last 300-500 years only. There are relatively few that extend back prior to 1000AD. This lack is very apparent in [START_REF] Jansen | Palaeoclimate[END_REF], which shows only 16 locations globally from which tree-ring data have been used in millennium long temperature reconstructions. It is related to the scarcity of wood material from the early part of the last millennium and it limits the assessment of climate variations during the putative Medieval Warm Period (MWP) (Lamb,1977). These difficulties may be overcome by combining data from living trees with those from dry dead and sub-fossil wood (Nicolussi and Schiebling, 2002) together with dendroarchaeological sources. This is done for the Briançonnais region (French Alps) using living trees and historical samples from elevations sites >1,500 m asl, where summer temperature is the main limiting factor for many Alpine tree species [START_REF] Frank | Temperature reconstructions and comparisons with instrumental data from a tree-ring network for the European Alps[END_REF]. The resulting composite chronologies cover the past 1250 years and allow the assessment of long-term temperature trends for the region especially during the MWP and the Little Ice Age (LIA, e.g. [START_REF] Grove | The little ice age[END_REF]. In Sect. 2.1, the compiled dataset from living trees and historical timber is introduced. Individual ring width measurements are standardized using the regional curve standardization method (RCS, e.g. [START_REF] Briffa | Increasing productivity of ''natural growth[END_REF], which may help to preserve long-term, millennial-length trends in resulting chronologies (Sect. 2.2). It is demonstrated that living trees and historical chronologies share a strong common signal though they are combined into a composite chronology. In Sect. 2.3, the climatic signal of the 1,244-year composite historical/living chronology is assessed by calibration with HISTALP instrumental temperature measurements. A discussion of the new record, and comparison with other alpine proxies and reconstructions are provided in Sect. 4.

Material, methods and results

The tree ring-data

The tree-ring widths (TRW) used in this study originate from the Briançonnais-Queyras-Embrunais region in the French Alps (44°45-45°00N, 6°30-6°45'E, Fig. 1). The sample material is derived from multi-centennial living trees (Fig. 2A) and also includes larch (Larix decidua Mill.) series from historical construction timbers [START_REF]Les très vieux arbres vivants et les arbres morts, témoins et vestiges des forêts du passé dans les Alpes du Sud : lecture dendrochronologique d'un patrimoine naturelle et humain[END_REF][START_REF] Bernardi | Données nouvelles sur l'histoire de la Chapelle Saint Hippolyte (Névache, Hautes-Alpes, France). Apport de la datation dendrochronologique[END_REF][START_REF] Thomas | Cernes d'arbres et chronologie holocène dans les Alpes françaises[END_REF], Fig. 2B). A total of 561 increment cores were collected, processed, and measured using standard dendrochronological techniques [START_REF] Fritts | Tree rings and climate[END_REF]. These cores include 185 samples from living or dead larch trees from 14 locations >2000 m, 336 samples from traditional buildings (mountain summer chalets and chapels) and 40 samples from an archaeological mining site (Freissinières-Grand Bois; Py, 2009). These historic larch timbers are located above 1900 m asl, excepted two chapels located at about 1600 m asl (Puy Saint Vincent).

From each living and dead tree, two to three cores were extracted with an increment borer at breast height on the cross-slope sides of the trunk. Sections were sampled on tumbledown or under restoration buildings in various parts of the wood structure: roof structure (beams, purlins, rafters), load-bearing walls, ceiling and floor joists. Individual series from a same tree (core and radius of a section) were first crossdated. Using both statistical and graphical programs, e.g. TSAP [START_REF] Rinn | TSAP Reference Manual[END_REF] and Dendron II [START_REF] Lambert | Dendrochronologie, Histoire et Archéologie, modélisation du temps. Le logiciel Dendron II et le projet Historic Oaks[END_REF], Fig. 3), the interannual ring-width variability within each mean-tree series was compared against all other series and against local or regional chronologies. Before the 14 th century, the millennium-long master chronology from the Merveilles valley (Mercantour, Serre 1979) was used as reference. This procedure is known as 'cross-dating' (a description of 'European style' cross-dating methodology is given in Schweingruber, 1988). In order to identify locally absent rings, eliminate measurement errors and ensure relative and absolute dating accuracy, the quality of the cross-dating was repeatedly checked with several statistical measurements (coefficient of correlation, coefficient of coincidence, Student's t test). Problematic samples, such as samples with missing rings or a distorted growth pattern, were dismissed from the data set. By including the dry dead and the samples from historical buildings, progressively older material could be cross-dated and, thus, the record of precisely dated tree rings could be extended back to 751AD while the outermost ring in the current sample of living trees was produced during the summer of 1995AD. However, the data are not continuous over the last 1250 years. A gap exists between 1304 and 1337AD.

First, series were combined on a site-by-site basis. Sample replication and chronology characteristics for each site and for both groups are detailed in table 1. The number of series per dataset varies considerably from 8 (Granon 3) to 323 (Névaches granges). The mean segment length (MSL) -the average number of rings per core or disc sample -ranges from 139 years (Puy-Saint-Vincent) to 458 (Echalp), with a mean of 270 years. These values are important since the segment lengths of the individual detrended series fundamentally limit the maximum wavelength of recoverable climatic information [START_REF] Cook | The 'segment length curse' in long treering chronology development for palaeoclimatic studies[END_REF]. Average growth rate of 0.82 mm/year varies between 0.61 mm/year (Oriol, Puy-Saint-André) and 2.57 mm (Pas d'Archail). Finally, the larch data were 'horizontally' split into two groups representing samples from living trees and historic construction wood, as these subsamples considerably differ in: (i) the temporal distribution of the constituent series with older series, before 1303AD sampled on historical buildings (Fig 4A); (ii) the average cambial age of the constituent tree rings for each year (Fig. 4B); (iii) their mean growth rates of 0.72 and 0.91 mm, (iv) and mean segment length of 290 and 178 years, respectively (Tab. 1). Hence, two standardization runs were calculated.

Standardization

The annual growth of a tree, manifest as tree-ring width (TRW), may in theory be attributed to a restricted number of different growth-controlling factors [START_REF] Cook | A conceptual linear aggregate model for tree rings[END_REF]. In this study, sampled trees live close to their climatological limit of distribution and tree growth is assumed to be forced by two dominant factors: the climate [START_REF] Tranquillini | Physiological ecology of the alpine timberline[END_REF] and the age/size of individual trees. Non-climatic environmental factors will certainly have been important for individual trees at different times, but these factors are regarded as haphazard events with random distribution over space and time. To preserve long-timescale climate variability in long tree-ring chronologies while removing most of the age related variance, RCS was applied (e.g. [START_REF] Briffa | Increasing productivity of ''natural growth[END_REF]. This is a so-called age-dependent composite detrending method, where (1) a given number of power transformed measurements are aligned by cambial age, (2) a single growth function [regional curve (RC), smoothed using a spline function of 10% the series length] is fit to the mean of all age-aligned measurements, and (3) the deviations of the measurements from this smoothed growth function are calculated as residuals (details in Esper et al. 2003). All detrended series are averaged to chronologies using the biweight robust mean [START_REF] Cook | Methods of dendrochronology: applications in the environmental science[END_REF]. The variance is stabilized to minimize the effect of changing sample size using a technique that considers the number of samples per year and the crosscorrelation between the single measurements [START_REF] Osborn | Adjusting variance for sample size in tree-ring chronologies and other regional-mean time series[END_REF]. The signal strength and confidence in the standardized chronologies is assessed by calculating R-BAR and expressed population signal (EPS) statistics in a 50-year window moved in 25-year steps over the total length of the series [START_REF] Wigley | On the average of value of correlated time series, with applications in dendroclimatology and hydrometeorology[END_REF]. The R-BAR statistic is a measure of the average inter-correlation of all overlapping series, while EPS denotes the percent common signal. EPS is related to RBAR and to the number of replicate series, with EPS values above 0.85 generally regarded as satisfactory [START_REF] Wigley | On the average of value of correlated time series, with applications in dendroclimatology and hydrometeorology[END_REF].

While estimated TRW growth trends are remarkably similar for both groups (Fig. 4 C,D), a slight but systematic level offset between the TRW growth trends is seen. The RC based on the historic series indicates generally higher width compared to the RC of the living tree series, with largest differences found during juvenile growth. Resulting RCS chronologies of the historic and living trees measurements cover the 751-1894AD (with a gap between 1304 and 1337AD) and 1338-1995AD periods (Fig. 5 A,B). Figure 5C, D shows that the R-BAR values for the historic (mean: 0.48) and living tree chronology (mean: 0.45) are quite similar and the confidence in both chronologies is high, with EPS values typically above 0.85. In the historic chronology, EPS is on average 0.96, but with critically low values in some periods, particularly before 901 AD and 1210-1304 AD. In the living trees chronology, EPS is <0.85 between 1338 and 1413 AD.

Finally, after standardization, both RCS chronologies correlate at 0.57 over the well-replicated 1413-1894AD period of overlap and are merged for further analysis. The composite chronology, average of the two larch RCS chronologies, weighted by the number of series, is shown in Fig. 5E. After truncation at EPS>0.85, it covers the period 900-1210AD and 1413-1995AD.

Tree-growth response to climate

Monthly-homogenised records of temperature and precipitation from the HISTALP database (Auer et al. 2007) were used in order to identify climatic signals in the composite chronology and to investigate the optimum season for a climate reconstruction. This dataset consists in point values, gridded at 1°× 1° (temperature) and 0.1° x 0.1° (precipitation), lat × long, representing deviations from the 1961-90 mean. It results from a dense network of 134 meteorological stations, extends back to 1760 and covers the Greater Alpine Region (GAR,(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(43)(44)(45)(46)(47)(48)(49). We used its latest version, adjusted to take into account warm bias in summer temperature related to the insufficient sheltering of thermometers before 1850 (Böhm et al. submitted). Monthly mean series are calculated from 4 HISTALP points values (44-46°N, 6-7°E) enclosing the Briançonnais. Correlation coefficients between these 4 points over the 1760-2007 common period range from 0.89 to 0.99, demonstrating a rather homogenous climate field over the study area.

Correlation analyses using previous-year April to -year October monthly data and the RCS composite chronology were carried out. They give an indication of the direction and relative strength of the climatic forcing. The results show that TRW chronology has a significant positive response to June to August temperatures (Fig. 6A). JJA seasonal mean reveals highest correlations of 0.39. Moving 31-yr correlation analysis considering June, July, August and JJA indicates variable relationships for the individual monthly mean temperatures (Fig. 6B). Weaker correlations are obtained before 1880, and for August until 1940. Results for JJA are temporally more stable after 1880 and persistently significant at p < 0.05. No correlation with previous-year temperatures is significant at p < 0.05 and correlations with precipitation are weak.

A linear regression model was thus calibrated between JJA temperatures and the RCS composite chronology for the period 1760-1995. The fidelity of the reconstruction equation was tested in a calibration / verification exercise [START_REF] Fritts | Tree rings and climate[END_REF][START_REF] Cook | Methods of dendrochronology: applications in the environmental science[END_REF] where the Histalp mean summer (JJA) temperature was split into two independent 117-year periods. A linear regression was fitted over the calibration period, and the derived coefficients applied to the tree ring data in the verification period (Fig. 7A,B). This scheme produces independent estimates of temperature that can be compared to the instrumental temperature data. The verification results were tested using the squared Pearson correlation (R 2 ), the reduction of error (RE) and the coefficient of efficiency (CE). The R 2 statistic, ranging from 0.0 to 1.0 is a measure of the proportion of variation explained by the regressor in the model. The RE and the CE statistics both range from minus infinity to +1.0 with RE > CE and where values greater than 0 give confidence to the model performance [START_REF] Fritts | Tree rings and climate[END_REF][START_REF] Briffa | Reconstructing summer temperatures in northern Fennoscandia back to A.D. 1700 using tree-ring data from Scots pine[END_REF]. Of these three statistics, CE is the most difficult to pass [START_REF] Cook | Spatial regression methods in dendroclimatology: a review and comparison of two techniques[END_REF].

When the composite chronology is used in a linear regression, the R2, RE, CE, calibration and verification statistics against JJA mean temperatures indicate reconstructive skill against instrumental temperature data (Tab. 2), say r=0.42 and 0.48, RE=0.21, 0.27 and CE=0.20, 0.16. These calibration and verification statistics demonstrate some useful information preserved by the model [START_REF] Cook | Spatial regression methods in dendroclimatology: a review and comparison of two techniques[END_REF]. The mean discrepancy between reconstructed and instrumental temperatures is -0.3•C (Fig. 7A) for the 1760-1995 period. Discrepancies are maximal at 1931 (+1.5•C) and 1909 (-1.29•C). When the curves are smoothed with a 10-year low-pass filter (Fig. 7B), r=0.62 and 0.71 for the calibration and verification periods, respectively. The maximal discrepancies occur around 1790, 1890 and 1970AD. Interestingly, no divergence with early instrumental data and during the last decades is found.

AD 751-1995 climate reconstruction

The reconstructions of summer temperature was made using regression weight for the Histalp June-August mean temperature and the full 1760-1995AD calibration period (Tab. 2). The reconstructions are expressed as temperature (1/10 °C) anomalies from the 1961-1990AD baseline period. The new reconstruction shows summer temperature variation in the Briançonnais region over the last 1,200 years with major confidence during the period's 900-1210AD and 1413-1995AD (Fig 8A). The maximum amplitude is about 5°C with the warmest summers occurring around 1990AD and the coldest summers occurring around 1810AD. Figure 6A shows that, on the decadal-to-centennial timescales, notably cold periods occur around 850, 1150850, , 1520850, , 1590850, , 1650850, , 1690850, , 1810850, , 1910 and 1970AD. and 1970AD. Notably warm periods occur around 1020, 1150, 1470 and 1990. The 60-year long warm period centred on 1150AD is especially noticeable while it appears as warm as the late twentieth century.

Discussion

Briançonnais climate history and possible bias in temperature reconstruction

A composite larch chronology is used to reconstruct past temperatures in the Briançonnais region until 751AD. This chronology yielded a high correlation to JJA summer temperatures (Tab. 2) as demonstrated by several authors (e.g. [START_REF] Serre | The dendroclimatological value of European larch (Larix Decidua Mill.) in the French Maritime Alps[END_REF][START_REF] Carrer | Age-dependent tree-ring growth responses to climate in Larix decidua and Pinus cembra[END_REF][START_REF] Rolland | Changes in radial tree growth for Picea abies, Larix decidua, Pinus cembra and Pinus uncinata near the alpine timberline since 1750[END_REF][START_REF] Frank | Temperature reconstructions and comparisons with instrumental data from a tree-ring network for the European Alps[END_REF].

The correlation significantly increases after smoothing the ring width record using a 10-year low pass filter. This shift is related to the periodic (8-9 years) population waves of the larch budmoth (LBM, Zeiraphera diniana Guénée) and their feeding on needles, which typically cause growth reductions in the year and generally subsequent years from the attack and affect the climatic signal particularly in the higher frequency, inter-annual domain. For this reason, the climate description was only conducted in the medium frequency domain.

The Briançonnais temperature history reports warm summers during much of the tenth, twelfth and late thirteenth centuries (Fig. 6E) commonly associated with the late Medieval Warm Period [START_REF] Hughes | Was there a ''Medieval Warm Period'', and if so, where and when?[END_REF][START_REF] Esper | Low-frequency signals in long tree-ring chronologies for reconstructing past temperature variability[END_REF]. This period, associated with the Medieval Warm Period, is characterized by significant inter-decadal temperature variations, with warm decades centered around 1020, 1160 and 1240. High summer temperatures are also inferred prior to AD 1000, but caution with this observation is warranted as this period is characterized by low sample replication (Fig. 8A) and low EPS values.

In our reconstruction, the Little Ice Age (LIA) is well marked and the cold period extends from 1500 to 1850AD but the transition MWP-LIA remains temporally imprecise due to the gap and the associated uncertainties in the reconstruction between 1220 and 1420AD. Cooler conditions during the LIA are accompanied by significant inter-decadal variations. Low temperatures are reconstructed in the 1520s, around 1600, around 1700 and 1815-1820. The inter-decadal variations seen since 1760AD are in line with JJA temperatures recorded in the Alps (Fig. 7A,B). We do not notice any misfit between colder reconstruction and warmer instrumental data before ~1840AD as related in several papers (Büntgen et al., 2005[START_REF] Büntgen | Summer temperature variations in the European Alps, AD 755-2004[END_REF][START_REF] Frank | Warmer early instrumental measurements versus colder reconstructed temperatures: shooting at a moving target[END_REF][START_REF] Corona | Millennium-long summer temperature variations in the European Alps as reconstructed from tree rings[END_REF][START_REF] Corona | Long-term summer (751-2008) temperature fluctuation in the French Alps based on tree-ring data[END_REF]. This improvement may be related to the adjustments made on the HISTALP dataset to take into account sheltering-related biases before 1850AD.

Since the 1820s, temperatures increased until the end of the record in 1995AD. Because of uncertainties in the proxy-instrumental temperature calibration, it is still difficult to unequivocally assert that the late 20 th century warming is significantly greater than the peak warmth of the Medieval Warm Period. But there is even less justification to assert the opposite-it is not possible to make a robust statement that the Medieval Warm Period was warmer than the last two decades of the 20 th century.

Comparison with other climate reconstructions

Among the previous reconstructions of late-Holocene climate changes available in the Alps, we selected, as reference climate records to validate the Briançonnais temperature reconstruction: (i) the Great Aletsch glacier history [START_REF] Holzhauser | Glacier and lake-level variations in west-central Europe over the last 3500 years[END_REF]; (ii) the Lake Anterne temperature reconstruction derived from chironomid assemblages [START_REF] Millet | Late-Holocene summer temperature reconstruction from chironomid assemblages of Lake Anterne[END_REF]; (iii) the speleothem-based temperature reconstruction in Spannagel Cave (Austria, [START_REF] Mangini | Reconstruction of temperature in the Central Alps during the past 2000 yr from a d18O stalagmite record[END_REF]; (iv) the reconstruction of River Rhône flooding activity from Lake Bourget [START_REF] Arnaud | 7200 years of Rhône river flooding activity recorded in Lake Le Bourget: a high resolution sediment record of NW Alps hydrology[END_REF]; (v) the summer temperatures reconstructed from larch TRW and Maximum Density (MXD) series from the Swiss and Austrian Alps (Büntgen et al., 2005[START_REF] Büntgen | Summer temperature variations in the European Alps, AD 755-2004[END_REF], Past variations in the length of the Great Aletsch glacier are well-documented over the last 3500 years on the basis of tree-ring studies of subfossil wood found in situ in the glacier forefield and its surroundings. This chronology is supported by radiocarbon and tree-ring dating, moraine studies and annual direct observations since 1892 [START_REF] Holzhauser | Glacier and lake-level variations in west-central Europe over the last 3500 years[END_REF]. In Lake Anterne (2060 m asl), [START_REF] Millet | Late-Holocene summer temperature reconstruction from chironomid assemblages of Lake Anterne[END_REF] presented a chironomid-based reconstruction of late-Holocene temperature. Chironomid assemblages were studied in 49 samples along an 8 m long sediment core covering the last 1800 years. July air temperatures were inferred using an inference model based on the distribution of chironomid assemblages in 100 Swiss lakes. The Mangini summer reconstruction is based on the precisely dated isotopic composition of a stalagmite from Spannagel Cave (2524 m asl) in the Central Alps of Austria. The stalagmite was dated using the U/Th method and δ18O profile was converted into yearly temperatures. Despite it concerns the whole year, this reconstruction is used for comparison because it is fully independent from our dataset. In Lake Bourget, [START_REF] Arnaud | 7200 years of Rhône river flooding activity recorded in Lake Le Bourget: a high resolution sediment record of NW Alps hydrology[END_REF] reconstructed the Holocene evolution of the River Rhône flooding activity by measuring the typical mineralogical fraction transported by the river and diluted by the autochthonous carbonate fraction within the lake. This series covers the last 7200 years and was dated with nine AMS 14C measurements. Finally, annually resolved summer temperatures were reconstructed for the European Alps at high temporal resolution by Buntgen et al. (2005Buntgen et al. ( , 2006)). The reconstructions repectively cover the periods 951-2003AD (755-2004AD) and are based on recent and historical larch ring width (maximum density) series.

Changes in tree-ring-inferred temperatures were within the prediction error of the model. However, the decadal to centennial temperature variability inferred from the Briançonnais tree-ring composite chronology matches the climate history inferred from other proxies in the northern and central Alps (Fig. 8) as follows:

• The high temperatures reconstructed in the Briançonnais between 900 and 1200AD match with a major glacier retreat in the northern Alps, a period of reduced Rhône flooding activity and high temperature reconstructed from chironomids, speleothems and larch series. In the detail, between 1020 and 1180AD, the tree-ring reconstructions are remarkably synchronous with a cold phase in the mid-11 th century and hot summers reconstructed around 1180AD. The sharp decline in the reconstructed temperatures around 1050AD is related to the Oort solar minimum [START_REF] Büntgen | Summer temperature variations in the European Alps, AD 755-2004[END_REF]. It is in phase with the beginning of a Great Aletsch glacier advance event reported by [START_REF] Holzhauser | Glacier and lake-level variations in west-central Europe over the last 3500 years[END_REF] but is observed neither in speleothems nor in lake sediments. Between 1150 and 1220AD, hot summers reconstructed from tree-rings coincide with a peak in July temperatures reconstructed from chironomid assemblages in lake Anterne [START_REF] Millet | Late-Holocene summer temperature reconstruction from chironomid assemblages of Lake Anterne[END_REF]) and predate a major phase of glacier advance (1250AD).

• The LIA lasting from 1500 to 1850AD in the Briançonnais coincides with an Alpine glacier readvance, a period of increased sediment supply in Lake Bourget and a decrease in summer temperature reconstructed from other speleothems, chironomid assemblages and alpine larch series (Fig. 8). The cold period that occurs around 1690AD is also visible in other proxies and particularly with a major period of glacial advance and flooding activity. This so-called Maunder minimum period coincides with an enhanced concentration in atmospheric 14C [START_REF] Stuiver | Sun, ocean, climate and atmospheric 14CO2: an evaluation of causal and spectral relationships[END_REF], several large volcanic eruptions [START_REF] Briffa | Influence of volcanic eruptions on Northern Hemisphere summer temperature over the past 600 years[END_REF], a reduced solar activity, and a low number of sunspots [START_REF] Spörer | Über die Periodizität der Sonnenflecken seit dem Jahre 1618, vornehmlich in Bezug auf die heliographische Breite derselben, und Hinweis auf eine erhebliche Störung dieser Periodizität während eines langen Zeitraumes[END_REF][START_REF] Maunder | The Prolonged Sunspot Minimum 1675-1715[END_REF][START_REF] Eddy | The Maunder Minimum[END_REF][START_REF] Beer | Reconstruction of Solar Irradiance since 1610: Implications for Climate Change[END_REF]. During the Dalton solar minimlum (1810), a series of tropical eruptions (1808-15), likely resulted in an aerosol-accumulated summer cooling effect [START_REF] Chenoweth | Two major volcanic cooling periods derived from global marine air temperature, AD 1807-1827[END_REF]Dai et al.1991). This volcanic forcing has been demonstrated to be largely responsible for the drop recorded in tree-ring series and glacier advance whereas changes in solar forcing and the increasing atmospheric CO2 concentrations were of minor importance [START_REF] Wagner | The influence of volcanic, solar and CO2 forcing on the temperatures in the Dalton Minimum (1790-1830) : A model study[END_REF].

• Since the mid-18 th century, the Briançonnais temperature reconstruction tracks the post LIA warming recorded in meteorological series (Fig. 7) and attributed to the increasing forcing induced by human activities [START_REF] Wanner | Mid-to Late Holocene climate change: an overview[END_REF], such as rapid land cover change, increase of greenhouse gases and aerosols, and stratospheric ozone depletion. This evolution is synchronous with the Great Aletsch Glacier retreat and is recorded in many other climatic records from the Alpine region (e.g., [START_REF] Holzhauser | Glacier and lake-level variations in west-central Europe over the last 3500 years[END_REF][START_REF] Mangini | Reconstruction of temperature in the Central Alps during the past 2000 yr from a d18O stalagmite record[END_REF].

Conclusion

In the Briançonnais, the tree-ring reconstruction of summer temperature indicates a high variability of summer temperature in the northern French Alps during the last 1250 years. At a hemispheric or global scale the existence of the LIA and MWP have been questioned [START_REF] Jones | Climate over past millennia[END_REF]. The present temperature reconstruction, however, together with the reconstruction of Alpine glacier fluctuations, major Rhône river flooding activity, summer temperature derived from Larix decidua, speleothems and chironomid assemblages confirm the existence of several significant climatic changes during the last 1250 years in the French Alps and more specifically the MWP and the LIA.

Moreover, the Briançonnais reconstruction indicates a relatively high magnitude of change in the inferred summer temperature. The temperatures reconstructed during the hottest periods of the MWP reach values close to the 1980-1995 mean JJA temperature. In comparison with conditions prevailing during MWP, the LIA coldest events were characterized by a cooling of about 1.5 (Maunder Minimum) and 1.3°C (Dalton Minimum). These results match well with the magnitude of temperature changes reconstructed in the Alps for the LIA cooling from other proxies, such as oxygen isotopic composition in speleothems [START_REF] Mangini | Reconstruction of temperature in the Central Alps during the past 2000 yr from a d18O stalagmite record[END_REF] and the Larix decidua ring width and density series (Büntgen et al., 2005[START_REF] Büntgen | Summer temperature variations in the European Alps, AD 755-2004[END_REF].

A detailed discussion of the relationship between human society and climate, the complexity of which has been stressed by previous studies (e.g., [START_REF] Berglund | Human impact and climate changes -synchronous events and a causal link?[END_REF], lies well beyond the scope of this article. However, the results obtained in the Briançonnais provide the opportunity to briefly point to the possibility of impacts of climate variations on human activities in such an elevated Alpine area. Several authors have already discussed the coincidence of climatic reversals and reduced human activity in sensitive regions (e.g. [START_REF] Berglund | Human impact and climate changes -synchronous events and a causal link?[END_REF]. In the Briançonnais, large fluctuations in summer temperature as shown in this study might have at least partly influenced the land-use activities in alpine highelevation sites. Thus, the lack of archeological wood around 1300AD coincides with the transition between MWP and LIA recorded in other proxies, says a major advance in Grindewald glacier length and a strong cooling recorded in chironomid assemblages from Lake Anterne, tree-ring series in Central Europe and oxygen isotopic composition in speleothems. Interestingly, this period also corresponds to the end of mining activities in the high altitude mine of Brandes (neighbouring Oisans massif) probably triggered by increasing environmental and climatic constraints (Bailly-Maître, 1996). Table 1. Characteristics of the tree-ring sites: site name, location, period covered by the chronology, exposition, number of individual tree-ring width series (Nb series), Mean Serie Length (MSL in years), Average Growth Rate (AGR in mm) and source (DENDRODB: http:// dendrodb.cerege.fr).

Captions

Figure 3. Correlation matrix of the tree ring series from the Buffères buildings and local or regional master chronologies (Dendron II software). . [(B), inset to a] Moving 31-year correlations with current-year June (light grey), July (dark grey), August (black), and JJA (dotted black) mean temperatures. Horizontal lines denote 95% and 99% significance levels.

Table 2. Calibration and verification results

Figure 7. The regression equation was tested in a cross calibration/verification exercise using the RCS-chronology and the HISTALP mean summer (June-August) temperature. The data were split in two 118-year periods: 1760-1877 and 1878-1878. The regression coefficients derived from the first period were used to reconstruct temperatures in the second period, and vice versa, thus producing two sets of independent data. The instrumental temperature is shown in blue and estimates in red for unfiltered and 20-year low pass dataset.

Figure 8. Comparison of the Briançonnais reconstruction with regional temperature reconstructions. All reconstructions were 20-yr low-pass filtered. (A) our raw (grey) and low-pass filtered (black) reconstruction is compared with : (B) the advance/retreat history of the Great Aletsch glacier (Switzerland, Holzhauser et al., 2005); (C) the chironomid-based temperature reconstruction at Lake Anterne (France, Millet et al., 2009) and the speleothem-based temperature reconstruction in Spannagel Cave (Austria, [START_REF] Mangini | Reconstruction of temperature in the Central Alps during the past 2000 yr from a d18O stalagmite record[END_REF]; (D) the reconstruction of River Rhône flooding activity from Lake Bourget (France, [START_REF] Arnaud | 7200 years of Rhône river flooding activity recorded in Lake Le Bourget: a high resolution sediment record of NW Alps hydrology[END_REF]; (E) the Alpine temperatures derived from Larix decidua ring width and wood density series [START_REF] Büntgen | Summer temperature variations in the European Alps, AD 755-2004[END_REF]. All tree-ring temperature reconstructions were transformed to z scores over the 951-1995AD common period.

Figure 1 .

 1 Figure 1. Localisation map of the tree-ring sites used for the reconstruction. All tree-ring sites are above 1900 m asl excepted Puy-Saint-Vincent.

Figure 2 .

 2 Figure 2. Multi-centennial living larches in the upper part of the Oriol forest (A) and bearing wall of a mountain summer barn composed of a typical assemblage of larch boles, wood pieces on wood pieces (Névache, Hautes-Alpes).

Figure 4 .

 4 Figure 4. A. Temporal distribution of the 376 historical (blue) and 185 living trees (red) samples ordered by calendar age of their innermost ring. B. Mean cambial age of the historical and living trees samples for each calendar year. C, D. Regional curves (RCs) of the age-aligned living trees (C) and historical (D) TRW series. The number of replicate series for each year (the sample depth) is shaded in grey.

Figure 5 .

 5 Figure 5. Fig. 5 RCS standardized chronologies for historic (blue, A), living trees (red, B) and composite TRW (black, E), with their corresponding EPS (C) and r-bar (D) values plotted for 50-year windows with 25 years overlap. The number of replicate series for each year (the sample depth) is shaded in grey.

Figure 6 .

 6 Figure6. Climate response of the composite RCS chronology. (A) Correlations with the HISTALP monthly mean temperatures (black) and precipitation (gray) of the previous and current year(1760- 1995). [(B), inset to a] Moving 31-year correlations with current-year June (light grey), July (dark grey), August (black), and JJA (dotted black) mean temperatures. Horizontal lines denote 95% and 99% significance levels.

Apr (n-1)


Louis Lagardère
email: louis.lagardere@sorbonne-universite.fr

Félix Aviat

Jean-Philip Piquemal
email: jean-philip.piquemal@sorbonne-universite.fr

Pushing the limits of Multiple-Timestep Strategies for Polarizable Point Dipole Molecular Dynamics

published or not. The documents may come   L'archive ouverte pluridisciplinaire

approaches by coupling these splittings to Langevin dynamics and to Leimkuhler's BAOAB integrator in order to reach larger time-steps (6 fs) for long-range forces. We further increase sampling efficiency by: i) accelerating the polarization evaluation using a fast/non-iterative Truncated Conjugate Gradient (TCG-1) as short-range solver; ii) pushing the outer time-step to 10 fs using hydrogen mass repartitioning. The new BAOAB-RESPA1 integrators demonstrate up to a 7-fold acceleration over standard 1 fs (Tinker-HP) integration and reduce the performance gap between polarizable and classical force fields while preserving static and dynamical properties.

The most straightforward way to speedup molecular dynamics (MD) [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Mccammon | Dynamics of Proteins and Nucleic Acids[END_REF] is to use larger time-steps. In this context, multi-timestep schemes emerged 3 , but the largest usable timestep is limited by resonance effects [START_REF] Biesiadecki | Dangers of Multiple Time Step Methods[END_REF][START_REF] Ma | Verlet-I/R-RESPA/Impulse is Limited by Nonlinear Instabilities[END_REF] . As pointed out by various authors, it is possible to overcome these effects by using modified dynamics that still sample the correct measure, but these solutions alter the dynamical properties (Generalized Langevin Equation (GLE) [START_REF] Morrone | Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances[END_REF] , stochastic isokinetic extended phase-space algorithm [START_REF] Leimkuhler | Stochastic, resonance-free multiple time-step algorithm for molecular dynamics with very large time steps[END_REF][START_REF] Margul | A Stochastic, Resonance-Free Multiple Time-Step Algorithm for Polarizable Models That Permits Very Large Time Steps[END_REF][START_REF] Albaugh | Combining Iteration-Free Polarization with Large Time Step Stochastic-Isokinetic Integration[END_REF] ). However, in practice, one would like to accelerate MD while also preserving the dynamic [START_REF] Morrone | Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances[END_REF][START_REF] Feenstra | Improving efficiency of large time-scale molecular dynamics simulations of hydrogen-rich systems[END_REF] . This letter addresses this problem in the particular context of polarizable force fields (PFF) 11,[START_REF] Gresh | Anisotropic, Polarizable Molecular Mechanics Studies of Inter-and Intramolecular Interactions and Lig-andâˆ'Macromolecule Complexes. A Bottom-Up Strategy[END_REF] . This class of methods is more computationally expensive than classical force fields (FF) because of the need to evaluate a many-body polarizable energy [START_REF] Lipparini | Scalable evaluation of polarization energy and associated forces in polarizable molecular dynamics: I. Toward massively parallel direct space computations[END_REF][START_REF] Lagardère | Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: II. Toward Massively Parallel Computations Using Smooth Particle Mesh Ewald[END_REF] . Multi-timesteping is therefore essential. The general consensus to ensure conserved properties is to limit the use of reversible Reference System Propagator Algorithm (RESPA) integrators [START_REF] Tuckerman | Reversible multiple time scale molecular dynamics[END_REF] to a a bonded/non-bonded forces split and to use a 2 fs time-step for the non-bonded forces. Further splitting of the non-bonded forces is not straightforward, because of the many-body nature of polarization, but has been shown to be applicable [START_REF] Margul | A Stochastic, Resonance-Free Multiple Time-Step Algorithm for Polarizable Models That Permits Very Large Time Steps[END_REF][START_REF] Masella | The multiple time step r-RESPA procedure and polarizable potentials based on induced dipole moments[END_REF] . Indeed one can define a short-range polarization energy and evaluate, at an outer-timestep, the slowly varying difference between the actual polarization energy (and forces) and the short-range ones. More precisely, one has to evaluate both the short-range and total polarization terms at these outer time-steps. The reduced computational cost of the short-range polarization contribution and the less frequent evaluation of the total one effectively reduces the computational effort. Since the upper limits of these strategies have not yet been evaluated by the community, we will, in this letter, assess this frontier to improve simulation performances while respecting two important constraints: i) the mandatory need to preserve static and dynamical properties; ii) the possibility of a black-box implementation allowing for strong computational speedups without dependence to the studied system. In everything that follows tests have been made using the AMOEBA polarizable force field [START_REF] Ren | Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation[END_REF] and the Tinker-HP software [START_REF] Lagardère | polarizable force fields[END_REF] . Technical details as well as various algorithmic setups are provided in Supplementary Information (see section S1). A summary of our incremental strategy is depicted on Figure 1. Interested developers can also look at the code that will be available on the Tinker-HP website [START_REF]Tinker-HP's Website[END_REF] and later on Github [START_REF]Tinker-HP's Github[END_REF] .

Velocity-Verlet (reference)

Split forces (bonded/non -bonded)

V-RESPA

Faster computation

Change integrator

BAOAB-RESPA1 Larger timestep

Conserved properties

Split non-bonded contributions

V-RESPA1

Larger timestep Degraded properties

Less expensive polarization solver

BAOAB-RESPA1 TCG

Faster computation

Hydrogen Mass Repartitioning

BAOAB-RESPA1 TCG + HMR

Larger timesteps

Figure 1: 'V-' (respectively 'BAOAB-') indicates that the numerical integration scheme is Velocity-Verlet 20 (respectively BAOAB [START_REF] Matthews | Rational Construction of Stochastic Numerical Methods for Molecular Sampling[END_REF][START_REF] Leimkuhler | Robust and efficient configurational molecular sampling via Langevin dynamics[END_REF] ). 'RESPA' and 'RESPA1' respectively mean the RESPA single-split (bonded vs. non-bonded) strategy 3 and the RESPA1 double-split (bonded, short-range non-bonded, long-range non-bonded) one [START_REF] Zhou | Efficient multiple time step method for use with Ewald and particle mesh Ewald for large biomolecular systems[END_REF] . 'TCG' is the acronym for Truncated Conjugate Gradient, a fixed-cost non-iterative polarization solver 24 . 'HMR' stands for Hydrogen Mass Repartitioning 10 , implemented to avoid high-frequency motions.

A popular integrator: V-RESPA. Let us first evaluate the limits of the bonded/nonbonded RESPA integrators for which all the bonded terms are evaluated within a Velocity-the non-bonded terms (van der Waals, electrostatics and polarization) are evaluated first at 2 and then at 3 fs. To assess the accuracy of the associated integrators we ran simulations of 2 ns in the NVT ensemble at 300 K with two test systems: a cubic box of 500 water molecules, with a 24.66 Angstroms edge, and a 9737 atoms box with edges of 54.99×41.91×41.91 Angstroms containing a solvated protein (the ubiquitin). In both cases, periodic boundary conditions for electrostatics and polarization were evaluated with Smooth Particle Mesh Ewald (SPME) [START_REF] Lagardère | polarizable force fields[END_REF][START_REF] Essmann | A smooth particle mesh Ewald method[END_REF][START_REF] Stamm | A coherent derivation of the Ewald summation for arbitrary orders of multipoles: The self-terms[END_REF] with standard parameters (see SI) as we chose a Preconditioned Conjugate Gradient (PCG) polarization solver using a diagonal preconditionner and a 10 -5 convergence threshold [START_REF] Lipparini | Scalable evaluation of polarization energy and associated forces in polarizable molecular dynamics: I. Toward massively parallel direct space computations[END_REF][START_REF] Lagardère | Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: II. Toward Massively Parallel Computations Using Smooth Particle Mesh Ewald[END_REF] . For each of these systems and for each integrator, we computed various static observables: average potential energy, average polarization energy and for the bulk water system , the radial distribution functions. In this last case, we also computed the self-diffusion coefficient, a dynamical property evaluated with the Einstein formula by averaging over a series of time origins [START_REF] Maginn | Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics [Article v1.0[END_REF] . The self-diffusion coefficient is known to have a size dependency vanishing at the infinite size limit [START_REF] Maginn | Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics [Article v1.0[END_REF] , but here these values are only used as means of comparison between integrators, hence these corrections were not applied.

These tests are performed in the canonical ensemble for which the choice of the thermostat impacts the dynamics of the system. We ran these tests using the global velocity rescaling thermostat developed by Bussi et al. with a relaxation constant of 0.2 ps, for which the dynamical properties are close to the one obtained with pure hamiltonian dynamics [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] . These values have been compared to the ones obtained with a Velocity-Verlet integrator used at a 0.5 fs time-step, which can be considered as a reference. In the rest of the text, we will be denoting the different time-step lengths as a/b, a being the bonded terms one and b the non-bonded terms one (both in fs).

For both systems, the V-RESPA integrator where non-bonded forces are evaluated at 2 fs gives similar results as the reference (within statistical uncertainty) with a difference of less than two percents for average energies (see Tables 1-3 of the SI). With an outer time-step the error on the polarization energy grows significantly (more than 2.5%). This advocates for a careful use of this setup.

Concerning O-O radial distribution function for water, no significant differences with the reference Velocity-Verlet (0.5fs) ones can be observed among these different methods (see Figure 2). The self-diffusion coefficients of water is also nicely preserved for the V-RESPA integrator with a 2 fs outer time-step, though it is slightly off (around eight percent) with a 3 fs outer time-step (see Table 2, SI).

Further range-separation in polarizable force fields: V-RESPA1. As a second part, we will now evaluate the limits of other RESPA integrators, for which the non-bonded terms are further split in two parts, the short-and long-range. We are now considering three terms:

the bonded, the non-bonded short-range, and the non-bonded long-range terms. Regarding the split of the electrostatics and the polarization energies, we chose to use the RESPA1 logic [START_REF] Zhou | Efficient multiple time step method for use with Ewald and particle mesh Ewald for large biomolecular systems[END_REF] , where the short-range part of the electrostatic (and polarization) energy is defined as the short-range part of the real space SPME energy. In this case, it has been shown that the stability of the integrator is less dependant on the smoothing parameters used to switch between short-and long-range [START_REF] Morrone | Molecular Dynamics with Multiple Time Scales: How to Avoid Pitfalls[END_REF] . Details on the definition of these short and long range forces as well as these smoothing parameters can be found in SI. We test various setups within this context: the bonded forces are always evaluated every 0.25 fs, but the short-range non-bonded ones are either evaluated every 2 or 2.5 fs, and the time-step of the long-range forces (that has to be a multiple of the previous one) is either 4, 5 or 6 fs. In the rest of the text, these integrators will be denoted as V-RESPA1. We will be denoting the different time-steps lengths of the integrators as a/b/c, a being the bonded terms time-step length, b the short-range non-bonded and c the long-range ones (all in fs). For the bulk water system (see Table 4 in SI), we observe that both the average potential and polarization energies are preserved within two percents of the reference value for the 0.25/2/4 and the 0.25/2.5/5

setups, but that the average polarization energy is more than two percents off for the 0.25/2/6 setup. Concerning the radial distribution functions of water, it is clear that only the 0.25/2/4

integrator gives satisfactory result as other choices diverge from the reference, as can be seen in Figure 2. Furthermore, if the self-diffusion coefficient is stable for the 0.25/2/4 integrator (see Table 5, SI), it exhibits a dramatic decrease for the other ones (falling at 1.34 instead well preserved with these setups, but also that the computational gains expected due to the use of a larger time-step are counterbalanced by a lower sampling rate [START_REF] Morrone | Efficient multiple time scale molecular dynamics: Using colored noise thermostats to stabilize resonances[END_REF][START_REF] Feenstra | Improving efficiency of large time-scale molecular dynamics simulations of hydrogen-rich systems[END_REF] .

Indeed, as pointed out by Berendsen 6,10 , such a decrease in the self-diffusion coefficient is expected to reduce the sampling efficiency by a similar amount because it is associated to an increase of water viscosity and thus to a slowing down of large scale motions. For the solvated ubiquitin, it is also clear that only the 0.25/2/4 setup corresponds to a satisfactory accuracy, as the other ones give average potential and polarization energies off by more than three percents (see Table 6, SI). ensemble is Langevin dynamics where the coupling to a heat bath is made through additional local frictions and dissipative stochastic terms. The dynamics is then known to be altered compared to pure Hamiltonian one, but this impact is expected to be small for relatively small friction constant (less than 1 ps -1 ). In this context, Leimkuhler and collaborators proposed an integrator for Langevin dynamics based on operator splitting and a particular ordering of the terms of the equations of motion named BAOAB [START_REF] Matthews | Rational Construction of Stochastic Numerical Methods for Molecular Sampling[END_REF][START_REF] Leimkuhler | Robust and efficient configurational molecular sampling via Langevin dynamics[END_REF] . They showed in various contexts [START_REF] Matthews | Rational Construction of Stochastic Numerical Methods for Molecular Sampling[END_REF] that this integrator has an improved accuracy for configurational sampling compared to other ones. We thus also tested the previously presented splittings using this new integrator, with a 1 ps -1 friction constant (they will be denoted BAOAB-RESPA and BAOAB-RESPA1 in the rest of the text ; see SI for an additional description of these integrators) and noted a significant improvement in terms of accuracy (reported on Tables 7-9 of SI). Indeed, for the bulk water system (Table 7 of SI), the errors for the BAOAB-RESPA integrator with a bonded/non-bonded split and a 3 fs outer time-step is limited to less than one percent for the average total potential energy and two percents for the average polarization energy, staying in both case within statistical error, compared to 1.1 and 2.8%

for the Velocity-Verlet based integrator. Figure 2 also shows an improved agreement with the reference for the water radial distribution compared to RESPA. The same behaviour can be observed on the solvated ubiquitin for which both these values stay respectively below one and two percents -within statistical error (see Table 9, SI). Because the dynamics is modified when running NVT trajectories with Langevin, even if the differences are expected to be small for a small friction, comparing values of the self-diffusion coefficients obtained with these integrators only makes sense by taking as a reference a numerical scheme integrating Langevin dynamics with conservative parameters. This is why we chose, as a reference for these values, the ones obtained with a 1 fs BAOAB integrator and a 1 ps -1 friction, which as expected gives a self-diffusion close to the reference Velocity-Verlet one (1.82 versus 2.08). For the BAOAB-RESPA integrators, we see that errors on the self-diffusion coefficients (see Table 8, SI) are limited to six percents with a 3 fs outer time-step compared to eight percents with a similar time-step and a Velocity-Verlet inner loop. The better performances of BAOABbased integrators with respect to the Velocity-Verlet ones becomes obvious within a RESPA1 split. Indeed, we computed the same observables for the BAOAB-RESPA1 integrators (see Tables 10-12, SI and Figures 34), equivalent to the V-RESPA1 integrators, and we see that the average potential and polarization energies are strikingly more stable and always within two percent error with respect to their reference value, that is to say within statistical uncertainty. Similar comments can be made on the radial distribution functions: unlike the Velocity-Verlet integrators, they almost perfectly overlap with their references even with a 5 or 6 fs outer time-step. When using the Velocity-Verlet based RESPA1 integrators with an outer-time-step larger than the most conservative (4 fs) one, the diffusion coefficient showed a dramatic decrease (see Figure 3, and Table 5 in SI). Yet in the BAOAB-RESPA1 case, this dynamical observable remains far more stable for all integrator setups: even for the choice of evaluating long-range non-bonded forces every 6 fs, the error is less than eight percents, whereas it exceeds thirty-five percents in the equivalent V-RESPA1 setup. This highlights again that BAOAB-RESPA1 integrators are not only more accurate, but also ensure a conserved sampling rate, which is not the case for the V-RESPA1 ones.

affected by this choice [START_REF] Lipparini | Scalable evaluation of polarization energy and associated forces in polarizable molecular dynamics: I. Toward massively parallel direct space computations[END_REF][START_REF] Lagardère | Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: II. Toward Massively Parallel Computations Using Smooth Particle Mesh Ewald[END_REF][START_REF] Aviat | Truncated Conjugate Gradient: An Optimal Strategy for the Analytical Evaluation of the Many-Body Polarization Energy and Forces in Molecular Simulations[END_REF] .

A this point, we reached the performance limits if one wants to preserve a tight accuracy on the dynamics. One of the most natural ways to further increase the size of the usable time-step when simulating a large biological system is to redistribute the mass of the heavy atoms on the hydrogens they are carrying (a method named Hydrogen Mass Repartioning, or HMR 10 ), thus limiting the high frequency stretching motions of these atoms while keeping the same configurational potential energy surface. In the following, we show that this redistribution allows to use even larger time-steps while maintaining satisfactory accuracy with a BAOAB-RESPA1 integrator and the same TCG-1 short-range polarization solver as before (Table 16-18, SI). As can be seen in Table 1, the approach appears to be a very good compromise: large speedups can be obtained by pushing the bonded force time-step to 1 fs, the short-range non-bonded forces time-step up to 10 3 fs and the outer one up to 10 fs. Details on how the mass repartitioning is done can be found in SI (section S1). A resulting acceleration of 4.72 (6.8 without ASPC) is obtained, keeping the errors on the average energies below two percents, maintaining an accurate evaluation of radial distribution functions and a good enough evaluation of the self-diffusion coefficient so that sampling efficiency is preserved.

Since PCG, as a Krylov method, is systematically improvable [START_REF] Lipparini | Scalable evaluation of polarization energy and associated forces in polarizable molecular dynamics: I. Toward massively parallel direct space computations[END_REF] , additional small speedups can be obtained by focusing on the long range PCG solver performances. For example, besides using a diagonal preconditioner, one could use more advanced techniques such as those proposed by Skeel [START_REF] Wang | Fast evaluation of polarizable forces[END_REF] or by Beran 34 . Improved performances of three to four percents are observed, reaching a global acceleration of more than seven, with the same accuracy as a 1 fs Velocity Verlet scheme without ASPC. Finally, beside the net acceleration, another advantage of the TCG use lies in the absence of use of a dipole history (as in predictor-correctors such ASPC), leading to a method free of time-reversibilty and volume preservation issues [START_REF] Aviat | Truncated Conjugate Gradient: An Optimal Strategy for the Analytical Evaluation of the Many-Body Polarization Energy and Forces in Molecular Simulations[END_REF] .

Finally, the small decrease (8 %) of the self-diffusion constant (see Figure 4) observed in the most aggressive setup has to be compared with actual available large step methods 9 that, despite their qualities, are not able to maintain accuracy on dynamical properties provid-ing diffusion constants reduced by a factor of 5 9 . Our approach does not suffer from these problems: it remains operational, maintaining sampling efficiency.

To illustrate the robustness of these approaches, we performed several tests taking advantage of our massively parallel AMOEBA production implementation in Tinker-HP [START_REF] Lagardère | polarizable force fields[END_REF] . First, we checked the stability of the dynamics using the fastest available setup. We provide a 15 ns simulation of ubiquitin (see SI, section S3): the potential and polarization energies normally fluctuate around their mean values, demonstrating the stability of the approach. Furthermore, we computed the average molecular dipole moments for the bulk water systems and confirmed their full stability (see SI, section S4). Second, we ran simulations on large systems of biological interest, namely the solvated dihydrofolate reductase protein (DHFR, 23358 atoms) and the solvated Satellite Tobacco Virus (STMV, 1066628 atoms). The discussed speedup of 7 (vs. a 1 fs/Velocity Verlet/PCG-10 -5 ) is conserved as we obtained a production of 22.2 ns/day on 680 cores for DHFR and 1.2 ns/day on 10800 cores for STMV. Such results are of major interest, as a 7-fold acceleration will enable to save millions of hours of computing time while enabling long and accurate polarizable molecular dynamics studies on large systems. Finally, we computed a more involved property which is of key importance in biological simulations: hydration free energies. We applied the Bennett Acceptance Ratio method [START_REF] Bennett | Efficient estimation of free energy differences from Monte Carlo data[END_REF] , a commonly used approach to compute free energies differences [START_REF] Rackers | TINKER 8: A Modular Software Package for Molecular Design and Simulation[END_REF] to evaluate the solvation free energy of a sodium cation in water. Results are shown in Table 2), and practical details on the choice of alchemical free energy difference windows can be found in SI (section S1). Even for the fastest setup, the values obtained are within 0.1 kcal/mol over 89.7 kcal/mol for the reference 37 , demonstrating the validity of these acceleration schemes and their capability to preserve accuracy.

To conclude, after examining the limits of a standard Velocity-Verlet integrator for Polarizable Force Fields used in combination with a RESPA1 split, we introduced new BAOAB-RESPA1 Langevin integrators coupled to fast short-range non-iterative TCG-1 polarization solver and Hydrogen Mass Repartitioning, achieving all together large computational speedups. Two optimal BAOAB-RESPA1 setups were presented and compared to a 1 fs Velocity-Verlet reference: i) one (namely 0.25/2(TCG)/6) for which all properties are preserved while providing a global speedup of more than four; ii) a second (1/ 10 3 (TCG)/10+HMR) for which dynamical properties are slightly affected but where sampling remains efficient, offering a strong acceleration up to a seven-fold. As accuracy is maintained and sampling efficiency is preserved while being system independent, the proposed methodology can be used as a black-box in our Tinker-HP framework, benefiting from its massive parallelism implementation and offering therefore further computational gains [START_REF] Lagardère | polarizable force fields[END_REF] . Such findings are game-changing, as they extend the applicability of polarizable Molecular Dynamics to longer timescale simulations and larger systems. In practice, the resulting performance gain helps reducing the computational gap between point dipole polarizable force fields such as AMOEBA and more tractable models such as Drude [START_REF] Jiang | High-Performance Scalable Molecular Dynamics Simulations of a Polarizable Force Field Based on Classical Drude Oscillators in NAMD[END_REF] or even non-polarizable force fields such as CHARMM [START_REF] Huang | CHARMM36m: an improved force field for folded and intrinsically disordered proteins[END_REF] or AMBER [START_REF] Maier | ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB[END_REF] .

Table 1: Speedup of BAOAB-RESPA and BAOAB-RESPA1 integrators calculated with respect to the Velocity-Verlet integrator at 1fs. The types of RESPA integrators are defined by: R1=RESPA1 and R=RESPA. Speedups obtained with V-RESPA and V-RESPA1 integrators are identical.*=replacement of the PCG diagonal preconditioner by an improved technique [START_REF] Wang | Fast evaluation of polarizable forces[END_REF][START_REF] Nocito | Fast divide-and-conquer algorithm for evaluating polarization in classical force fields[END_REF]

Figure 2 :

 2 Figure 2: Oxygen-oxygen radial distribution function for various integrators (see Text for notations). Radial distributions appear correct with most of the setups. However, degraded results are obtained with the V-RESPA1 integrators using large outer-timesteps beyond 4 fs.

Figure 3 :

 3 Figure 3: Average potential and polarization energies (in kcal/mol) for a 500-molecule water box and solvated ubiquitin computed using various integrators.

Figure 4 :

 4 Figure 4: Self-diffusion (D) coefficient for various integrators (in 10 -5 cm 2 s -1 ).

Table 2 :

 2 , see text. Hydratation free energies for the Na + cation

	Splits	0 .2 5 / 2	0 .2 5 / 3	0 .2 5 / 2 / 4	0 .2 5 / 2 .5 / 5	0 .2 5 / 2 / 6	0 .2 5 / 2 (T C G )/ 4	0 .2 5 / 2 .5 (T C G )/ 5	0 .2 5 / 2 (T C G )/ 6	1 / 10 3 / 1 0	1 / 10 3 (T C G )/ 1 0
	ASPC	1.75 2.37 1.72 2.43 2.53 2.32	2.7	2.9	4	4.72 (4.91*)
	No ASPC	2.53 3.42 2.5	3.5	3.7	3.4	3.9	4.2	5.8	6.8 (7.0*)
	RESPA-type	R	R	R1	R1	R1	R1	R1	R1 R1(HMR) R1(HMR)
	∆G hydrat Na + (kCal/mol)	Velocity-Verlet 0.5fs 1/ 10 3 /10-HMR 1/ 10 3 (TCG1)/10-HMR 89.7 (+/-0.13) 89.7 (+/-0.13) 89.6 (+/-0.13)

Acknowledgement

This work was supported in part by French state funds managed by CalSimLab and the ANR within the Investissements d' Avenir program under reference ANR-11-IDEX-0004-02. We thank GENCI for computer time within the Grand Challenge on the Iréne Joliot-Curie (TGCC, Bruyères le Châtel, France).

Concerning the effective speedups in our implementation, Table 1 displays the gains obtained for the BAOAB-RESPA and the BAOAB-RESPA1 integrators (compared to a regular 1 fs Velocity-Verlet). They are the same as the one obtained for the V-RESPA and the V-RESPA1 integrators. We show two entries in table 1: one where a guess based on Kolafa's Always Stable Predictor Corrector 30 (ASPC) is used for the induced dipoles (standard Tinker-HP setting), but only at short-range for the RESPA1 integrators, and one where the "direct field" guess is used [START_REF] Lipparini | Scalable evaluation of polarization energy and associated forces in polarizable molecular dynamics: I. Toward massively parallel direct space computations[END_REF][START_REF] Lagardère | Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: II. Toward Massively Parallel Computations Using Smooth Particle Mesh Ewald[END_REF] , showing that up to a 2.53 speedup (3.7 without ASPC) is achieved. Note that for the RESPA1 schemes, additional gains are made in the long-range polarization solvers by using, at the same time-step, the short-range dipoles obtained as guess for the long range ones, effectively reducing the number of iterations required to converge. For the BAOAB based integrators, the benefits of using the RESPA1 splitting are clearly demonstrated as the 0.25/2.5/5 and 0.25/2/6 frameworks are both faster than the 0.25/3 RESPA integrator for a similar accuracy.

Speeding up BAOAB-RESPA1: TCG-1 solver for the short-range polarization and Hydrogen Mass Repartioning. When using a RESPA1 multi-timestep integrator and a polarizable force field, the sole purpose of the short-range polarization energy is to eliminate the high-frequency part of the total polarization energy [START_REF] Margul | A Stochastic, Resonance-Free Multiple Time-Step Algorithm for Polarizable Models That Permits Very Large Time Steps[END_REF][START_REF] Liberatore | A Versatile Multiple Time Step Scheme for Efficient ab Initio Molecular Dynamics Simulations[END_REF] . This is why an approximate but less computationally expensive and non-parametric expression of the polarization energy can be used to fill this role and provide an additional speedup. In that context we decided to use the recently introduced Truncated Conjugate Gradient (TCG) [START_REF] Aviat | Truncated Conjugate Gradient: An Optimal Strategy for the Analytical Evaluation of the Many-Body Polarization Energy and Forces in Molecular Simulations[END_REF][START_REF] Aviat | The truncated conjugate gradient (TCG), a non-iterative/fixed-cost strategy for computing polarization in molecular dynamics: Fast evaluation of analytical forces[END_REF] as short-range solver. TCG can be chosen to be minimal in cost (TCG-1 with a diagonal preconditioner, but without any guess and without peek step) to be the fastest possible. This coupling provides an additional computational gain at a conserved accuracy (see Tables 131415, SI), corresponding to a final speedup of more than four times compared to a regular molecular dynamics of 1 fs with a Velocity-Verlet integrator (see Table 1). Such polarization setup offers full energy conservation and the static and dynamical properties are marginally

Supporting Information Available


François Javier Anez

Xavier Demoulin
email: demoulin@coria.fr

Nicolas Hecht

Julien Reveillon

François Xavier Demoulin

A general purpose LES Model for Atomization

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

Introduction

Fluid Catalytic Cracking (FCC) is a process where the crude oil is transformed into Gasoline, olefins and distillates. In this process, the liquid oil feed is atomized through a set of injectors connected circumferentially to a riser. Moreover, the oil droplets are vaporized and cracked within the riser upon contacting hot gases and catalyst, respectively [START_REF] Sadeghbeigi | Fluid catalytic cracking handbook[END_REF]. The main goal in an injector design is to produce small size droplets in order to ensure quick vaporization and intimate contact with the catalyst [START_REF] Grace | Circulating fluid beds[END_REF].

Most studies have assumed instant vaporization and uniform catalyst/oil ratio at the bottom of riser crosssectional area, based on the assumption that the cracking time scale is much longer than the time scale of interface transport in the feed injection zone [START_REF] Gan | Numerical simulation of hydrodynamics and cracking reactions in the feed mixing zone of multiregime gas-solid riser reactor[END_REF][START_REF] Rajesh | Effect of injection zone cracking on fluid catalytic cracking[END_REF]. Some studies considered the vaporization rate of the liquid oil using two approaches, namely homogeneous and heterogeneous mode, in which both methods set the limits for the actual vaporization time [START_REF] Tuyen | Comparison of vaporization models for feed droplet in fluid catalytic cracking risers[END_REF]. When it comes to FCC riser modeling many researchers have focused on the upper and middle zone of the riser and just a few researchers on nozzles and feed atomization zone [START_REF] Qilong | Experiment and Numerical studies on the atomization of a swirl feed nozzle[END_REF], mainly due to the high constraints using multidimensional models with chemical reactions, vaporization, and atomization within an industrial scale. Theologos [START_REF] Theologos | Feedstock atomization effects on FCC riser reactors selectivity[END_REF] studied atomization effects on reactor performance and found that smaller droplets produce higher vaporization rates. Therefore, the need for a deeper understanding of the physics and atomization mechanisms involved within a FCC injector coupled with the feed injection zone urged to be tackled. The aim of this project is to numerically model the oil-steam mixing inside the FCC injector.

Material and methods

The objective of this work is to extend LES-ELSA method [START_REF] Chesnel | Large Eddy Simulation for atomization: application to automotive engines injection[END_REF] into a general purpose solver within OpenFOAM® framework. A postulated transport equation of the interface density is implemented to describe the subgrid spray formation from interface wrinkling, ligament breakup, sheets to the breakup droplets. The interface capturing method is based on the transport of a liquid volume fraction (i.e. ). Accordingly, an additional compressive term is added to the equation to treat sharpened gradients. In LES, a subgrid-scale (SGS) model accounts for the dynamics of the unresolved scales of motion which induces a subgrid scale (i.e.

) that is not compatible with the numerical procedure used to capture the interface. The advantages of the proposed method is to disable the compressive term when becomes important (i.e. when the interface fluctuations become significant at subgrid-scale). Two criterions are used to switch from interface capturing method to subgrid-scale or vice versa, based on Interface Resolution Quality (i.e.

, which is derived from the transport equation of the interface density, ELSA algorithm).

In order to upgrade LES-ELSA method a coupling with Lagrangian particle was performed. Lagrangian particles are created when liquid structures become droplets. The diameter and the number of particle are determined via liquid volume fraction and via the interface density equation. The Eulerian model has been linked to the Lagrangian phase using the liquid turbulent diffusion flux closure.

Numerical Model

Atomization by definition is the process of converting a liquid form to a free gaseous atom. In other words, it is the transformation of a bulk liquid into a spray of liquid droplets in a surrounding gas or vacuum [START_REF]Corrosionpedia[END_REF]. Three different zones can be easily defined as shown in figure 1. The first region near the tip of the injector is called primary atomization. In this area the liquid volume fraction is close to one and the liquid surface topology is very complex. Here the liquid sheets experience longitudinal instabilities based on liquid and gas interactions (Kelvin-Helmholtz instabilities) which disturb the plane sheet in a sinusoidal stream-wise oscillation mode. the secondary atomization begins downstream of the flow, where instabilities turn into three dimensions, and the sheet breaks into smaller liquid packs, ligaments and bag-like structures. This zone ends with the formation of a spray of droplets (Dispersed zone) [START_REF] Davide | Numerical simulation of primary and secondary atomization[END_REF]. The ELSA model, which is devoted to describe atomization flows with high Weber and Reynolds values, has been adapted to account for both primary and secondary atomization zone by using an Euler approach. In this model the two phase flow is considered as single phase flow composed of two species, liquid and gas, with highly variable density [START_REF] Vallet | Euler model for liquid jet atomization[END_REF][START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modeling study[END_REF]. The dispersion of the liquid phase is calculated by a liquid volume fraction balance equation ( ) using a pondered inclusion of a 'turbulent diffusion liquid flux' [START_REF] Hecht | Large Eddy Simulation of liquid flow: Spray application[END_REF] and interface capturing method (ICM) used in OpenFOAM®, as shown in equation 1. Such ponderation is ruled by an arbitrary parameter ( ) which depends on the interface resolution scale.

(1)

The interface capturing method consists in defining an artificial and supplementary velocity field, ( ), in the vicinity of the interface, in such a way that the local flow steepens the gradient of the volume fraction and the interface resolution is improved [START_REF] Raees | Evaluation of the interface-capturing algorithm of OpenFOAM for the simulation of incompressible immiscible two-phase flow[END_REF]

Where '"#$' is the Schmidt number. It is important to note that this approach is valid only in the absence of a mean slip velocity between the phases [START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modeling study[END_REF]. The arbitrary parameter ( ) was set zero when the interface is poorresolved (Secondary atomization) and set to one otherwise (Primary atomization) as shown in figure 2. Two criterions have been used to properly set, , based on the Interface Resolution Quality (i.e. , which is derived from the transport equation of the interface density [START_REF] Hecht | Large Eddy Simulation of liquid flow: Spray application[END_REF]). The breakup process and the polydisperse spray distribution are accounted through the quantity, Σ & (' ( ), mean interfacial area per unit of volume, as shown in equation 3.

) & * + ) & = , ! ) & -+ ) & . /1 - ) & 2 * 4 (3)
The first criterion to define the weighted parameter, , is given by the ratio of the minimum (resolved) interface ) over the actual interface calculated from Σ & equation. It can be seen on figure 3-a, the higher the interface fluctuations, the lower 2 , which means instabilities have turned into three dimensional effect and secondary atomization zone is considered. The second criterion is grid-dependant, especially if LES turbulence modelling is applied using as a cutoff width ∆, cubit root of the grid cell volume.

; , is defined taking into account the curvature of the interface, k. Figure 3-b displays the relationship between interface radius and cell size. The higher the interface radius, the better resolution of the interface, thus is set to zero. (curvature resolution). Source: [START_REF] Hecht | Large Eddy Simulation of liquid flow: Spray application[END_REF] In this work, the two phases are assumed at constant density, <, so that turbulent fluctuations of mixture density < = are only due to volume fraction, as displayed in equation 4.

< = = < ( = + < > (1 -( = ) (4)
Finally the filtered unsteady Navier-Stokes equations are presented. Instead of time-averaging, LES uses a spatial filtering operation to separate the larger and smaller eddies [START_REF] Versteeg | An Introduction to Computational Fluid Dynamics, The finite volume method[END_REF]. In equation 6, the last term on RHS are the subgrid scale stresses, just like the Reynolds stresses in the RANS momentum equations that arise as a consequence of time-averaging.

? & @ * + ? & @ = A̅ + ,C @ -+ . D , E FG = H I -H & I & = -C > , @ + D -+ ( J E FF . K FG (SGS Stresses) (5)
These stresses are modelled using Smagorinsky-Lilly SGS model which are taken to be proportional to the local rate of strain of the resolved flow. The Smagorinsky-Lilly model builds on Prandtl's mixing length model, and assumes that we can define a dynamic SGS viscosity, C > , which can be described in terms of one length scale ∆, and one velocity scale and it is related to the kinematic SGS viscosity L > . As in the length mixing model, the velocity scale is expressed as the product of the length scale and the average strain rate of the resolved flow, ∆. | ̅ | [START_REF] Versteeg | An Introduction to Computational Fluid Dynamics, The finite volume method[END_REF], Thus, the SGS viscosity is evaluated as follows:

a) b) C > = <( > ∆) 5 | ̅ | (6)
It is worth mention that the previous variables such as liquid volume fraction, , and the mean interface area per unit of volume Σ & , were filtered as well, and the previous approach was implemented.

Results and discussion

Validation

Before proceeding and applying LES-ELSA-Lagrange model to the FCC oil injector is necessary first to validate the algorithm by means of the experimental dataset developed by the ECN Group (Engine Combustion Network). Specifications of the injector are detailed in figure 4-a. In this work the injector chamber was chosen and meshed using OpenFOAM®'s utility, BlockMesh. The mesh consists of 1.4 million un-structured hexahedral cells, as shown in figure 4-b and figure 4-c. Table 1 displays the liquid properties used to simulate the atomization process within the ECN injector chamber.

In order to reach convergence and hence, get time-independent results, the liquid was allowed to cross four times the injector chamber completely in the axial direction. Once the results reach convergence, the liquid volume fraction is plotted on plane slice and as iso-surface coloured by velocity as displayed in figure 5-a, and figure 5-b, respectively. It can be seen on both pictures the difference between the resolved scale near the exit of the injector (inlet of the injector chamber) and the poorlyresolved scales downstream of the flow where the subgrid fluctuations become important. At the inlet there is a clear definition of the interface, thanks to the interface capturing method. On the other hand, downstream of the flow, the turbulent diffusion liquid flux dominates, and it is perceived the beginning of the mixing of the phases as the liquid volume fraction decreases. For validation purpose, two radial planes were located at 500 and 1500 microns from the inlet of the injector tip as displayed in figure 5-c, in order to calculate the surface liquid velocity, ( F | 2 ), as it moves away radially from the centerline of the injector. Firstly the difference between DNS and New DNS is based on how the convective term is solved in the volume fraction and momentum equations. Since the latter relies on the mixture density (equation 4), the momentum equation can be recast, so the exact convective term can also appear, plus additional terms. This convective term is treated differently for DNS cases. On the contrary, New DNS solves consistently [START_REF] Vaudor | Atomisation assistée par un cisaillement de l'écoulement gazeux Développement et validation[END_REF]. Secondly the experimental results were taken using kind of PTV measurements, based on the structure detected techniques. Those velocities cannot be captured by the PTV instrument within the liquid jet core, as shown in Figure 5-c

Application

As stated above, FCC risers involved complex physics, such atomization, vaporization and chemical reactions. In this work, we study the atomization process within a FCC injector as displayed in figure 7-a. In figure 7-b it can be seen the boundary conditions applied to the FCC injector model, hence making it suitable for ELSA algorithm based on high Reynolds and Weber numbers. The mesh was done with an open-source library for mesh generation implemented within OpenFOAM® framework, named cfMesh. The library supports generation of meshes of arbitrary cell types, and the currently implemented workflows generate Cartesian type of polyhedral in both 2D and 3D space, tetrahedral and arbitrary polyhedral [START_REF]cfMesh, meshing tool[END_REF] as exhibited in figure 8-a. In figure 8-b, it is shown the unstructured hexahedral mesh with inflation layers adjacent to the wall, with a total of 5 million cells.

Firstly to evaluate the mesh and get quickly robust and computationally less expensive results than LES, a timeaveraging of the governing equations of the Eulerian solver in ELSA was used. Figure 9-a shows, the results of gradients in the Stream-wise direction of the axial velocity which reveals wall-bounded effects in the phase mixing zone, with high gradients near the wall. Phenomena never encountered in ECN injectors. Based on the fact, the FCC injector model was meshed using the same height of the first cell adjacent to wall, in figure 9-b exposes the values of Y+, which indicates that highest velocities in the FCC injector are located near the wall, therefore requiring a more accurate turbulence model with a better mesh quality (sensitivity analysis) Figure 10-a illustrates the liquid volume fraction using an improved mesh and ELSA-LES model without the interface capturing method. The two-phase flow within the FCC injector demonstrates that there is a liquid core flow and annulus gas flow compressed by the latter inside the injector. This concentric behavior pushes the gas to wall, hence increasing its stream-wise velocity near the wall. Figure 10-b shows the velocity profile in the transverse cross-section which verifies the highest velocities near the wall. based on experimental dataset and DNS simulation. Results have not yet been conclusive due to certain parameter sensitivity when it comes to creating droplets. Authors strongly advise further research in the spray zone as the one being held in our facilities. Having modelled FCC injectors, it has been demonstrated the wallbounded turbulent effects on the two-phase flow, where the highest velocities are located near the wall, thus increasing the phase mixing.

Figure 1 .

 1 Figure 1. Atomization zones. Source: [10]

  . The turbulent diffusion flux represents the transport of the liquid volume fraction induced by velocity fluctuations as shown in equation 2 by analogy with Fick's law.[START_REF] Andreini | Development of a turbulent liquid flux model for Eulerian-Eulerian multiphase flow simulation[END_REF]

Figure 2 .

 2 Figure 2. Ponderation parameter ( ). Source [10]

(

Figure 3 .

 3 Figure 3. a) Interface resolution quality, Σ & . b) Interface resolution quality,; (curvature resolution). Source:[START_REF] Hecht | Large Eddy Simulation of liquid flow: Spray application[END_REF]

Figure 4 .

 4 Figure 4. a) ECN injector geometry and dimensions. b) Mesh view at inlet plane of the injector. c) 3-D Mesh of the injection chamber. Source: [10]

Figure 5 .Figure 6

 56 Figure 5. a) Liquid volume fraction. b) Iso-Surface of 0.5 liquid volume fraction coloured by Velocity. c) Radial probe planes located in longitudinal axis (source [17]).

  Figure 6 exhibits the comparison of the averaged surface liquid velocity, ( F | 2 ), on the two radial planes highlighted above, from Experiments, DNS, New DNS, and ELSA-LES + Lagrange equation for droplets velocity.Firstly the difference between DNS and New DNS is based on how the convective term is solved in the volume fraction and momentum equations. Since the latter relies on the mixture density (equation 4), the momentum equation can be recast, so the exact convective term can also appear, plus additional terms. This convective term is treated differently for DNS cases. On the contrary, New DNS solves consistently[START_REF] Vaudor | Atomisation assistée par un cisaillement de l'écoulement gazeux Développement et validation[END_REF]. Secondly the experimental results were taken using kind of PTV measurements, based on the structure detected techniques. Those velocities cannot be captured by the PTV instrument within the liquid jet core, as shown in Figure5-c, and figures 6-a/6-b. Consequently in the absence of experimental dataset from the centerline of the liquid jet until 100-120 µm, New DNS results were used. Figure6displays good agreement between New DNS and Experiments from 100 to 200 µm (figure6-a) and from 200 to 400 µm (figure6-b), which set the limits of New DNS validated results. Comparison above upper limits fails to predict dispersed zone flow field mainly due to unsufficient mesh refinement using New DNS. Therefore the available validation range is increased from the centerline up to 700 µm approximately by combining both Experimental and New DNS results. Quantitatively speaking on one hand, figure6-a clearly exhibits the differences from the centerline of the liquid jet up to 100 µm radially, between the results from ELSA's and New DNS. On the other hand, beyond 100 µm (200 µm, figure6-b), ELSA tends to underestimate the liquid velocity as compared to New DNS and Experimental results. Such logical results arrive from the formation of tiny droplets beyond 100-200 µm, thus un-resolved scales become relevant as we move outward radially. ELSA-LES + Lagrange model accounts properly for the droplets formation in the dispersed spray zone until certain extent. However, the results presented so far are not conclusive, and they need further research. It is believed when droplets are created, there is a relative velocity between the phases at the interface, instead of an independent flow field for each phase.

Figure 6 .

 6 Figure 6. Left, liquid velocity located 250 microns from the inlet. Right, liquid velocity located 1500 microns from the inlet.

Figure 7 .Figure 8 .

 78 Figure 7. a) FCC injector geometry. b) FCC injector model and boundary conditions.

Figure 9 .

 9 Figure 9. a) Gradients in the Stream-wise direction of the axial velocity. b) yPlus.

Figure 10 .

 10 Figure 10. a) Liquid Volume Fraction. b) Velocity profile in the 'A' transverse cross-section.ConclusionsAn enhanced ELSA-LES model implemented within OpenFOAM® framework was presented. The liquid interface capturing method has been handled a long with the turbulent diffusion liquid flux by a pondered parameter, , which depends on the interface resolution and curvature of the interface. Moreover, a Lagrangian approach was studied to account for the spray atomization zone, where droplets are created. Validation has been carried out

Table 1 .

 1 Liquid properties. Source[START_REF] Hecht | Large Eddy Simulation of liquid flow: Spray application[END_REF]

		Symbols	Value
	Liquid Density Liquid Viscosity Surface Tension	< (NO. ' J ) C (NO. ' ( P ( ) σ (V. ' ( )	803 3,2x10 UJ 2,54x 10 U5

Acknowledgements

This work was partly supported by VINCI TECHNOLOGIES and CNRS. Simulations were carried out at the CRIHAN (Centre de Ressources Informatiques de Haute Normandie) and CINES (Centre Informatique National de l'Enseignement Supérieur).

Nomenclature