Olivier Cotto
email: olivier.cotto2@gmail.com

Luis Miguel Chevin
email: luis-miguel.chevin@cefe.cnrs.fr

FLUCTUATIONS IN LIFETIME SELECTION IN AN AUTOCORRELATED ENVIRONMENT

Most natural environments vary stochastically and are temporally autocorrelated. Previous theory investigating the effects of environmental autocorrelation on evolution mostly assumed that total fitness resulted from a single selection episode. Yet organisms are likely to experience selection repeatedly along their life, in response to possibly different environmental states. We model the evolution of a quantitative trait in organisms with non-overlapping generations undergoing several episodes of selection in a randomly fluctuating and autocorrelated environment. We show that the evolutionary dynamics depends not directly on fluctuations of the environment, but instead on those of an effective phenotypic optimum that integrates the effects of all selection episodes within each generation. The variance and autocorrelation of the integrated optimum shape the variance and predictability of selection, with substantial qualitative and quantitative deviations from previous predictions considering a single selection episode per generation. We also investigate the consequence of multiple selection episodes per generation on population load. In particular, we identify a new load resulting from within-generation fluctuating selection, generating the death of individuals without significance for the evolutionary dynamics. Our study emphasizes how taking into account fluctuating selection within lifetime unravels new properties of evolutionary dynamics, with crucial implications notably with respect to responses to global changes. 2 𝑖 ) -1 is the total strength of selection (that is the harmonic mean of the width of the selection function at each selection episode), 𝜃 𝑡𝑜𝑡 = ∑ 𝑒 𝑖 𝑖 𝜃 𝑖 is the per-generation effective optimum corresponding to the sum of the phenotypic optima at each selection episode, weighted by the contribution of these episodes to the total strength of selection. We denote this contribution as 𝑒 𝑖 = 𝜔 𝑡𝑜𝑡 2 𝜔 𝑖 2 , by reference to elasticities in stage-structured population models (Caswell 2001). The maximum absolute fitness over a generation is 𝑊 𝑚𝑎𝑥,𝑡𝑜𝑡 = 𝐵 exp (-1 2𝜔 𝑡𝑜𝑡 2 (∑ 𝑒 𝑖 𝑖 𝜃 𝑖 2 -(∑ 𝑒 𝑖 𝑖 𝜃 𝑖 ) 2 )) ∏ 𝑊 𝑚𝑎𝑥,𝑖 𝑖

INTRODUCTION

The study of adaptation to changing environments has received renewed attention in the context of current global changes induced by anthropogenic activities. Many theoretical studies on adaptation to changing environments have focused on cases where environmental change follows a deterministic trend [START_REF] Lynch | Adaptive and demographic responses of plankton populations to environmental-change[END_REF][START_REF] Lynch | Evolution and extinction in response to environmental change[END_REF][START_REF] Gomulkiewicz | Demographic and Genetic Constraints on Evolution[END_REF][START_REF] Cotto | Maladaptation as a source of senescence in habitats variable in space and time[END_REF], motivated by major tendencies such as climate warming. Yet, most variation on short timescales results from stochastic fluctuations, and climatic time series generally include stochastic variations around their main trends [START_REF] Stocker | Climate change 2013: The physical science basis[END_REF]. The variance and autocorrelation of these fluctuations need to be modeled to accurately predict the impact of climate warming on biosystems (e.g. [START_REF] Katz | Use of conditional stochastic models to generate climate change scenarios[END_REF][START_REF] Rowell | A scenario of European climate change for the late twenty-first century: seasonal means and interannual variability[END_REF]. Global climatic dynamics have recently been suggested to affect local selection pressures [START_REF] Siepielski | Precipitation drives global variation in natural selection[END_REF], and time series of climatic variables are generally autocorrelated in time, notably as a result of thermic inertia [START_REF] Hasselmann | Stochastic climate models part I[END_REF][START_REF] Rowell | A scenario of European climate change for the late twenty-first century: seasonal means and interannual variability[END_REF]. As most climatic (and other ecologically relevant) environmental variables are autocorrelated in time [START_REF] Vasseur | The color of environmental noise[END_REF], environmental selective pressures are likely to also be temporally autocorrelated, an aspect that is often neglected in models of adaptation that do include stochasticity in the environment (e.g. [START_REF] Lynch | Evolution and extinction in response to environmental change[END_REF][START_REF] Engen | Evolution of a plastic quantitative trait in an agestructured population in a fluctuating environment[END_REF], Engen et al. 2012).

Theoretical models of adaptation to stochastically changing environments often assume that the fitness of an individual depends on the match between its phenotype and an optimum phenotype influenced by the environment (so-called "moving optimum" models, reviewed by [START_REF] Kopp | Rapid evolution of quantitative traits: theoretical perspectives[END_REF]. When the environment undergoes stationary stochastic fluctuations, the effect of evolutionary responses on long-term fitness and expected population growth depends on the variance and autocorrelation of changes in the optimum phenotype (Charlesworth 1993, Lande and Shannon 1996[START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF]. If the stochastic fluctuations of the optimum are not autocorrelated, then responses to selection in a given generation might increase maladaptation in the next, with detrimental consequences for population growth (Lande and Shannon 1996). In contrast, temporal autocorrelation of the optimum makes the environment more predictable, such that the evolutionary response in a given generation is more likely to be beneficial in the next, and that any factor improving the response to selection leads to a higher expected long-term growth rate of the population (Lande andShannon 1996, Chevin 2013). The autocorrelation of the optimum also affects the variance and overall shape of the probability distribution of population size of an evolving population, with consequences for extinction risk [START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF]. Lastly, temporal autocorrelation of the environment is a major driver of the evolution of plasticity, with higher reaction norm slopes evolving under more predictable (autocorrelated) environments of selection [START_REF] Gavrilets | The genetics of phenotypic plasticity. V. Evolution of reaction norm shape[END_REF][START_REF] Lande | Adaptation to an extraordinary environment by evolution of phenotypic plasticity and genetic assimilation[END_REF][START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF].

On the empirical side, a number of studies undertook to determine the temporal characteristics of selection in natural populations. Temporal variation in selection is observed in classic long-term surveys, like beak shape in Darwin's finches [START_REF] Grant | Unpredictable evolution in a 30-year study of Darwin's finches[END_REF], banding patterns in Cepea snails [START_REF] Cain | Population size and morph frequency in a long-term study of Cepaea nemoralis[END_REF] or spine number in threespined sticklebacks [START_REF] Reimchen | Temporal variation in divergent selection on spine number in threespine stickleback[END_REF]Nosil 2002, reviewed in Bell 2010). Consistent with these observations, studies analyzing several datasets with long-term selection estimates using a common framework also found that selection varies in strength and direction [START_REF] Kingsolver | Phenotypic selection in natural populations: what limits directional selection?[END_REF]Diamond 2011, Siepielski et al. 2011, but see [START_REF] Morrissey | Directional selection in temporally replicated studies is remarkably consistent[END_REF]. Further, the development of a statistical framework that models temporal fluctuations in phenotypic selection as a random process, characterized by the variance and autocorrelation of an optimum phenotype, allowed to show that the optimum laying date in a population of Great tits undergoes autocorrelated temporal fluctuations [START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF][START_REF] Gamelon | Environmental drivers of varying selective optima in a small passerine: A multivariate, multiepisodic approach[END_REF]). This framework rests on a model of a Gaussian fitness peak, as in many theoretical studies [START_REF] Kopp | Rapid evolution of quantitative traits: theoretical perspectives[END_REF], so that the estimated parameters should have a direct theoretical interpretation.

However, the comparison of empirical measures of fluctuating selection with theoretical results on adaptation to autocorrelated environment is an arduous task, because the definitions of fitness and selection may not be easily transposed from one to the other. Even regardless of theoretical considerations about the definition of long-term expected fitness in a stochastic environment, understanding how variation in selection integrates within lifetime is not straightforward. Most theoretical studies that included autocorrelation in the environment neglected the dynamics of selection within generations, by assuming (often implicitly) that selection occurs instantaneously, in discrete non-overlapping generations. In other word, lifetime fitness is assumed to result from a single instantaneous episode of selection; or at least, the unfolding of selection within lifetime is not modeled explicitly, making it ambiguous whether the per-generation optimum represents instantaneous or repeated selection. More realistically, selection occurs in several episodes along a generation (Lande 1982, Arnold andWade 1984b). In fact, the vast majority of long-term estimates of selection gradients focused on a single fitness component (~ 98% among 2819 selection estimates in [START_REF] Kingsolver | Phenotypic selection in natural populations: what limits directional selection?[END_REF], and therefore correspond to a single episode of selection, among several possible others in a generation or lifetime. This has several related consequences for the interpretation of measurements of selection. First, the strength of selection over a single episode cannot necessarily be compared to predictions based on overall selection per generation, because this selection episode needs to be weighted by its contribution to the lifetime (or per-generation / per-unit-time) fitness. Second, both the magnitude and predictability of temporal variation in selection is probably not well captured by treating each selection episode in isolation. The reason is that in organisms that undergo multiple selection episodes, environmental variability gets integrated over lifetime. This has important consequences, notably regarding the contribution of extreme events to selection, a topic of recent interest [START_REF] Grant | Evolution caused by extreme events[END_REF][START_REF] Marrot | Multiple extreme climatic events strengthen selection for earlier breeding in a wild passerine[END_REF].

The study of demography has long questioned the relationship between variation in vital rates that contribute to elements of the transition matrix (e.g. fecundities, survival rates) and various measures of the population growth rate [START_REF] Tuljapurkar | The many growth rates and elasticities of populations in random environments[END_REF][START_REF] Morris | Buffering of life histories against environmental stochasticity: accounting for a spurious correlation between the variabilities of vital rates and their contributions to fitness[END_REF][START_REF] Engen | Estimating the effect of temporally autocorrelated environments on the demography of density-independent age-structured populations[END_REF]. This theory has been used to analyse fluctuating selection in stochastic environments that are not autocorrelated (Engen et al. 2012[START_REF] Engen | Estimating phenotypic selection in age-structured populations by removing transient fluctuations[END_REF]. However, how stochastic autocorrelated fluctuating selection mediated by different fitness components integrates over lifetime is an underinvestigated topic, limiting our understanding of adaptation to randomly changing environments.

Making progress on this question requires extending previous approaches that integrate the multiplicity of selection episodes along life. Of particular note is the approach by Arnold andWade (1984b, 1984a), who proposed and applied a decomposition of lifetime selection into a series of selection episodes with multiplicative effects on fitness, resulting in selection gradients (dependent on log-fitness) that sum along a generation (see also [START_REF] Mcglothlin | Combining selective episodes to estimate lifetime nonlinear selection[END_REF]. Integrating the effect of successive selection episodes on lifetime fitness is appropriate to organisms with non-overlapping generations (Arnold and Wade 1984a), but has not been used to investigate the effect of environmental stochasticity.

We here build a model of a phenotypic trait exposed to multiple selection episodes, in nonoverlapping generations. We assume that the phenotype of an individual at the focal trait does not change between selection episodes. This could correspond for example to a morphological trait that only influences selection after it is fully developed, or a coloration trait that affects survival probability (due to e.g. predation, or thermoregulation) prior to reproduction. Under these assumptions, we detail and clarify how the temporal structure of selection during lifetime influences the predictability of lifetime viability selection, and the load caused by maladaptation in a stochastic environment. Our model applies directly to organisms that are semelparous with non-overlapping generations, including univoltine insects such as cicadas, and species with economic or cultural significance such as salmons or bamboos. Our analysis further applies to cohort analysis (see discussion) when considering more general life histories with overlapping generations. As such, our study provides a necessary step toward a complete understanding of how life history affects evolutionary dynamics in the context of autocorrelated fluctuating selection.

Overall, we find that the way selection operates during a lifetime strongly determines the impact of environmental fluctuations on variation in selection across generations, and the incurred maladaptation load in a stochastic environment.

MODEL

Selection model

We investigate the evolution of a single quantitative trait z in an organism with non-overlapping generations. Along their life cycle, individuals undergo several selection episodes that determine their probability to survive between stages before reproduction, with the latter occurring only at the end of each generation. We further assume that there is no selection on fecundity, which we denote B. For simplicity, we focused on a single trait, but our framework can be readily extended to the multivariate case, with different (correlated) traits undergoing selection at different stages (Arnold and Wade 1984b, Arnold and Wade 1984a[START_REF] Cotto | Maladaptation as a source of senescence in habitats variable in space and time[END_REF].

We further assume that the phenotype z at the beginning of a generation can be decomposed as z = x + e, where x is the additive genetic value and e is a residual component of variation independent of the additive genetic value, and with normal distribution of mean 0 and variance E [START_REF] Falconer | Introduction to quantitative genetics[END_REF]. We assume that many loci contribute to the additive genetic value, such that x is normally distributed with mean 𝑥̅ and variance G. Under these assumptions, z is also normally distributed with mean 𝑧̅ and variance P = G + E.

The survival probability of an individual with phenotype z at age i in generation g depends on the match between its phenotype z and an optimum phenotype 𝜃 𝑔,𝑖 that changes with the environment, 𝑊 𝑖,𝑔 (𝑧) = 𝑊 𝑚𝑎𝑥,𝑖 exp(-

(𝑧-𝜃 𝑖,𝑔 ) 2 2𝜔 𝑖,𝑔 2 ).
(

) 1
where 𝑊 𝑚𝑎𝑥,𝑖 is the maximum survival rate at age i and 𝜔 𝑖 measures the width of the Gaussian survival function (inversely related to the strength of selection). We assume that selection occurs through all transitions preceding reproduction, so that the components of fitness are multiplicative (Arnold and Wade 1984a). If reproduction could occur at several stages, then the components of fitness would no longer be multiplicative, and another framework (e.g. overlapping generations, [START_REF] Charlesworth | Evolution in age-structured populations[END_REF]) should be used. The total absolute fitness of an individual with phenotype z is then (where the subscript g has been dropped for simplicity)

𝑊(𝑧) = 𝐵 ∏ 𝑊 𝑖 (𝑧) 𝑛 𝑖=1
(2)

Under the assumption that selection is Gaussian at each episode (eq. 1), the total fitness function is also Gaussian,

𝑊(𝑧) = 𝑊 𝑚𝑎𝑥,𝑡𝑜𝑡 exp (- (𝑧-𝜃 𝑡𝑜𝑡 ) 2 2𝜔 𝑡𝑜𝑡 2 ), (3)
. ( 4)

With the further assumption that the strength of selection (here measured as a width of the fitness function) is the same at all selection episodes, we have 𝜔 𝑖 = 𝜔 and 𝑒 𝑖 = 1 𝑛 ⁄ , so the equations above further simplify as

𝜔 𝑡𝑜𝑡 2 = 𝜔 2 /𝑛 (5)
𝜃 𝑡𝑜𝑡 = 1 𝑛 ∑ 𝜃 𝑘 𝑘 = E ̂[𝜃] (6)
𝑊 𝑚𝑎𝑥,𝑡𝑜𝑡 = 𝐵 exp (-

𝑛V ̂[𝜃] 2𝜔 2 ) ∏ 𝑊 𝑚𝑎𝑥,𝑖 𝑖 , (7)
where E ̂[𝜃] andV ̂[𝜃] are the "sample" expectation and variance (respectively) of θ within a generation.

The change in the mean phenotype over an entire generation, after reproduction, is [START_REF] Lande | Natural selection and random genetic drift in phenotypic evolution[END_REF])

∆𝑧̅ = ℎ 0 2 (𝑧̅ 𝑛 -𝑧̅ 0 ) = 𝐺 0 𝜕 ln 𝑊 ̅ (𝑧) 𝜕𝑧̅ , ( 8
)
where 𝑧̅ 𝑛 is the mean phenotype after the last episode of selection, and additive genetic variance 𝐺 0 and heritability ℎ 0 2 = 𝐺 0 /𝑃 0 are measured at the beginning of a generation, prior to any selection. Indeed, each selection event can generate a covariance between P and E that needs to be accounted for if dealing with selection episode by episode (as in Arnold and Wade 1984a), resulting in complications (see Supp. Mat.) that can be overcome by focusing on the effect on total selection on the phenotype distribution prior to any selection. With Gaussian stabilizing selection (eq. 1) and a normally distributed trait, the mean fitness in the population is 𝑊 ̅ (𝑧̅ ) = ∫ 𝑝(𝑧)𝑊(𝑧)𝑑𝑧 = 𝑊 𝑚𝑎𝑥,𝑡𝑜𝑡 √𝑆 𝑡𝑜𝑡 𝜔 𝑡𝑜𝑡 2 exp (-

(𝑧̅ -𝜃 𝑡𝑜𝑡 ) 2 2(𝜔 𝑡𝑜𝑡 2 +𝑃 0 ) ), (9)
where

𝑆 𝑡𝑜𝑡 = 1 𝑃 0 +𝜔 𝑡𝑜𝑡 2 .
From equation 8, it follows that the per-generation change in the mean trait value in the population is

∆𝑧̅ = -𝐺 0 𝑆 𝑡𝑜𝑡 (𝑧̅ -𝜃 𝑡𝑜𝑡 ) = -𝐺 0 𝛽 𝑡𝑜𝑡 , ( 10
)
where 𝛽 𝑡𝑜𝑡 is the total selection gradient, and which has a similar form as in models with a single selection episode (e.g. [START_REF] Lande | Natural selection and random genetic drift in phenotypic evolution[END_REF]).

Fluctuations of the environment and timing of selection

We are interested in cases where the optimum phenotype fluctuates because of random fluctuations in the environment, i.e. fluctuations of the optimum phenotype represent those of the environment.

For simplicity, we assume that the optimum phenotypes at all selection events in a life cycle respond to the same environmental variable, so that fluctuations in the optimum during a lifetime can be modeled as a temporal sampling of the same stochastic process. We consider the case where the optimum follows a first-order autoregressive process with mean 0 (without loss of generality), defined by

𝜃 𝑡+1 = 𝜌 𝜃 𝑡 + 𝜎 𝜃 √1 -𝜌 2 𝜁 𝑡 ( 11
)
where 𝜌 is the autocorrelation of the optimum over one (absolute; e.g. one year) time unit (which we assume is positive), 𝜎 𝜃 2 its stationary stochastic variance, and 𝜁 𝑡 a standard normal random deviate. For such a process, the autocorrelation of the optimum over k units of time is 𝜌 𝑘 = exp(-𝑘/𝑇 𝜃 ), where 𝑇 𝜃 = -1/ln(𝜌 ) is the characteristic timescale of the autocorrelation of the optimum. Note that the autoregressive process converges to white noise when 𝜌 tends to 0, i.e.

when the optimum phenotype is correlated over an infinitesimally small timescale.

We assume that n episodes of selection occur before reproduction. These selection episodes all take place over a time-window that spans a fraction α of generation time T. The time-window for selection remains unchanged across generations. This models the fact that all traits in an organism need not be exposed to selection throughout life. Instead, it is more likely that there is a time-window within a generation during which a given trait may be sensitive to the environment and exposed to selection (e.g. during the first stages). When exactly the time-window of selection occurs within a generation (e.g. at the beginning or at the end) has no effect for the remaining of the analysis. Under our assumption of non-overlapping generations, only the duration of this window matters. For simplicity, we assume that the episodes of selection are evenly spaced along the time-window for selection, so that they occur every αT/n units of time. Using the index g for the number of generations, we denote as 𝑧̅ 𝑔,0 the mean trait value before the first episode of selection of generation g. The timescale of life history relative to that of the environment is represented in figure 1.

RESULTS

Predictability of selection

The predictability of selection can be defined in several ways. For our purpose, we will focus on the autocorrelation of the effective optimum per generation, and that of the total selection gradient, because these relate to measurable quantities (Lande and Arnold 1983[START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF][START_REF] Gamelon | Accounting for interspecific competition and age structure in demographic analyses of density dependence improves predictions of fluctuations in population size[END_REF]) known to affect adaptation and population dynamics in a stochastic environment (Lande and Shannon 1996[START_REF] Chevin | Genetic constraints on adaptation to a changing environment[END_REF][START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF]. We first calculate the autocorrelation of 𝜃 𝑡𝑜𝑡 , before investigating the properties of the fluctuations in the selection gradient 𝛽 𝑡𝑜𝑡 , which despite being less directly connected to theoretical predictions, are closer to what is commonly estimated empirically (Lande and Arnold 1983).

Autocorrelation of the effective optimum.

The autocorrelation of the effective optimum at a lag of one generation is (see also Supp. Mat.)

𝜌 𝜃𝑡𝑜𝑡 = 𝜌 𝑇 ∑ 𝑒 𝑘 2 𝑛 𝑘=1 +∑ ∑ 𝑒 𝑖 𝑒 𝑗 𝜌 (𝑗-𝑖)𝛼𝑇 𝑛 ⁄ (1+𝜌 2(𝑖-𝑗)𝛼𝑇 𝑛 ⁄ ) 𝑛 𝑗=𝑖+1 𝑛-1 𝑖=1 ∑ 𝑒 𝑘 2 +2 ∑ ∑ 𝑒 𝑖 𝑒 𝑗 𝜌 (𝑗-𝑖)𝛼𝑇 𝑛 ⁄ 𝑛 𝑗=𝑖+1 𝑛-1 𝑖=1 𝑛 𝑘=1 = 𝜌 𝑇 𝜑, (12)
where 𝜌 𝑇 is the autocorrelation of the optimum over a generation time, and its autocorrelation between two consecutive selection episodes is 𝜌 𝛼𝑇/𝑛 (which differs fromthe autocorrelation 𝜌 per time-unit). In the following we detail the inflating factor φ corresponding to the complex fraction that appears in the first equality in equation ( 12). The first sum in the numerator results from the covariance between the optimum in a given episode of selection and the same episode in the next generation. The second sum results from the covariances of the optimum between two different selection episodes in two consecutive generations. The denominator has a similar structure, and results from similar effects (variance per episode vs covariance across episodes) within a generation, contributing to the variance of the effective optimum (see Supp. Mat.). It can be shown (by noting that 0 < 𝜌 (𝑗-𝑖)𝛼𝑇 𝑛 ⁄ ≤ 1, so that 𝜌 -(𝑗-𝑖)𝛼𝑇 𝑛 ⁄ > 1) that the ratio in equation ( 12) is always larger than 1, and is thus an inflating factor, which always increases the autocorrelation of the effective optimum relative to the autocorrelation of the environment, as a results of the averaging of selection between episodes within generations. Remarkably, this result extends to the autocorrelation over 𝜏 > 1 generations when the autocorrelation of the environment decreases exponentially with time, as occurs under a Markovian process such as our first-order autoregressive process. We thus have for the autocorrelation of the effective optimum over 𝜏 ≥ 1 generations 𝜌 𝜃𝑡𝑜𝑡,𝜏 = 𝜌 𝜏𝑇 𝜑.

(13) Equation ( 13) shows that, once the inflation factor has been accounted for, the autocorrelation of 𝜃 𝑡𝑜𝑡 declines exponentially at the same rate as the autocorrelation of θ (fig. 2A). In principle, then, information on the autocorrelation of the effective optimum between two generations allows to determine its autocorrelation over many generations, providing that the autocorrelation of the environment (here also corresponding to the autocorrelation of θ) driving selection is known.

More insights arise from the simpler case where the strength of selection is the same at all selection episodes (eq. 5-6), wherein 𝜑 simplifies to

𝜑 = 𝜓+𝜓 -1 -2 2(𝜓-1)+𝑛(𝜓 - 1 𝑛-𝜓 1 𝑛) , (14)
where 𝜓 = 𝜌 𝛼𝑇 is the autocorrelation of θ over the time window for selection. Equation ( 14) shows that when variance in the strength of stabilizing selection within generation can be neglected, then the inflation factor does not depend on the strength of selection, but only on environmental autocorrelation 𝜓 over the time window for selection, and on the number of selection episodes.

The inflation factor φ increases when the autocorrelation of the environment decreases (fig. 2B,C), and when the number of selection episode increases (fig. 2B). In other words, weakly autocorrelated environments can actually correspond to much larger autocorrelations of variables relevant to selection. For example, when there are many episodes of selection per generation and the autocorrelation of the environment is 0.2, the autocorrelation of 𝜃 𝑡𝑜𝑡 can be increased by more than 50% as compared to the autocorrelation of θ (fig. 2B). When the number of episodes of selection is large, φ converges to

𝜑 𝑛→∞ → (1-𝜓 ) 2 2𝜓(𝜓-ln 𝜓-1) , ( 15
)
which only depends on environmental autocorrelation 𝜓 over the time window for selection. We found that equation ( 15) provides a good approximation for equation ( 13) as soon as there are more than a few episodes of selection per generation, and even if the actual strength of selection varies across episodes (fig. 2B). Interestingly, equation ( 15) does not depend on the strength of selection at each episode, and proposes a simple transformation between the autocorrelation of the environment and the autocorrelation of the effective optimum over a generation time. Substituting equation ( 15) into equation ( 13) thus provides, in principle, a simple way to estimate the autocorrelation of the effective optimum in natural populations from a measure of the autocorrelation of the environment, over the duration of selection within a generation on one hand (to compute the inflation factor), and across generations on the other hand. The longer the organism is sensitive to selection within a generation, the larger the inflation factor (fig. 2C). Lastly, when the environment is not autocorrelated, the effective optimum is not autocorrelated either, showing that within-generation selection alone does not generate autocorrelation of the effective optimum.

Fluctuations in directional selection.

Variance of the selection gradient

Ultimately, evolutionary change in the mean phenotype in response to selection is directly related to the total per-generation selection gradient 𝛽 𝑡𝑜𝑡 (eq. 10). At stationarity after many generations, we find that the variance of the selection gradient is approximately (see Supp. Mat.)

Var[𝛽 𝑡𝑜𝑡 ] ≈ 2Var[𝜃 𝑡𝑜𝑡 ]𝑆 𝑡𝑜𝑡 2 1-𝜌 𝑇 (1+𝑆 𝑡𝑜𝑡 𝐺 0 (𝐼 𝑓 -1)) (2-𝑆 𝑡𝑜𝑡 𝐺 0 )(1-𝜌 𝑇 (1-𝑆 𝑡𝑜𝑡 𝐺 0 )) . ( 16
)
The variance of the selection gradient increases with the variance of fluctuations in the effective optimum (and of the environment, see Supp. Mat.) and with the total strength of selection (eq. 16 and fig. 3). Equation ( 16) further shows that the variance of the selection gradient decreases when the inflation factor increases. Correspondingly, Var[𝛽 𝑡𝑜𝑡 ] decreases with α (because the inflation factor increases with α, fig. 2C) and with the number of selection episodes (fig. 2B and3A). In other words, recurrent selection along a generation tends to decrease the variance of the selection gradient, as expected since environmental variation is averaged among selection episodes.

Less expected is the effect of environmental autocorrelation. We find that with several selection episodes per generation, the variance of the selection gradient depends non-monotonically on the autocorrelation of the environment, with a maximum when the environment is moderately autocorrelated (fig. 3B). This result contrasts with the case with a single selection episode per generation, where the variance in the selection gradient decreases with increasing autocorrelation of the environment, as a result of better adaptive tracking of the optimum (Lande andShannon 1996, Chevin andHaller 2014). This difference arises because, with several episodes of selection, positive autocorrelation in the environment not only produces autocorrelation of the total optimum 𝜃 𝑡𝑜𝑡 across generations, but also increases the variance of 𝜃 𝑡𝑜𝑡 (see eq. A1 and fig. A1 in Supp.

Mat.), with antagonistic effects on adaptive tracking by genetic evolution, and hence on the variance of the total selection gradient across generations (Lande andShannon 1996, Chevin andHaller 2014).

We also investigated the influence of variable selection strength, by comparing the variance of the selection gradient obtained for a fixed strength of selection per episode (eq. 5-6) with the case where the strength of selection is drawn randomly at each episode from a distribution with the same expected total selection Stot (details in fig. 2

Autocorrelation of the selection gradient

For the autocorrelation function of the total selection gradient across generations, we start by analyzing its expression in the case of a single selection episode per generation (i.e. when φ = 1).

We obtain (see Supp. Mat.) for selection gradients 𝜏 generations apart

𝜌 𝛽𝑡𝑜𝑡,𝜏|𝑛=1 = 𝜌 𝜏𝑇 -𝜅(1-𝑆 𝑡𝑜𝑡 𝐺 0 ) 𝜏 1-𝜅 , ( 17
)
where 𝜅 ≈ 𝑆 𝑡𝑜𝑡 𝐺 0

(1+𝜌 𝑇 ) 2(1-𝜌 𝑇 )
(see eq. A15). Equation ( 17) is similar in form to an expression obtained previously using a continuous-time approximation (eq. 5a in Chevin and Haller 2014), and leads to similar conclusions. The autocorrelation of the selection gradient with one selection episode is a linear combination of two exponentially decreasing functions, corresponding on the one hand to the autocorrelation of the environment (first term in numerator), and to evolutionary responses to deviations from the optimum (second term in the numerator) on the other hand. Consistent with the result of Chevin and Haller (2014), we find that for large autocorrelation in the environment, equation ( 17) converges to (1 -𝑆 𝑡𝑜𝑡 𝐺 0 ) 𝜏 , demonstrating that inertia in responses to selection sets an upper limit to the autocorrelation of the selection gradient.

In the general case where selection can occur in several episodes within a generation, we find a more complex expression for the autocorrelation of the selection gradient, which can nevertheless still be written in the form

𝜌 𝛽𝑡𝑜𝑡,𝜏 = 𝜑𝜀𝜌 𝜏𝑇 +𝛾(1-𝑆 𝑡𝑜𝑡 𝐺 0 ) 𝜏 𝜗 , ( 18
)
where the detailed expressions of ε, γ and ϑ depend on the per-generation evolutionary potential 𝑆 𝑡𝑜𝑡 𝐺 0 , the autocorrelation of the environment over a generation 𝜌 𝑇 , and γ and ϑ depend additionally on the inflation factor φ. The details of these terms can be found in the calculations provided in supplementary material (eq. A14). Equations ( 17) and ( 18) provide general expressions for the autocorrelation of the selection gradient with discrete generations and multiple selection episodes, extending previous results that used a continuous-time approximation with a single episode (Chevin and Haller 2014).

Comparing equations ( 17) and ( 18 2B and eq. 18). We also find that for a given number of selection episodes per generation, the autocorrelation of the selection gradient slightly decreases when the total strength of selection increases (fig. 4B).

Over time intervals longer than a generation, the autocorrelation function of selection gradients for a single selection episode (eq. 17) and for multiple episodes (eq. 18) have a similar shape (fig. 5), because they depend on the autocorrelation of the environment and of the effective optimum (respectively), which have the same exponential rate of decrease (fig. 2A). Importantly however, the autocorrelation of 𝛽𝑡𝑜𝑡 in the first few generations -where the effect of multiple episodes is most pronounced -weights most in determining the expected maladaptation in a given generation (see e.g. eq. A6 in Supp. Mat.), so that differences in autocorrelation over short intervals of generations can have profound effects. It is worth mentioning that equations ( 17) and ( 18) can produce negative autocorrelations (fig. 5) when the effect of the evolutionary inertia (second term at the numerator) grows larger that the effect of the autocorrelation of environment (see also Chevin and Haller 2014).

Lastly, and similarly to the variance of the selection gradient, we found that when the strength of selection is drawn randomly at each selection episode (see previous related paragraph for variance), the autocorrelation of selection is below its expectation under constant selection strength, when there are only few selection episodes per generation (fig. 4A). This difference is especially marked when the environment is very autocorrelated and when selection is strong (fig.

4A, B) and vanishes as the number of selection episodes per generation increases.

Maladaptation load

Much of the work on adaptation of quantitative traits to autocorrelated stochastic environments has focused on the resulting expected lag load caused by deviations of the mean phenotype from the optimum, as this directly connects to long-term population growth and extinction risk (Lande andShannon 1996, Chevin et al. 2017). From equation ( 10), the log mean fitness in the population (which relates to the relative increase in population size) is

ln𝑊 ̅ = 𝑟 𝑚𝑎𝑥 - 1 2𝜔 𝑡𝑜𝑡 2 [∑ 𝑒 𝑘 𝑘 𝜃 𝑘 2 -(∑ 𝑒 𝑘 𝑘 𝜃 𝑘 ) 2 ] - 1 2 ln (1 + 𝑃 0 𝜔 𝑡𝑜𝑡 2 ) - 𝑆 𝑡𝑜𝑡 2 (𝑧̅ -𝜃 𝑡𝑜𝑡 ) 2 (19)
where

𝑟 𝑚𝑎𝑥 = ln𝐵 + ∑ 𝑟 𝑚𝑎𝑥,𝑘 𝑛 𝑘=1
corresponds to the maximum potential fitness, for a hypothetical population that would be monomorphic for the optimal phenotype in each selection episode. The last two terms describe reductions in fitness caused by maladaptation of the mean phenotype and phenotypic variance in the population, as described in the literature (Lande and Shannon 1996).

When considering several selection episodes per generation, there is an additional load that is directly linked to the fluctuations in the optimum (second term in eq. 19). Indeed, this load vanishes when the optimum does not fluctuate across episodes (𝜃 𝑘 = 𝜃 for all k). It thus corresponds to the deaths of individuals induced by environmental fluctuations within a generation, which reduce population growth without evolutionary significance. Note that when the trait has phenotypic variance, the mean phenotype does change under selection within a generation, tracking the optimum phenotype to some extent following each selection episode. But each of these selection episodes requires a number of selective deaths, and the resulting total cost of natural selection (sensu Haldane 1957) over a generation is unaffected by phenotypic tracking: the load is exactly the same if there is no phenotypic variation, and no within-generation phenotypic tracking. When the strength of selection at each episode ω is constant (corresponding to eq. 5), we further have

1 2𝜔 𝑡𝑜𝑡 2 [∑ 𝑒 𝑘 𝑘 𝜃 𝑘 2 -(∑ 𝑒 𝑘 𝑘 𝜃 𝑘 ) 2 ] = 𝑛V ̂[𝜃] 2𝜔 2 , ( 20
)
such that the intra-generational load is proportional to the sample variance V ̂[𝜃] of the optimum across selection episodes within a generation.

In models of selection with a Gaussian fitness peak, there is a simple connection between the temporal distribution of the directional selection gradients and the expected maladaptation 6A), contrary to the case of a single selection episode. Furthermore, the expected maladaptation decreases with increasing number of selection episodes within generation (fig. 6B), because this decreases the variance of the effective optimum, and in turn the variance of the selection gradient (see Supp. Mat. and fig. A1). However, note that this decreased variance of the effective optimum come at the expense of an increased load of fluctuating selection within a generation, which cannot be overcome by adaptive tracking.

DISCUSSION

Most theoretical studies investigating the effect of an autocorrelated environment on evolutionary dynamics and their consequences for population dynamics assumed simple life cycles, where the timing of selection within a generation was left implicit (e.g. Charlesworth 1993, Lande and Shannon 1996[START_REF] Chevin | Genetic constraints on adaptation to a changing environment[END_REF]. In models of non-overlapping generations, the timing of life-history events is generally not addressed, and selection is most often modeled as a single, effectively instantaneous process. This contrasts with the fact that many traits are likely to affect fitness at different times in life, and thus undergo several episodes of selection (Arnold andWade 1984b, Marshall and[START_REF] Marshall | Ecological and evolutionary consequences of linked lifehistory stages in the sea[END_REF]. Another difficulty at the interface of theory and empirical work on these questions is that measuring fitness is in practice challenging, such that natural selection is mostly estimated in the form of individual selection episodes acting via specific fitness components in the life cycle [START_REF] Kingsolver | Phenotypic selection in natural populations: what limits directional selection?[END_REF]Diamond 2011, Siepielski et al. 2017) To this aim, we developed a model of an organism with non-overlapping generations, where several transitions between stages (survival) prior to reproduction depend on the match between a phenotypic trait and a fluctuating optimal value directly related to the environment, consistent with recent empirical estimates of fluctuating selection in the form of a moving phenotypic optimum [START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF][START_REF] Gamelon | Accounting for interspecific competition and age structure in demographic analyses of density dependence improves predictions of fluctuations in population size[END_REF]. Our analysis of this model demonstrates that knowledge of fluctuation patterns for the environment affecting selection is not sufficient to infer the predictability of selection across generations, when there are multiple selection episodes per generation. Selection across generations depends on an integrated optimum that includes the relative contribution of each selection episode to total selection along the life cycle, similar to recent formulations for age-structured populations but without autocorrelation [START_REF] Engen | Evolution of a plastic quantitative trait in an agestructured population in a fluctuating environment[END_REF][START_REF] Cotto | Maladaptive Shifts in Life History in a Changing Environment[END_REF]. Therefore, fluctuations of selection across generations depend on the pattern of variance and autocorrelation of the effective optimum, rather than of the environment per se. One of our main results is that, when the environment undergoes red noise (first-order autoregressive process), fluctuating selection within a generation causes the autocorrelation of the effective optimum to be inflated relative to the autocorrelation of the environment. We show that as the number of selection episode per generation increases, the inflation factor rapidly converges to an asymptotic value that depends neither on the number of episodes of selection nor on the strength of selection, but only on the autocorrelation of the environment over a time window of selection. This asymptotic value can be used in practical situations to estimate the autocorrelation of the effective optimum directly from measurement of the autocorrelation of the environment. Importantly, however, this result depends on the assumption that total selection (or its expectation) does not depend on the number of selection episodes per generation. This assumption allows us to disentangle the effect of selection acting repetitively along a generation from the effect of stronger total selection. We emphasize that this independence may not be satisfied in natural populations, where more selection episodes may often cause stronger overall selection.

The inflation factor for the optimum directly influences the autocorrelation of directional selection gradients, which relates more directly to empirical measurements of the predictability and consistency of directional selection through time (Grant andGrant 2002, Morrissey and[START_REF] Morrissey | Directional selection in temporally replicated studies is remarkably consistent[END_REF]. However, an important point when considering the connection with empirical work is that, instead of tracking phenotypic change under selection episode per episode as done empirically (Arnold andWade 1984b, McGlothlin 2010), we combined all episodes to derive a We also found that allowing for multiple selection episodes complicates the connection between environmental autocorrelation and the variance of selection gradients, which is directly linked to the expected maladaptation load reducing population growth [START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF]. With multiple selection episodes, both these metrics vary non-monotonically with the autocorrelation of the environment, with a maximum at an intermediate value. This contrasts with previous results from models with a single episode of selection per generation, where the expected maladaptation load decreases monotonically with increasing autocorrelation of the environment, because this allows for closer adaptive tracking of the moving optimum via genetic evolution [START_REF] Chevin | Genetic constraints on adaptation to a changing environment[END_REF][START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF]. That this may not hold with multiple episodes of selection crucially revises previous understanding of the role of environmental autocorrelation on maladaptation: maladaptation can actually increase with increasing autocorrelation of the environment.

Lastly, we found that selection that fluctuates within a generation results in a specific load, which adds up to the maladaptation and variance loads classically described in the literature on adaptation to changing environments (see [START_REF] Lynch | Evolution and extinction in response to environmental change[END_REF][START_REF] Burger | Evolution and extinction in a changing environment: a quantitative-genetic analysis[END_REF], Lande and Shannon 1996[START_REF] Kopp | Rapid evolution of quantitative traits: theoretical perspectives[END_REF]. This load only exists when the direction and/or strength of selection varies within a generation, and also increases with the number of selection episodes and the total per-generation strength of selection. The reason for this load is that each selection episode results in selective deaths that reduce the demographic output of the population, as previously discussed in extensions of [START_REF] Haldane | The cost of natural selection[END_REF] famous "cost of natural selection" to the context of adaptation to a changing environment, with a single selection episode [START_REF] Nunney | The cost of natural selection revisited[END_REF]. Within-generation fluctuations in selection increase this number of selective deaths without any net beneficial effect on adaptation: from an evolutionary perspective, individuals are essentially wasted by within-generation fluctuation in selection. An important implication of this result is that disturbances that increase temporal variation in the environment can lead to the decline of a population, regardless of its phenotypic state (mean phenotype, additive genetic variance). This bears particular importance in light of the increasing frequency of extreme climatic events [START_REF] Coumou | A decade of weather extremes[END_REF], and accumulating evidence of their effect on evolutionary and demographic dynamics [START_REF] Moreno | Extreme climatic events in relation to global change and their impact on life histories[END_REF]Møller 2011, Marrot et al. 2017). This load caused by fluctuating selection within a generation may only be alleviated by mechanisms that cause withingeneration adaptive changes in the mean phenotype without requiring any selective death. This suggests that a powerful mechanism by which labile, reversible plasticity can evolve (as investigated theoretically in a few studies, [START_REF] Gabriel | How stress selects for reversible phenotypic plasticity[END_REF][START_REF] Lande | Evolution of phenotypic plasticity and environmental tolerance of a labile quantitative character in a fluctuating environment[END_REF][START_REF] Ratikainen | The coevolution of lifespan and reversible plasticity[END_REF] is by reducing the load caused by fluctuating selection within a generation.

Our non-overlapping generation model applies directly to univoltine species. An interesting extension would be to investigate the case with overlapping generations. The most closely related life cycle with overlapping generations is that of perennial semelparous organisms, encompassing a large diversity of organisms (Young and Augspurger 1991) spanning from microbial parasites [START_REF] Ebert | Optimal killing for obligate killers: the evolution of life histories and virulence of semelparous parasites[END_REF] to long lived plants (e.g. semelparous yucca) and vertebrates (e.g. pacific salmon). This type of life cycle can lead to asynchrony in reproduction among age classes (Caswell 2001 p.81-88), such that the population at any time is composed of cohorts with different ages that rarely interbreed, in which case our conclusions directly apply to each of these cohorts. Perennial semelparous life cycles can also produce instability in age structure, which has been proposed to contribute to periodicity in species such as cicadas [START_REF] Bulmer | Periodical insects[END_REF] or bamboos [START_REF] Keeley | Mast flowering and semelparity in bamboos: the bamboo fire cycle hypothesis[END_REF]. This instability may even lead to a collapse of age structure and to a synchrony in reproduction [START_REF] Mjølhus | On synchronization in semelparous populations[END_REF], in which case our model would apply to the entire population. In contrast, iteroparity in perennial organisms would add considerable complications to the investigation of the effects of fluctuations in the environment.

Indeed, each cohort of offspring is contributed from cohorts of reproductive individuals having and variance 4𝜔 𝑡𝑜𝑡 2 -2𝑛. For simulations, we kept the expected total strength of selection constant (independent on the number of selection episode) to isolate the effect of the number of selection episodes from that of the total strength of selection. For A and B, α = 0.8. For all panels: light gray to black: 𝜌 𝑇 = 0.2, 0.5 and 0.9 respectively, 𝜔 ̅ 𝑡𝑜𝑡 2 = 25, The dashed lines correspond to the predictions from eq. ( 16). Continuous lines represent predictions derived from a continuous time approximation as in Lande and Shannon (1996, eq. 7), resulting in equation 5a in Chevin and Haller (2014):

𝜎 𝜃 2 = 1.
Var[𝛽 𝑡𝑜𝑡 ] = Var[𝜃 𝑡𝑜𝑡 ]𝑆 𝑡𝑜𝑡 2 (1 + 𝑆 𝑡𝑜𝑡 𝐺 0 𝑇 𝜃 ) ⁄
, where 𝑇 𝜃 is the characteristic time scale of the autocorrelation of the environment (from eq. 13, and see eq. 11 and below). For both panels: α = 1 and 𝜎 𝜃 2 = 1. single selection episode (eq. 17). From light gray to black, 𝜌 𝑇 = 0.1, 0.5 and 0.9. Panel A: 𝜔 𝑡𝑜𝑡 2 = 10. Panel B: n = 5.

For both panels: α = 1 and 𝜎 𝜃 2 = 1. 18). Continuous lines: expectation for a single selection episode (eq. 17). Note that for parameters shown here (number of selection episode and strength of selection) there is almost no difference gray to black, 𝜌 𝑇 = 0.1, 0.5 and 0.9, 𝜔 𝑡𝑜𝑡 2 = 10, α = 1 and 𝜎 𝜃 2 = 1.

Variance of the selection gradient

The variance of 𝛽 𝑡𝑜𝑡 is proportional to the variance of the distance of the mean phenotype to the optimum ) autocorrelation of the selection gradient is described in Chevin and Haller (2014, appendix for the autoregressive optimum). This approach cannot a priori be used in the context of the present model because the inflation factor is 1 for 𝜏 = 0 and is not defined for 0 < 𝜏 < 1. However, we found that the naïve approach to use equation ( 12) in the main text directly in the derivations from Chevin and Haller (2014) matches closely the exact discrete time approach that we develop in the next section. As we could not demonstrate this observation formally, we chose not to present this result in our main text. We nevertheless provide this result here as a heuristic:

Var[𝛽 𝑡𝑜𝑡 ] = (𝑆 𝑡𝑜𝑡 ) 2 Var[(𝑧̅ -𝜃 𝑡𝑜𝑡 )] = (𝑆 𝑡𝑜𝑡 ) 2 E[(𝑧̅ -𝜃 𝑡𝑜𝑡 ) 2 ]. ( A4
𝜌 𝛽𝑡𝑜𝑡,𝜏 = 𝜑 𝜌 𝜏𝑇 -𝑆 𝑡𝑜𝑡 𝐺 0 𝑇 𝜃 𝑇 𝑒 -𝑆 𝑡𝑜𝑡 𝐺 0 𝜏 1-𝑆 𝑡𝑜𝑡 𝐺 0 𝑇 𝜃 𝑇 , ( A9
)
where

𝑇 𝜃 𝑇
measures the characteristic timescale of the autocorrelation of the environment in units of generation time.

Discrete time approach

In contrast with previous studies (Lande andShannon 1996, Chevin andHaller 2014), we calculate the autocorrelation of the selection gradient in discrete time. Similarly to the calculation for the variance of the selection gradient, we use equation (A5) assuming that many generations occurred so that the first term can be neglected. We also use equation ( 12 (A11)

The expectation for second cross product is

𝐸[𝜃 𝑡𝑜𝑡,𝑔-𝜏 𝑧̅ 𝑔 ] = 𝑆 𝑡𝑜𝑡 𝐺 0 ∑ (1 -𝑆 𝑡𝑜𝑡 𝐺 0 ) 𝑖-1 𝐸 (𝜃 𝑡𝑜𝑡,𝑔-𝑖 𝜃 𝑡𝑜𝑡,𝑔-𝜏 ) 𝑔 𝑖=1
.

We need to distinguish 3 cases, i < τ, i > τ and i = τ : .

𝐸[𝜃 𝑡𝑜𝑡,𝑔-𝜏 𝑧̅ 𝑔 ] = 𝜑Var[𝜃 𝑡𝑜𝑡 ]𝑆 𝑡𝑜𝑡 𝐺 0 (∑ (1 -𝑆 𝑡𝑜𝑡 𝐺 0 ) 𝑖-1 𝜌 (𝜏-𝑖)𝑇 𝜏-1 𝑖=1 + ∑ (1 - 𝑔 𝑖=𝜏+1 𝑆 𝑡𝑜𝑡 𝐺 0 ) 𝑖-1 𝜌 (𝑖-𝜏)𝑇 + 𝜑 -1 (1 -𝑆 𝑡𝑜𝑡 𝐺 0 ) 𝜏-1 ). ( A12
To calculate the double sum, it is convenient to perform a change of variables

∑ ∑ (1 -𝑆 𝑡𝑜𝑡 𝐺 0 ) 𝑖+𝑘-2 𝐸[𝜃 𝑡𝑜𝑡,𝑔-𝑖 𝜃 𝑡𝑜𝑡,𝑔-𝜏-𝑘 ] 𝑔-𝜏 𝑘=1 𝑔 𝑖=1 = ∑ ∑ (1 - 𝑔 𝑗=𝜏+1 𝑔 𝑖=1 𝑆 𝑡𝑜𝑡 𝐺 0 ) 𝑖+𝑗-𝜏-2 𝐸[𝜃 𝑡𝑜𝑡,𝑔-𝑖 𝜃 𝑡𝑜𝑡,𝑔-𝑗 ].
We then need to distinguish 4 cases, . Equation (A15) is very similar to equation 5b in Chevin and Haller (2014). Comparing equation (A14) and (A15) shows that the inflation factor, i.e. the autocorrelation of selection within generation, complicates the analysis by introducing new terms at the interface between the autocorrelation of the environment and the dynamics of selection within and across generations.

A Taylor expansion performed on the reciprocals of equation (A15) and 5b in Chevin and Haller (2014) shows that both equations are identical to at least the second order. This result can be confirmed by representing both functions for different autocorrelation of the environment. Both functions match exactly on most of the autocorrelation range with minor deviation at very weak autocorrelation (not shown).

Changes in (co)variances under selection

Here we derive the changes in the variances and covariances of phenotypic components after each selection episode. Following classical quantitative genetics, we assume that the phenotype 𝑧 of

  ) highlights that with multiple selection episodes, the autocorrelation of the environment 𝜌 and the rate of response to selection 𝑆 𝑡𝑜𝑡 𝐺 0 do not have in general an additive influence on the autocorrelation of selection gradients across generations, because the factor φ multiplying the former depends on the latter, and reciprocally for γ. This occurs because with multiple episodes of selection, the relevant autocorrelation is that of the effective optimum, which integrates the effect of selection across and within generations. The autocorrelation of the selection gradient increases most with the number of selection episodes when the autocorrelation of the environment is low (fig. 4A), consistent with what we found for the inflation factor (fig.

  affecting populations growth, since both depend on the expected squared mismatch of the mean phenotype from the optimum(Chevin and Haller 2014), 𝑦 = 𝑧̅ -𝜃 𝑡𝑜𝑡 . The expected maladaptation load Lm, defined as the decrease in log mean fitness induced by a mismatch between the population mean and the optimum phenotype, is 𝐸[𝐿 𝑚 ] ] in a stationary stochastic environment, and the variance of the selection gradient is 𝑉𝑎𝑟[𝛽 𝑡𝑜𝑡 ] = 𝑆 𝑡𝑜𝑡 2 𝑉𝑎𝑟[𝑦] (eq. A5), leading to the simple relationship E[𝐿 𝑚 ] = Var[𝛽 𝑡𝑜𝑡 ] (2𝑆 𝑡𝑜𝑡 ) ⁄ . The maladaptation load in a fluctuating environment thus increases with the variance of the selection gradient and decreases with the total strength of selection. In particular, in an autoregressive environment where several selection episodes can occur within generations, we find that E[𝐿 𝑚 ] = Var[𝜃 𝑡𝑜𝑡 ]𝑆 𝑡𝑜𝑡 1-𝜌 𝑇 (1+𝑆 𝑡𝑜𝑡 𝐺 0 (𝐼 𝑓 -1)) (2-𝑆 𝑡𝑜𝑡 𝐺 0 )(1-𝜌 𝑇 (1-𝑆 𝑡𝑜𝑡 𝐺 0 )) (21) The first term in equation (21) shows that the expected component of the maladaptation load affected by responses to selection increases with the variance of the effective optimum (weighted by the strength of selection). The maladaptation load changes non-monotonically as a function of the autocorrelation of the environment (fig. 6A), as does the variance of the total selection gradient (fig 3A). Total maladaptation is strongest at intermediate autocorrelation of the environment (fig.

  . Recent refinements (e.g. McGlothlin 2010) of theory on the measurement of selection gradient (Lande and Arnold 1983, Arnold and Wade 1984a, Wade and Kalisz 1989) investigated how estimates of selection at different episodes within a generation can be integrated to infer overall selection, and the resulting evolutionary response. The present study extends this theory by investigating how fluctuating selection resulting from changes in the environment within and across generations impacts evolutionary trajectories, and the potential for population persistence.

  total per-generation individual fitness function. This allowed us to elude the difficulties arising from the dynamics of variances and covariances of breeding and phenotypic values at each successive episode of selection (see Supp. Mat. for the calculations per-episode). Indeed, the quantitative genetics framework assumes that genetic and environmental values are independent prior to selection. It is no longer the case after the first episode of selection within a generation (prior to reproduction), which generates a covariance between the above components. Tracking changes in the mean phenotype and breeding values across selection episodes therefore requires taking into account the change in covariance between breeding value and environmental values (See. Supp. Mat. for a demonstration). The buildup of this covariance complicates the investigation at the scale of individual selection episodes, from both a theoretical and empirical point of view[START_REF] Wade | The additive partitioning of selection gradients[END_REF][START_REF] Cam | Individual covariation in life-history traits: seeing the trees despite the forest[END_REF][START_REF] Mcglothlin | Combining selective episodes to estimate lifetime nonlinear selection[END_REF], because the selection gradients are no longer sufficient to describe accurately the change in the mean breeding value.

FIGURE CAPTIONS Figure 1 :

 CAPTIONS1 FIGURE CAPTIONS

Figure 2 :

 2 Figure 2: Autocorrelation of the effective optimum and inflation factor. Panel A: Autocorrelation of the effective optimum as a function of the number of generations for 10 selection episodes per generation. Panel B: Inflation factor as a function of the number of selection episodes. Panel C: Inflation factor as a function of the duration of the selection window relative to the generation time α. Dotted lines (panel A) correspond to the autocorrelation of the environment (θ) alone over g generations: 𝜌𝑔𝑇 .Crosses show results from simulations with a strength of stabilizing selection that varies across selection episodes, while lines show the expectation assuming constant selection strength with same mean (full lines: exact prediction using eq. 12, dashed lines: prediction for continuous selection along each generation infinite number of selection episodes, eq. 15). Numerical simulations are performed using eq. (10), with 𝜃 𝑡𝑜𝑡 and 𝜔 𝑡𝑜𝑡 2 calculated as described in eq. 3 from phenotypic optima 𝜃 𝑖 and corresponding strengthes of selection 𝜔 𝑖 2

Figure 3 :

 3 Figure 3: Variance of the selection gradient as a function of the number of selection episode per generation (A) and of the autocorrelation of the environment (B). Panel A: environmental autocorrelation over a generation time is 𝜌 𝑇 = 0.5 (dark gray) and 𝜌 𝑇 = 0.2 (light gray), for two total strengths of selection per generation (𝜔 ̅ 𝑡𝑜𝑡 2 = 25 [lower two sets of simulations/expectations] and 10 [upper sets]). Panel B represents the variance of the total selection gradient as a function of the autocorrelation of the environment, for the same two total strengths of selection, and n = 10. Crosses represent results from numerical simulations where the strength of stabilizing selection varies across selection episodes, whereas circles represent constant strength of stabilizing selection, as described in the caption of figure 2.

Figure 4 :

 4 Figure 4: Autocorrelation of 𝛽𝑡𝑜𝑡 over a single generation as a function of the number of selection episodes (A) and of the total strength of selection per generation 𝜔 𝑡𝑜𝑡 2 (B). Crosses represent results from numerical simulations with variable strength of stabilizing selection, whereas circles represent those with constant stabilizing selection, as described in the caption of figure 2. Dashed lines: expectations using equation 18. Continuous lines: expectation for a

Figure 5 :

 5 Figure 5: Autocorrelation function of 𝛽𝑡𝑜𝑡 across generations. The embedded graph shows a detailed view of the first five generations. Crosses represent results from numerical simulations with variable strength of stabilizing selection whereas circles represent those with constant stabilizing selection, as described in the caption of figure 2. Dashed lines: expectations using equation (18). Continuous lines: expectation for a single selection episode (eq. 17). Note that

Figure 6 :

 6 Figure 6: Maladaptation load as a function of the autocorrelation of the environment 𝜌 𝑇 (A) and as a function of the number of selection episodes (B). Crosses represent results from numerical simulations with variable strength of stabilizing selection whereas circles represent those with constant stabilizing selection, as described in the caption of figure 2. The gray scale represents different number n of selection episodes per generation in panel A, and different autocorrelation of the environment in panel B (light to dark gray, 𝜌 𝑇 = 0.2, 0.5 and 0.9). For both panels, 𝜔 ̅ 𝑡𝑜𝑡 2 = 25, 𝜎 𝜃 2 = 1, 𝛼 = 1.

Figure A1 :

 A1 Figure A1: Variance of the effective optimum as a function of the number of selection episodes (A) and as a function of the autocorrelation of the environment (B). Crosses represent the results from simulations where the strength of selection is drawn from a probability distribution as described in fig. 2 caption. Full lines represent the expectation from eq. A3 and dashed lines represent the limit for a large number of selection episodes. Panel A: from light grey to dark 𝜌 𝜃,𝑇 = 0.2, 0.5 and 0.9 respectively. Panel B: the gray scale represents different number of selection episode per generation (see within the panel). For both panels, 𝜔 ̅ 𝑡𝑜𝑡 2 = 25, 𝜎 𝜃 2 = 1, 𝛼 = 1.

  ) in the main text for the autocorrelation of the effective optimum 𝜌 𝜃𝑡𝑜𝑡,𝜏 = 𝜌 𝜏𝑇 𝜑, reminding that 𝜑 = 1 for τ = 0. The general procedure follows Chevin and Haller (2014, appendix for the autoregressive optimum) by calculating each expectation from equation (A7). Even though we kept notations concise, each sum in the expectations below can be resolved exactly. The first expectation is straightforward E[𝜃 𝑡𝑜𝑡,𝑔 𝜃 𝑡𝑜𝑡,𝑔-𝜏 ] = 𝜑Var[𝜃 𝑡𝑜𝑡 ]𝜌 𝜏𝑇 . (A10) The expectation for the first cross product is E[𝑧̅ 𝑔-𝜏 𝜃 𝑡𝑜𝑡,𝑔 ] = 𝑆 𝑡𝑜𝑡 𝐺 0 ∑ (1 -𝑆 𝑡𝑜𝑡 𝐺 0 ) 𝑖-1 𝐸[𝜃 𝑡𝑜𝑡,𝑔 𝜃 𝑡𝑜𝑡,𝑔-𝑖-𝜏 ] 𝑔-𝜏 𝑖=1 , taking into account the discontinuity of the inflation factor 𝐼 𝑓 in 𝑖 = 𝑔 -𝜏, it results in 𝐸[𝑧̅ 𝑔-𝜏 𝜃 𝑡𝑜𝑡,𝑔 ] = 𝜑Var[𝜃 𝑡𝑜𝑡 ]𝑆 𝑡𝑜𝑡 𝐺 0 0 ) 𝑔-𝜏-1 Var[𝜃 𝑡𝑜𝑡 ].

  𝑔 𝑧̅ 𝑔-𝜏 ] = (𝑆 𝑡𝑜𝑡 𝐺 0 ) 2 ∑ ∑ (1 -𝑆 𝑡𝑜𝑡 𝐺 0 ) 𝑖+𝑘-2 𝐸[𝜃 𝑡𝑜𝑡,𝑔-𝑖 𝜃 𝑡𝑜𝑡,𝑔-𝜏-𝑘 ]

  i >  and i = j, 𝐵 = Var[𝜃 𝑡𝑜𝑡 ] ∑ (1 -𝑆 𝑡𝑜𝑡 𝐺 0 ) 2𝑖-𝜏-2 𝑔 𝑖=𝜏+1 . for i >  and j > i, 𝐶 = 𝜑Var[𝜃 𝑡𝑜𝑡 ]  and i > j, 𝐷 = 𝜑Var[𝜃 𝑡𝑜𝑡 ] D, and resulting in 𝐸[𝑧̅ 𝑔 𝑧̅ 𝑔-𝜏 ] = (𝑆 𝑡𝑜𝑡 𝐺 0 ) 2 (𝐴 + 𝐵 + 2𝐶).(A13)Adding equations (A10) to (A13), finding the limit of this sum for large g, dividing by the variance of the selection gradient (eq. 16) and rearranging leads to the autocorrelation of 𝑆 𝑡𝑜𝑡 𝐺 0 )(𝜌 𝑇 -1 + 𝜌 𝑇 (𝜑 -1)𝑆 𝑡𝑜𝑡 𝐺 0 ).Equation (A14) is written in a form similar to equation (A9), decomposing the effects of the autocorrelation of the environment and of the deviations of the mean phenotype to the optimum.An interpretation for equation (A14) is proposed in the main text.In the special case where there is a unique episode of selection per generation, i.e. 𝜑 = 1, equation (A14) is 𝜌 𝛽𝑡𝑜𝑡,𝜏 = 𝜌 𝜃 𝜏𝑇 -𝜅(1-𝑆 𝑡𝑜𝑡 𝐺 0

  any individual prior to any selection in a generation is the sum of two normally distributed components: the breeding value 𝑥 and environmental effect 𝑒. At the beginning of a generation, these phenotypic components have means 𝑥̅ and 𝑒̅ = 0, variances 𝐺 and 𝐸, and are uncorrelated, such that their covariance is 𝐶 𝑥,𝑒 = 0 and the total phenotypic variance is 𝑃 = 𝐺 + 𝐸. Changes in the mean breeding value under selection are obtained from changes in mean phenotype using the regression of breeding values on phenotypes, which at the beginning of a generation has slope following the usual notation for heritability. phenotype 𝑧 have fitness (𝑧) ∝ exp (-(𝑧-𝜃) 2 2𝜔 2 ) (where we omit indices for selection episodes and generation for simplicity). It is informative to start by investigating the first selection episode in a generation, where 𝑥 and 𝑒 are independent and 𝑒̅ = 0, 𝑧̅ = 𝑥̅ . The fitness function on breeding values is obtained by integrating over the distribution of environmental effects, 𝑊 ̃(the fitness function on environmental effects is obtained similarly by integrating over breeding values, 𝑊 𝑒 (𝑒) ∝ exp (-(𝑥̅ +𝑒-𝜃) 2 2(𝜔 2 +𝐺)). The fitness function on breeding values (respectively, on environmental effects) is thus also Gaussian, with squared width 𝜔 2 + 𝑉 𝑒 (respectively 𝜔 2 + 𝐺) and optimum 𝜃 (respectively 𝜃 -𝑥̅ ). With a Gaussian fitness peak, the distribution of phenotypes, breeding values and environmental effects after selection remain Gaussian. Denoting the per-episode strength of stabilizing selection as𝑆 𝜃) = 𝜃 + (1 -𝑆𝐺)(𝑥̅ -𝜃) = 𝑧̅ -𝑆𝐺(𝑧̅ -𝜃) 𝜃) = 𝜃 -𝑥̅ + (1 -𝑆𝐸)(𝑥̅ -𝜃) = -𝑆𝐸(𝑧̅ -𝜃) ,(A16c)where we have used the fact that 𝑒̅ = 0 and thus 𝑧̅ = 𝑥̅ prior to any selection. Note that 𝑥̅ * -𝑥̅ = ℎ 2 (𝑧̅ * -𝑧̅ ) and 𝑒̅ * -𝑒̅ = (1 -ℎ 2 )(𝑧̅ * -𝑧̅ ), as required. The variances after selection are

  caption). We found that when there are only few episodes of selection per generation, variance in total selection per generation increases the variance of the selection gradient relative to what is expected in the deterministic case. However, with many selection episodes per generation, the variance in total selection decreases, resulting in a good fit between the simulations and the deterministic expectation (fig.3, compare crosses with dashed lines and circles). The effect of random selection per generation increases when selection is strong (fig.3). Lastly, we compared predictions from equation (16) calculated from a discrete generation model to previous results using a continuous time approximation(Lande and Shannon 1996 eq. 7 and table 1, and[START_REF] Chevin | Genetic constraints on adaptation to a changing environment[END_REF] Haller 2014 eq. 5a). The continuous time approximation provides more compact equations but tends to underestimate the variance of the selection gradient when selection increases (fig.3, compare dashed and continuous lines).

Young, T. P., and C. K. Augspurger. 1991. Ecology and evolution of long-lived semelparous plants. Trends in Ecology & Evolution 6:285-289.

SUPPLEMENTARY ANALYSIS

Variance of the effective optimum

By definition, the variance in the per-generation effective optimum is

where 𝜎 𝜃 2 is the stationary variance of the per-episode optimum. Rewriting (A1) using the properties of the autocorrelation of the environment (eq. 11) leads to

The variance of the effective optimum 𝜃 𝑡𝑜𝑡 is always lower than the variance of the environment (set to 1 in figure A1), and decreases with increasing number of selection episodes (fig A1 -A).

However, Var[𝜃 𝑡𝑜𝑡 ] does not decrease without bounds, instead reaching a limit value for a large number of episodes of selection, lim By rearranging equation 10 in the main text (Charlesworth 1993), the mean phenotypic value at generation g is

After a sufficiently large number of generations, the first term in equation (A5), corresponding to the initial conditions, can be neglected. We thus obtain

By developing terms, using the autocorrelation of the effective optimum (eq. 12) and assuming that 𝑔 → ∞, we obtain equation ( 16) in the main text (see also Charlesworth 1993;eq. 19a)

.

Autocorrelation of the selection gradient

The autocovariance of the selection gradient over τ generations is

(A7)

Continuous time approach

Following Lande and Shannon (1996), the dynamics can be approximated in continuous time. The continuous time approximation eases the analysis. After long enough time so that the we can neglect the initial conditions, the mean phenotype at time t is

Equation ( A8) is the continuous time equivalent to equation (A6) once the effect of initial conditions has been erased. From equations (A7) and (A8), the method to derive the

which satisfy the relationships 𝐺 * -𝐺 = ℎ 4 (𝑃 * -𝑃) and 𝐸 * -𝐸 = (1 -ℎ 2 ) 2 (𝑃 * -𝑃). The phenotypic variance after selection must also satisfy

so the covariance of breeding values and environmental effects after selection is

Using 𝑃 = 𝐺 + 𝐸 this yields

which satisfies 𝐶 𝑥,𝑒 * = 𝐶 𝑥,𝑒 * -𝐶 𝑥,𝑒 = ℎ 2 (1 -ℎ 2 )(𝑃 * -𝑃).

These formulas apply to the first selection episode, but in further episodes the change in the mean breeding value and environmental effects also depend on their covariance. More general formulas for change at any selection episode can be derived through a multivariate approach. Let us denote the vector of breeding values and environmental effects as 𝐲 = (𝑥, 𝑒) T (where superscript T denotes a transposition), with mean 𝐲 ̅ and covariance matrix 𝐘. Substituting 𝑥 + 𝑒 for 𝑧 in the fitness function 𝑊(𝑧), and integrating over the joint distribution of 𝑥 and 𝑒, the mean fitness is

Then using standard results from evolutionary quantitative genetics [START_REF] Lande | Quantitative genetic analysis of multivariate evolution, applied to brain: body size allometry[END_REF], 1980, Lande and Arnold 1983 ) and properties of the Gaussian function, the change in the vector of mean breeding values and environmental effects under selection within a generation is

where 𝜕 𝐲 ̅ is the vector of partial derivatives with respect to each element of 𝐲 ̅, yielding


Nelson Valdez

Corine Lacour

Bertrand Lecordier

Armelle Cessou

David Honore
email: david.honore@coria.fr

Experimental study of turbulent Bluff-Body flames stability by simultaneous high speed flame imaging and Particle Image Velocimetry

Development of innovative combustion regimes improving energy efficiency and limiting pollutant emissions requires better knowledge of the limits of turbulent flame stability and of the physical phenomena controlling them. The present work is focused on the description and the quantification of intermittent phenomena occurring in non-premixed turbulent flames generated by a Bluff-Body burner. The central jet of methane and the annular co-axial air flow are separated by a cylindrical Bluff-Body with a large blockage ratio. The aim is to generate a recirculation zone in the wake of the Bluff-Body in order to induce a partial mixing of reactants with trapped combustion products close to the burner exit. In such configuration, the modes of stabilization of turbulent flames are often controlled by intermittent unsteady phenomena in the recirculation zone which are complex and still poorly understood. Their study requires the ability to make time-resolved measurements by the use of high-speed optical and laser diagnostics. For this purpose, high speed flame imaging (1 kHz) is performed for three operating conditions. For a constant central methane jet velocity (i.e. constant thermal power), the variation of annular air velocity allows us to cover different flames modes from a fully stabilized turbulent flame to the limit of flame blow-out. The typical structure of the turbulent flame consists in two adjacent parts: a torus flame stabilized in the recirculation zone followed by a jet-like flame developing from the end of the recirculation zone. Time-resolved flame imaging allows to reveal the occurrence of intermittent ejections of reactive pockets from the recirculation zone to the base of the main jet-like flame, ensuring its stabilization. Image processing of time-series is performed to evaluate the characteristic frequency of this periodic ejection and its variation with the annular air flow velocity. Time-series of instantaneous velocity fields in the flames are obtained by 5 kHz Particle Image Velocimetry (PIV). Their analyses show the original structure of the flow downstream the Bluff-body burner, consisting in a highly turbulent internal recirculation zone surrounded by the centripetal annular air flow and pierced by the methane central jet. Comparison between high speed flame imaging and PIV measurements allow to determine correlations between the periodic burning pockets ejection and some features of the instantaneous velocity fields in this region, in order to fully describe this physical phenomenon controlling the stabilization of the turbulent flame.

Introduction

The requirement of the increase of efficiency and the pollutants reduction in combustion systems is still an essential task today for industrial burner manufacturers and furnaces users. Recent innovative combustion regimes with low pollutants emissions and high energy efficiency usually operate with turbulent flames close to limits of stability. It is then essential to have a deeper understanding of physical phenomena leading to flame stabilization in order to be able to master, control and enhance it. Bluff-body configuration is well adapted for the stabilization of non-premixed flames in industrial applications, such as in gas turbines, ramjets and furnaces [START_REF] Beér | Combustion Aerodynamics[END_REF][START_REF] Tong | Experimental and numerical study on bluff-body and swirl stabilized diffusion flames[END_REF]. The internal recirculation zone generated by the wake effect of the flow downstream of the bluff-body induces then the partially mixing of the fuel with oxidant. Caetano et al. and Masri et al. in their studies have observed two vortices in the recirculation zone: the inner recirculation vortex (or central jet driven) is located on the side of the fuel jet whereas the outer recirculation vortex (or air driven) is located at the edge of the recirculation zone, radially further away from the axis, in contact with the annular flow [START_REF] Caetano | A comparative experimental study of turbulent non premixed flames stabilized by a bluff-body burner[END_REF][START_REF] Masri | Turbulent diffusion flames of hydrocarbon fuels stabilized on a bluff body[END_REF]. Trapping of high-temperature combustion products in the recirculation zone brings the necessary heat for the stabilization of the flame [START_REF] Susset | Convenient laser diagnostics for aerodynamic and chemical study of axisymetric nonpremixed Bluff-Body burner flames[END_REF][START_REF] Susset | Experimental investigation of spatio-temporal correlation between aerodynamic and flame front location in an axisymetric non premixed bluff body burner flame[END_REF][START_REF] Mokaddem | Contribution à la validation expérimentale de deux modèles de combustion turbulente[END_REF].

The understanding of these stabilization phenomena in turbulent flames requires the ability to perform time-resolved measurements by use of high-rate (HR) optical and laser diagnostics. Morales et al. have investigated experimentally lean flame blowout within a high-speed combustor to analyze the temporal extinction dynamics of turbulent premixed bluff-body flames [START_REF] Morales | Mechanisms of flame extinction and lean blowout of bluff body stabilized flames[END_REF]. The lean blowout process is captured temporally using simultaneous high-speed PIV and CH*chemiluminescence. The evolution of the flame structure, flow field, and the resulting flame strain rate are analyzed throughout extinction. This study shows that global extinction can be predicted with a Karlovitz number, only after Ka exceeds 2.4 does the flame structure undergo noticeable necking and the flame begins to extinguish. Similarly in his work, Mokaddem used the flame imaging technique with a high acquisition rate (2000 images/s), to visualize the displacement of burning pockets of the Bluff-Body burner flame [START_REF] Mokaddem | Contribution à la validation expérimentale de deux modèles de combustion turbulente[END_REF]. Their results show that these structures provide the heat and reactive species necessary for the downstream jet-like flame stabilization. The flame-vortex dynamics are found to be the main driving mechanism of flame extinction. Furthermore, Tong et al. have studied flame stabilization mechanisms of a turbulent flame stabilized by combinations of swirl flow and bluff-body [START_REF] Tong | Experimental and numerical study on bluff-body and swirl stabilized diffusion flames[END_REF]. The flow is captured using high-speed PIV, while the flame structures were visualized by high-speed CH2O PLIF and CH * chemiluminescence. The global CO emissions from the flames were collected as well. They have found that the position of the outer recirculation zone is affected by the size of the bluff-body and the swirl strength. Also the recirculation zone determines the flame structures and the global CO emission levels. With a larger bluff-body, the air driven recirculation zone is located closer to the burner exit. The stabilization phenomenon is then controlled by heat transfers occurring in the recirculation zone. However, even if heat transfers play an important role, they cannot be separated from the wake aerodynamics, especially in the case of non-premixed flames, for which fluid mechanics will manage reactants residence times in the recirculation zone and thus influence the mixing process [START_REF] Nguyen | Etude expérimentale de l'influence de la geometrie du stabilisateur sur le développement d'une flamme non premlangée[END_REF][START_REF] Esquiva-Dano | Influence of a bluff-body's shape on the stabilization regime of non-premixed flames[END_REF]. In this configuration, the modes of stabilization of turbulent flames consist usually in unsteady phenomena in the internal recirculation zone which are complex and still poorly understood. Ejection of burning pockets from recirculation zone may be periodic and could be correlated with the complex aerodynamic structure of the flow induced by the bluff-body and the periodic burning pockets ejection. In the present study aerodynamics and flame structures are correlated for a turbulent non-premixed methane-air flame stabilized on a bluff-body burner, by means of coupling HR PIV with flame imaging.

Experimental setup 1.Burner configuration

To investigate turbulent flame stabilization, a bluff-body burner consisting of two coaxial jets is set up. The experimental device is presented in Figure 1. Methane is injected in the central jet (Dj = 8 mm) while air is injected in the annular one (Da = 55 mm). The Db = 48 mm-dia. bluff-body induces a large blocking ratio of 0.83.

Stability diagram

The stability diagram is established as a first step, it provides information on the operating range of the burner as well as on the overall structures of the flames stabilized on this burner. Figure 2 shows the stability diagram of the bluff-body burner. It has been established by fixing a constant methane jet velocity and increasing the annular air velocity up to the extinction limits. Several structures of the flames have been identified in the stability diagram and their global topologies are presented in Figure 3. Flame (a), occurring for low air velocity is dominated by the central jet of fuel. Flame (b) represents a transition flame characterized by the apparition of a small recirculation zone. Flame (c) reveals a blue cone with a base consisting of a laminar blue ring lifted from the burner exit. Both are broad and surround a laminar soot diffusion core. In flame (d) for low methane jet velocity and large annular air velocity, the fuel is mixed with the recirculating air and the whole is burnt before leaving the recirculation area. The flame structure (e) shows a flame leading edge anchored to the aerodynamic internal recirculation zone induced by the bluff-body wake. A jet-like flame is stabilized downstream of this leading edge. This flame structure, especially investigated in this work, is particularly of interest because it corresponds to typical operating conditions in practical burners and it may present intermittent phenomena. Flame (f) represents a lifted flame completely detached from the burner nozzle when the jet velocity of the flame is over a critical value, around 20 m/s in the present case. Upper limit of the stability diagram corresponds to the limit of flame blow-out (g). The operating conditions chosen for detailed characterization are established for a constant central methane jet velocity (Uj=9.6 m/s) and a variation of annular air velocity from Ua= 5.0 to 7.5 m/s (Figure 2). Three "type (e)" flames are retained from a fully stabilized turbulent flame (F1) chosen as reference case to a flame close to the blow-out limit (F3) passing by another stable flame (F2). The corresponding operating conditions are presented in Table 1. The air flow Reynolds number Rea is determined from the diameter of the bluff-body. A vertical laser sheet is formed by combining a spherical lens (f = 1000mm) and a cylindrical lens (f = 40mm). The laser sheet is carefully adjusted to cross the median plane of the burner and its height is around 130 mm from the burner exit. The air and methane flows are independently seeded with solid zirconium oxide particles (5 µm mean diameter) that survives the reaction zone (melting point of 2988 K). The reactants flowrates are controlled by mass flowmeters before flowing into the two fluidized beds seeding machines, allowing a careful and balanced adjustment of the particle density in each flow.

Particle images are collected with a Phantom V25 CMOS camera (1280 x 800 pixels) equipped with a Nikkor 50mm (f/1.2) and an interferential filter (centered on 527 nm). The measurement analysis area corresponds to a surface of 120 mm × 150 mm and the spatial resolution is 93 µm/pix. The acquisition rate is fixed to 5 kHz. Velocity fields are obtained by a laboratory property PIV cross-correlation algorithm [START_REF] Lecordier | Advanced PIV algorithms with image distortion validation and comparison using synthetic images of turbulent flow in 'Particle Image Velocimetry: Recent Improvements[END_REF][START_REF] Lecordier | Estimation of the accuracy of PIV treatments for turbulent flow studies by direct numerical simulation of multiphase flow[END_REF]. A multi-pass subpixel shift correlation algorithm has been used. Initial pass is computed with an interrogation window size of 16 × 32 pixels the final pass is fixed to 8 × 32 pixels (0.74 × 2.96mm²) with an overlap of 50% × 25%. Intermediate filtering is based on a minimum value of Signal to Noise Ratio (SNR) and a median filter that rejects most of non-valid vectors and resulting in instantaneous fields with more than 95% validated vectors.

The flame emission is collected through a Phantom V9 high-speed camera equipped with a double stage HiCATT 25 intensifier and a UV Cerco 105 mm lens (f/2.8). An acquisition rate at 1 kHz with a gate of 140 µs are chosen as the best compromise between signal level and time integration. The spatial resolution is 113.6 µm/pix. The high speed camera and intensifier are synchronized with the PIV laser and camera by means of a BNC-555 delay generator. Simultaneous time series of 1.5 seconds duration are recorded.

Aerodynamic characterization

The mean velocity field obtained for F1 flame is presented in Figure 4. The average of 7500 instantaneous velocity fields ensures the convergence of the statistics. Three regions can be distinguished in the flow. A first zone corresponds to the recirculation zone from the nozzle of the burner up to z = 40mm. In a second zone, between 40 mm and 60 mm, the air flow streamlines converge towards the centerline. The last zone located from 60 mm and downstream corresponds to the development of the whole flow as a jet. The recirculation zone is of particular importance in such a configuration. It is established between the central methane jet and the annular air jet and is defined by two toroidal vortices. The first large centripetal vortex is piloted by air jet and the second, smaller and centrifugal is piloted by the methane jet. These two vortices can be found in various non-premixed bluff-body configurations whatever the blocking ratio, the momentum ratio, in confined or free flow [START_REF] Susset | Convenient laser diagnostics for aerodynamic and chemical study of axisymetric nonpremixed Bluff-Body burner flames[END_REF][START_REF] Morales | Mechanisms of flame extinction and lean blowout of bluff body stabilized flames[END_REF]. The stagnation point gives the downstream boundary of the recirculation zone (respectively at r = -8 mm z = 40 mm and r = 9 mm z = 41 mm on left and right sides), which is close to the bluff-body diameter. Short-gated high-rate flame imaging is then used to get the turbulent flame structures and follow their evolution with time. Figure 6 shows a series of eighteen flames images, each separated with 2 ms. In the recirculation zone, the leading edge of the flame is seen to be lifted and anchored at the inner side of the annular air jet, as also noticed on the picture presented in Figure 3. Downstream, large regions of flame extinction are observed between the recirculation zone and the main trailing flame. On the eighteen successive flame images, a reactive pocket is detached from the flame leading edge, ejected from the recirculation zone and then transported by the flow, till the downstream zone of the depleted jet flame which is re-ignited by the flame pocket. This seems to be an intermittent phenomenon more or less observable on the time series depending on the operating conditions.

Frequency determination of reactive pockets ejection

The presence of a characteristic frequency of the reactive kernels ejection is evaluated by postprocessing of the flame emission images. Firstly, a standard deviation (std) image is obtained by subtracting the mean image on each instantaneous image of the time series. Then, the flame intensity is averaged in a specific area of 8 × 3.4 mm 2 (see Figure 6, yellow area) for each instantaneous image. This area is centered on point r=0 mm, z= 50 mm on top of the recirculation zone. A temporal evolution of F2 flame emission intensity as a function of time is shown in Figure 7. The evolution of the grey level extracted from the analysis window clearly shows the highly fluctuating nature of the turbulent flame. The characteristic frequency of a potential intermittent phenomenon is determined by means of a spectral analysis of the filtered signal. One example of FFT profile obtained for the F2 flame is presented in Figure 8. One clearly sees the presence of a peak which allows to quantify the frequency characterizing the periodic kernel ejection. This procedure is applied on three time series for each studied flame and the characteristic frequencies are gathered in Similar results are obtained for other position and size of the analysis window on the top of the recirculation zone. A typical frequency of 70 Hz is determined for F1 flame and a larger value of 90 Hz is observed for F2 flame due to the increase of air flowrate. This demonstrates the intermittent feature of the ejection of burning pockets from the recirculation zone to the main flame base. The periodic event seems to be correlated to the air flow aerodynamics since its characteristic frequency increases with air velocity. In the case of flame F3, close to the extinction limits, the flame structure does no more exhibits any periodic behavior as no characteristic frequency could be extracted from any flame emission images time-series.

Time resolved analysis of the aerodynamic structures

A sequence of successive instantaneous velocity fields obtained by PIV is presented Figure 9. The aerodynamics characteristic frequency is evaluated by extracting and averaging the radial component of the instantaneous velocity in a selected window positioned between the air and methane jets at z = 50 mm, r = 10 mm. The window area is 4.4 × 4.4 mm 2 for each instantaneous field. The windows size has been varied to ensure the reliability of the results. As previously done for flame imaging, the potential intermittent phenomena are determined by the spectral analysis of the time evolution. This procedure is applied on each operating condition, and the characteristic frequencies obtained for each flame type are presented in Table 2. The results show a perfect matching between the frequencies determined by PIV and by flame imaging for both stable F1 and F2 flames. This result points out the strong correlation of the aerodynamics with the periodic flame pocket ejection. As the air velocity increases (F1 to F3), the aerodynamics frequency increases. For the F3 flame, a high frequency is determined on each PIV time series whereas no frequency is obtained on flames images. The flame local detachment, in conditions close to blow-out, does not show any periodicity and the ejection of reactive zone cannot be correlated to aerodynamics structures, for these conditions. Nevertheless, the frequencies obtained from PIV measurements represent the intermittent aerodynamics features of the recirculation zone generated by the annular turbulent flow. The Strouhal number, characterizing the intermittent event of this flow is expressed as:

𝑆𝑆𝑆𝑆 ℎ = 𝑓𝑓•𝐷𝐷 𝑏𝑏 𝑈𝑈 𝑎𝑎
based on the air velocity Ua and bluff-body diameter Db.

Flames F1 and F2 presents similar Sth number about 0.67. This value can be linked to the natural instability of the annular air jet and corresponds to the jet preferred mode, resulting in the formation of regular vortical rings. Several studies have also shown that Sth number corresponding to the preferred mode of the jet lies in the range 0.3 -0.8 [START_REF] Danlos | Experimental characterization of the initial zone of an annular jet with a very large diameter ratio[END_REF][START_REF] Crow | Orderly structure in jet turbulence[END_REF][START_REF] Birbaud | Dynamics of free jets submitted to upstream acoustic modulations[END_REF][START_REF] Gutmark | On the preferred modes and the spreading rate of jets[END_REF]. These vortices, generated at the external base of the annular flow, feed the air driven recirculation zone and interact with the flame tube anchored downstream of the flame leading edge. The strong strain imposed by the vortex interaction on the flame leads to local extinction and periodic pockets ejection. In the case of F3 conditions, flame imaging does not present any characteristic frequency. The ratio

𝑓𝑓 𝑈𝑈 𝑎𝑎

and thus the Sth number, exceeds the value obtained for both F1 and F2 cases. The F3 flame does not seem to be controlled by aerodynamics. Under these conditions, the extinction limit is reached.

Conclusion

Turbulent Bluff-body flames stabilization phenomena were studied in the present study. The configuration studied is a non-premixed coaxial flow of air and methane. The results show that stabilization of combustion is controlled by the presence of the recirculation zone which induces some periodic ejection of burning pockets allowing ignition of the main trailing flame downstream in the flow. The phenomenon of ejection of flame pockets has a characteristic frequency, which depends on the annular air velocity injection. No specific frequency can be determined for flame close to the extinction limit.

The original flow structure downstream of the Bluff-Body burner is characterized by high-rate PIV. It consists primarily of an internal recirculation area which is defined by two vortices. The first one -a large centripetal vortex -is piloted by the annular air flow and the second one -smaller and centrifugal -is piloted by the central methane jet. Instantaneous velocity periodicity has led to the burning pockets frequency determinations which perfectly match with frequency values obtained from high-rate flame imaging for stable flames. This correlation is no more observed for a flame regime close to limit of extinction.

Figure 1 .

 1 Figure 1. The Bluff-Body burner.

Figure 2 .

 2 Figure 2. Stability diagram: (a) cylindrical flame, (b) Transition flame, (c) Ring laminar flame, (d) Recirculating flame, (e) long turbulent attached flame, (f) Lifted flame, (g) Extinction. F1, F2 and F3 refer to operating conditions studied in details in the present work.

Figure 3 .

 3 Figure 3. Direct flame visualizations of the different combustion regimes.

Figure 4 . 3 .

 43 Figure 4. Time averaged axial velocity field and mean streamlines, F1 flame.

Figure 5 .

 5 Figure 5. Averaged OH* chemiluminescence image from 1000 images, F1 flame.

Figure 6 .

 6 Figure 6. A short sequence of one time series of 1 kHz images of the F1 flame (1 image over 2 is shown).

Figure 7 .

 7 Figure 7. Temporal follow-up of the average flame intensity, F2 flame.

Figure 8 .

 8 Figure 8. FFT profile obtained for the F2 flame.

Figure 9 .

 9 Figure 9. Sequence of instantaneous axial velocity fields obtained by 5 kHz PIV (1 over 10 successive fields is shown).

Table 1 .

 1 Operating conditions.

	Uj (m/s)	Rej	Ua (m/s)	Rea	Flame
			5.0	16170	F1
	9.6	4293	6.5	21024	F2
			7.5	24258	F3

1.3. Optical diagnostics

The flow and flame structures are characterized by means of the coupling of high-rate PIV and timeresolved flame emission. The HR-PIV setup consists of a double cavity Nd:Ylf laser (Dual-DARWIN-527100-M 30 mJ @ 527 nm).

Table 2 .

 2

	0 ms	2 ms	4 ms	6 ms	8 ms	10 ms	12 ms	14 ms	16 ms
	18 ms	20 ms	22 ms	24 ms	26 ms	28 ms	30 ms	32 ms	34 ms

Table 2 .

 2 Characteristic frequencies obtained from flame imaging and PIV.

		Frequency (Hz)
	Flame Type	Flame Imaging	PIV
	F1	70 ± 3	67 ± 2
	F2	90 ± 3	91 ± 2
	F3	---	130 ± 2

Acknowledgments

The study is performed in the framework of the 'PASTEC' project with the financial support of ANR (ANR-16-CE22-0005-02).