F Schneefeld

S Schwartz

Giraud A-L

Arnal Lh

Salient 40 Hz sounds probe affective aversion and neural excitability

The human auditory system is not equally reactive to all frequencies of the audible spectrum. Emotional and behavioral reactions to loud or aversive acoustic features can vary from one individual to another, to the point that some exhibit exaggerated or even pathological responses to certain sounds. The neural mechanisms underlying these interindividual differences remain unclear. Whether distinct aversion profiles map onto neural excitability at the individual level needs to be tested. Here, we measured behavioral and EEG responses to click trains (from 10 to 250 Hz, spanning the roughness and pitch perceptual ranges) to test the hypothesis that interindividual variability in aversion to rough sounds is reflected in neural response differences between participants. Linking subjective aversion to 40 Hz steady-state EEG responses, we demonstrate that participants experiencing enhanced aversion to roughness also show stronger neural responses to this attribute. Interestingly, this pattern also correlates with inter-individual anxiety levels, suggesting that this personality trait might interact with subjective sensitivity and neural excitability to these sounds. These results support the idea that 40 Hz sounds can probe the excitability of non-canonical auditory systems involved in exogenous salience processing and aversive responses at the individual level. By linking subjective aversion to neural excitability, 40 Hz sounds provide neuromarkers relevant to a variety of pathological conditions, such as those featuring enhanced emotional sensitivity (hyperacusis, anxiety) or aberrant neural responses at 40 Hz (autism, schizophrenia).

Introduction

Beyond cultural variations in aesthetic tastes, there exists an important variability in the way individuals are sensitive and reactive to sounds, whether pleasant (e.g., baby laughs, ASMR) or unpleasant (e.g., screams, chalk on a blackboard). In extreme cases, sound hypersensitization can trigger pain [START_REF] Tyler | A Review of Hyperacusis and Future Directions: Part I. Definitions and Manifestations[END_REF], audiogenic seizures [START_REF] Italiano | Genetics of reflex seizures and epilepsies in humans and animals[END_REF] or excessive stress sometimes resulting in aggressive reactions such as in the shaken baby syndrome [START_REF] Li | Explaining individual variation in paternal brain responses to infant cries[END_REF]. Given sound potential pathogenicity, it is critical to understand the neural mechanisms underlying such variability between individuals. Yet, the sonic features that drive hypersensitization and how they affect neural circuitry remain unclear.

Recent works have shown that alarm auditory signals (e.g., screams, sirens) exploit a restricted 30-80 Hz range of temporal modulations known as roughness to communicate danger, trigger negative affective responses and elicit rapid reactions [START_REF] Arnal | Human Screams Occupy a Privileged Niche in the Communication Soundscape[END_REF]. Unlike conventional noise, roughness can trigger aversive responses at relatively mild sound intensity [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF]. Therefore, it is particularly well-suited to investigate the neural underpinnings of auditory aversion in humans. Mapping behavioral and brain responses in the roughness range and beyond, intracranial EEG studies recently revealed that the nonlinear pattern of aversion across frequencies strikingly resembles the spectral profile of a well-known auditory neural signal, the auditory steady-state response (ASSR) [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF]. This suggests that aversion to rough sounds might relate to the capacity of repetitive transients to induce sustained responses in the brain. However, subjective aversion profiles, though highly reliable across experiments when averaging data across participants, show substantial variability across participants. If the link between aversion and brain responses holds at the individual level, then we should expect similar variability in the ASSR signals of our participants. Here, we primarily sought to explore, using non-invasive surface EEG, whether aversive sound rating patterns relate to neural excitability at the individual level. More specifically, we tested the hypothesis that increased aversion to roughness could be underpinned by enhanced auditory steady-state responsiveness to sounds.

Electrophysiological brain signals are known to maximally entrain in a sustained manner to repetitive acoustic transients at 40 Hz. Since the study by Galambos and colleagues [START_REF] Galambos | A 40-Hz auditory potential recorded from the human scalp[END_REF], a wealth of experimental work has investigated this phenomenon. Although often used as an index of early auditory process integrity [START_REF] Korczak | Auditory Steady-State Responses[END_REF], it has been hypothesized that these signals more generally index long-range thalamocortical synchronization and information transfer in the brain [START_REF] Krishnan | Steady state and induced auditory gamma deficits in schizophrenia[END_REF][START_REF] Picton | Human Auditory Steady-State Responses: The Effects of Recording Technique and State of Arousal[END_REF]. The 40-Hz ASSR is often considered a promising biomarker of various neuropsychiatric [START_REF] Spencer | Gamma-Band Auditory Steady-State Responses Are Impaired in First Episode Psychosis[END_REF][START_REF] Thun | The 40-Hz auditory steady-state response in patients with schizophrenia a meta-Analysis[END_REF], neurodevelopmental [START_REF] Lehongre | Altered Low-Gamma Sampling in Auditory Cortex Accounts for the Three Main Facets of Dyslexia[END_REF][START_REF] Wilson | Children and Adolescents with Autism Exhibit Reduced MEG Steady-State Gamma Responses[END_REF] and neurodegenerative disorders [START_REF] Van Deursen | 40-Hz steady state response in Alzheimer's disease and mild cognitive impairment[END_REF]. However, the link between these disorders and auditory processing is unclear. While some disorders are comorbid with hearing problems [START_REF] Griffiths | How Can Hearing Loss Cause Dementia?[END_REF][START_REF] Livingston | Dementia prevention, intervention, and care: 2020 report of the[END_REF], the prevalence of ASSR impairments in these diseases might point to dysfunctions beyond or outside classical auditory pathways. Indeed, whereas previous EEG studies a priori locate the source of ASSR generators in the auditory cortex, recent findings using alternative reconstruction algorithms or spatially resolved intracranial recordings suggest that rough sounds additionally synchronize non-canonical auditory systems involved in salience and aversion processing [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF][START_REF] Farahani | Spatiotemporal reconstruction of auditory steady-state responses to acoustic amplitude modulations: Potential sources beyond the auditory pathway[END_REF][START_REF] Tada | Global and Parallel Cortical Processing Based on Auditory Gamma Oscillatory Responses in Humans[END_REF]. Therefore, impaired ASSR responses -as observed in neuropsychiatric disorders-might actually reflect dysfunctional sensory and emotional processing of exogenous stimulation in these networks.

To summarize, previous findings demonstrate (i) that 40 Hz sounds are perceived as aversive, (ii) that they recruit canonical and non-canonical auditory networks involved in aversion processing, and (iii) that the neural processing of these sounds might be altered in a variety of brain disorders featuring abnormal emotional responses. These converging observations suggest that the neural excitability of aversion networks -as revealed by the strength of ASSR-might influence emotional reactivity to these sounds, as well as indicate abnormal excitability in neuropsychiatric disorders. Here, we aimed to clarify the link between emotional reactivity and neural excitability to rough sound at the individual level in the normal population. In line with previous observations that ASSR are compromised in neuropsychiatric and neurodegenerative diseases classically featuring enhanced anxiety, we further tested how these responses interact with individual anxiety level.

Methods

Stimuli

All stimuli were digitally generated using MATLAB with a sampling rate of 96 kHz and presented in a pseudo-randomised order using Psychtoolbox (Version 3.0.12). Click trains (1 s duration, with 100 ms sine ramping onset and offset; click rise/fall time of 0 ms, plateau time of 1 ms, presented at ~60 dB SPL at frequencies varying between 10 and 250 Hz).

Paradigm

The first behavioral experiment (Exp. 1, pilot, N=12, 6 females) measured subjective aversion as a function of the frequency of click trains [START_REF] Spencer | Gamma-Band Auditory Steady-State Responses Are Impaired in First Episode Psychosis[END_REF][START_REF] Kim | Cortically projecting basal forebrain parvalbumin neurons regulate cortical gamma band oscillations[END_REF]40,60,80,100,120,140,170,200 and 250 Hz). After the presentation of each stimulus, participants reported their aversion to click trains of varying frequencies on a 5-level Likert-like subjective scale ranging from tolerable (1) to unbearable [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF]. The second experiment (Exp. 2, N=38, 25 females) replicated Exp. 1 while recording participants' EEG signals. In Exp. 2, stimuli (click trains at 10, 20, 30, 40, 60, 90, 130, 180 and 250 Hz) were presented in pairs, separated by a period varying from 800 to 1200 ms. Participants were required to rate the unpleasantness of the second stimulus of the pair only, to avoid memory effects, on an analog slider scale ranging from 1 (tolerable) to 10 (unbearable). Pairs were made of either the same (repeated, 80% of all pairs) or different (oddballs, 20% of all pairs) stimuli. Only pairs containing repeated stimuli were analyzed in the current study. All (oddball and repeated) pairs were presented in pseudorandomized order. Two participants did not complete the task correctly and were removed from subsequent analyses.

Behavioral data analysis

Individual aversion profiles in Exp. 1 and Exp. 2 were obtained by averaging aversion scores per stimulus frequency, which were subsequently z-scored per participant to account for baseline individual aversion biases (inter-individual variability of overall sound aversion across all frequencies). Qualitative inspection of individual profiles suggests the existence of two distinct aversion profiles: on the one hand, roughness responders exhibiting a non-linear response with heightened aversion to low frequencies (low, 30-60 Hz), and on the other hand, pitch responders, in whom aversion follows a linear profile with a maximal response for the highest frequencies (high, >130 Hz; Fig. 1). To quantify this dual patterning of aversion profiles across participants -namely, to test if those participants sensitive to high frequencies were less sensitive to low frequencies and conversely-, we computed the correlation of aversion between each tested frequency across participants. Computing the sign of the difference between aversion scores in low minus high [130-250 Hz] frequency ranges, we then categorized participants as roughness responders (low -high>0) or pitch responders (low -high<0).

To check the replicability of aversion profiles between Exp. 1 and Exp. 2 and to be able to subsequently pool the results across experiments, we then applied the following modelling procedure. As we selected slightly different stimulus frequencies (in a similar frequency range) for Exp.1 and Exp. 2, we interpolated the data at each frequency in the selected frequency range. Stemming from the hypothesis that non-linear profiles result from the differential weighting of non-linear roughness and pitch profiles (Figure 1), we then sought to fit individual behavioral data aversion scores. In each participant, we fit a summation between a gaussian distribution and a linear model (Supplementary Figure 1 and the equation below):

𝑃 ! = 𝐺 √2𝜋𝜎 " 𝑒 # (%#&) ! ( ! + 𝑚𝑓 + 𝑏
where 𝑃 ! is the predicted rating of a given trial and 𝑓 is the stimulus frequency of that trial. The model contains 5 parameters: 3 for the gaussian shape, a mean 𝜇, a standard deviation 𝜎 and a gain 𝐺 to normalize the height of the peak to fit in the scale of the subjective ratings; and 2 parameters for the linear portion, slope 𝑚 and intercept 𝑏. The parameters were fit by minimizing the squared error between the model predictions and the true sound rating taking the stimulation frequency as input. Minimization was performed using the fminsearch function in MATLAB. The model allowed us to interpolate our results to untested frequencies and better compare across experiments using different stimulus frequencies.

EEG analyses

EEG recording was performed using the DSI 24 EEG headset (Wearable Sensing, San Diego, USA). This headset is equipped with 21 dry electrodes, placed according to the 10/20 international standard and permits wireless data transmission at a 300 Hz sampling rate. Two additional reference electrodes were placed at the ears. The data was common-average re-referenced and chunked in 1500 ms epochs (500 ms before to 1000 ms after the click train onset). For each trial, we then converted the data into z-scores using the pre-stimulus baseline activity (from 200 ms before to stimulus onset) for each electrode. All EEG data analyses were performed in MATLAB (MathWorks) using the Fieldtrip (http://fieldtrip.fcdonders.nl) package and in-house custom code. Trials were visually inspected, and those with obvious artifacts were removed. An independent component analysis as implemented in FieldTrip was used to correct for eye blinks and eye movements. Adding to the two participants who did not perform the task correctly, three participants were excluded from the EEG analysis due to poor EEG signal quality, leaving thirty-three participants to be included in EEG analyses. Auditory ERPs and ASSRs were calculated for each individual, electrode and stimulus frequency. For ASSRs, we measured the power of neural responses filtered at the stimulus frequency. Amplitude time series were computed by band-pass filtering the EEG signal at the stimulus frequency F ± 0.5 Hz. The envelope of each narrow-band signal was obtained by taking the absolute value of the analytic signal resulting from a Hilbert transform. ASSR values were baseline-corrected by subtracting the power of the EEG signal at stimulus frequency during a 200-ms period preceding sound onset to capture endogenous activity. We then transformed these time-courses into z-scores per trial and electrode and participant. Given that the absolute magnitude of auditory responses (ERPs) substantially varied across participants (whether for SNR or physiological grounds), we quantified ASSR relative to the canonical auditory ERPs to high [130-250 Hz] sounds. We selected these frequencies because ASSR power is negligeable and minimally interacts with ERP components at these frequencies. As our goal was to explore and relate inter-individual brain responses to a variety of measurements (subjective aversion, anxiety) and given the relative SNR disparity across participants, we standardized ASSR responses in each participant as follows. We subtracted early evoked response to high frequencies (value of the N1 in the 3 central sensors) to the ASSRs at each frequency and electrode. Reasoning that roughness responders should exhibit larger responses to low than high sounds (and conversely for pitch responders), this procedure allowed standardizing these measurements across participants, and provides a relative metric to assess low versus high excitability at the individual level.

Finally, to avoid a potential contamination of entrained ASSR by ERPs to sound onsets (e.g., N1) and by motor preparation responses, we focused our analysis on a time window-defined by data visual inspection of the data-from 200 to 600 ms after stimulus onset.

Correction for Multiple Comparisons

To assess the statistical difference between the experimental conditions while controlling for multiple comparisons, we performed nonparametric cluster analyses [START_REF] Maris | Nonparametric statistical testing of EEG-and MEG-data[END_REF] for each statistical test. This method is based on comparing the cluster statistic of adjacent timepoints (significant at p<0.05) to the same statistic ran on randomly permuted data. The nonparametric statistic was performed by repeating 1000 times the calculation of a permutation test where the experimental conditions are randomly intermixed within each subject. For each of these permutations, we extracted the maximum (absolute value) cluster-level statistic. Finally, we calculated the corrected p-values by comparing the cluster-level statistics of our original data with those of all permutations. All the results presented are corrected for multiple comparisons using this method unless otherwise stated.

Results

Subjective aversion: roughness versus pitch responders

To investigate the inter-individual variability or subjective aversion across frequencies, we first analyzed the data from a previous behavioral experiment (Experiment 1, see [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF]. Twelve participants listened to a series of 1-s click trains of varying click frequency (ranging from 10 to 250 Hz) and rated the sounds' aversiveness from 1 (tolerable) to 5 (unbearable). Figure 1a (left panel) shows that aversion varied in a non-linear, non-monotonic manner with increasing sound rates. Qualitatively, individual rating profiles appeared quite variable across participants, some exhibiting higher aversion to low than high frequencies and others exhibiting the inverse pattern. To quantify this preliminary observation, we tested the correlation of ratings across participants and frequencies (Figure 1a, middle left panel). We observed that while aversion patterns were very similar across participants between adjacent frequencies (along the diagonal of the correlation matrix), aversion patterns were negatively correlated between low (<80 Hz) and high frequencies (>130 Hz). This effect reveals a dual pattern of auditory aversion across participants: while a subset of participants was most sensitive to high frequencies and relatively less to low (pitch responders, N=6), the other half (roughness responders, N=6) was more sensitive to lower, temporally salient frequencies but not to high frequencies. To categorize participants as roughness or pitch responders (Figure 1a, middle right panel), we measured the difference between averaged ratings of low (30-60 Hz) and high (130-250 Hz) sounds. This allowed us to assign each participant to either a roughness responders (low minus high > 0) or a pitch responders (low minus high < 0) group. Figure 1a (rightmost panel) reproduces Figure 1a (leftmost panel) color coded as a function of roughness (red) versus pitch (blue) responder groups.

To replicate these behavioral results and investigate their underlying neural correlates, we repeated the same experiment (Experiment 2) using slightly different stimulus rates (ranging from 10 to 250 Hz) in 33 participants while recording electroencephalographic (EEG) responses to these sounds. Focusing on behavioral data, we found that the pattern of ratings and ensuing analyses (Figure 1b) replicated the results obtained in Exp. 1, although the ratio of pitch (N=23) versus roughness (N=10) responders was slightly less balanced (Figure 1b). Finally, to pool the data from both experiments for subsequent analyses, we developed a non-linear curve fitting procedure (see Methods and supplementary data section) which reproduced the dual pattern of roughness versus pitch responsiveness to aversion (Figure 1c). To test the replicability of this non-linear pattern between experiments 1 and 2, we measured the Pearson's correlation between interpolated aversion profiles (averaged across participants within each experiment) and confirmed that these profiles were highly similar between the two experiments (Pearson's r = 0.79, p = 1e-53).

Event related potentials relate to subjective aversion

To analyze brain responses to click trains of varying rates, we first measured the magnitude of ERPs across channels and participants (Figure 2a). This analysis replicated most features of the main results previously observed in intracranial recordings [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF]. Focusing on the early, N1 component, we observed that the magnitude of this component linearly increased as a function of click trains frequency (Pearson's r=-0.91, p= 0.0007). Measuring the amplitude of ERPs in a late [200-600 ms] time-window, we observed that brain responses to low frequencies (low, 30-60 Hz) were more negative than responses to higher (high, >130 Hz) frequencies (t(32)=9.22, p=0.01) with a profile that qualitatively resembled the non-linear aversion rating profiles.

We first compared the magnitude of ERPs across all stimulus frequencies between the two groups of participants over central N1 channels (Figure 2c). This analysis did not yield significant differences between groups, suggesting similar levels of activity in responses to click trains. However, our primary hypothesis was that roughness and pitch responders might exhibit different neural responsiveness to low versus high frequencies, respectively. We therefore sought to compare the relative ERP magnitude to low versus high stimulus modulation frequencies between the two groups. To do so, we first computed the difference between low and high ERPs in each group (Figure 2d). This analysis revealed significant differences in early components (P50) as well as several significant clusters in the late [0.2-0.6 s] time-window. We then tested the Pearson's correlation between individual aversion ratings and the low minus high ERP difference across time, and found that the correlation was maximal and significant in the late time-window (Pearson's correlation r[200-600 ms] =0.48 , p= 0.004).

Auditory Steady-State Responses to low sounds reflect subjective aversion

Based on the hypothesis that temporally salient low sounds should entrain brain responses in a steadystate, sustained manner, we measured the power of steady-state responses (ASSR) at the rate of the exogenous stimulation in the late time window (restricted to [200-600 ms] to avoid potential contamination by onset and offset responses). ASSR power averaged across all electrodes showed a sustained response throughout the sound duration (Figure 3a, left panel). Focusing on brain responses to low sounds (30-60 Hz) across all participants (both groups together), we observed that ASSR power in this range was maximal in central and anterior contacts, a topographical pattern that is compatible with the recruitment of auditory and frontal regions as observed in intracranial recordings in humans [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF][START_REF] Tada | Global and Parallel Cortical Processing Based on Auditory Gamma Oscillatory Responses in Humans[END_REF] and murine models [START_REF] Kim | Cortically projecting basal forebrain parvalbumin neurons regulate cortical gamma band oscillations[END_REF].

We then aimed to assess the difference in ASSRs for roughness and pitch responders in centralanterior sensors (Figure 3b). Despite the qualitative difference between roughness and pitch responders, this contrast did not reach statistical significance. We then tested the correlation in time between low ASSRs and individual subjective ratings in the same set of electrodes and identified several correlation clusters in the late window. Averaged over the late time-window, the correlation proved to be maximal and significant at 40 Hz (Figure 3c).

Anxiety increases aversion and ASSRs to low frequency sounds

We then tested the hypothesis that anxiety might interact with the excitability of aversion networks. We tested the relationship between individual anxiety levels (as assessed using the STAI questionnaire) and (i) subjective aversion and (ii) brain responsiveness to low sounds. Focusing on subjective ratings from Exp. 1, we first determined that aversion ratings significantly correlated with Trait-anxiety scores (Figure 4a). In other words, while the most anxious participants rated rough, low (30-60 Hz) sounds as more unpleasant, the least anxious participants rated louder, high trains (>180 Hz) more negatively. Replicating this analysis in Exp. 2, we found a similar pattern of correlation with STAI-T scores (Pearson correlation between Exp. 1 and 2 correlation profiles: r = 0.929; p = 0.001), although non-significant. In addition, pooling interpolated data from Exp. 1 and 2 (see Fig. 1c) confirmed that across all participants and experiments, trait anxiety positively correlated with subjective aversion to low and negatively correlated with high sounds. Focusing on ASSR responses to low sounds, we then measured the relationship between interindividual anxiety (STAI-T) scores and cerebral responsiveness to aversive sounds. We observed that whereas anxiety does not interact with early auditory ERPs, this analysis revealed a significant correlation between anxiety and strength of late [200-600 ms] sustained auditory steady-state responses (ASSR) to low sounds (Fig. 4b). Focusing our analysis of ASSR in the late time-window, we found that the correlation was maximal and significant only at 30 and 40 Hz.

Discussion

Aversion to roughness links to steady-state excitability

Here, we investigated the relationship between individual subjective aversion and neural responsiveness. Replicating our previous work [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF], we show that aversion to repetitive sounds in the roughness range and beyond is not a linear function of stimulus frequency. When averaged across participants, aversion follows a non-linear pattern, peaking at 40 Hz in the roughness range and linearly increasing as a function of stimulus frequency above 130 Hz. However, further investigating inter-individual variability, we found that although highly reliable across experimental replications (see [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF] and current study), this pattern substantially varied across participants. While some participants' aversion profile followed a linear pattern seemingly determined by stimulus frequency (here proportional to sound intensity), other participants' aversion profile revealed a higher sensitivity to sounds in the roughness range, most prominent at 40 Hz. Exploring the variability of individual aversion profiles, we discovered that some participants were maximally sensitive to high frequency sounds (pitch responders), others to low frequencies (roughness responders) or exhibited a composite pattern between these two extremes. When studying the neurophysiological origin of the non-linear pattern of aversion across frequencies using intracranial recordings in epileptic patients, we recently suggested that unlike aversion to higher frequencies in the pitch range (>130 Hz), aversion to rough frequencies is best accounted for by the sustained entrainment of activity of brain regions involved in salience and aversion processing [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF]. In other words, beyond triggering classical onset related potentials (ERPs) in the auditory cortex, the nonlinear effect of subjective aversion in the roughness range is arguably induced by the additional steadystate neural entrainment (ASSR) of widespread non-canonical auditory networks [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF][START_REF] Farahani | Spatiotemporal reconstruction of auditory steady-state responses to acoustic amplitude modulations: Potential sources beyond the auditory pathway[END_REF][START_REF] Tada | Global and Parallel Cortical Processing Based on Auditory Gamma Oscillatory Responses in Humans[END_REF]. Considering the disparity of aversion profiles between participants, we then hypothesized that one possible origin of this inter-individual variability might be the relative excitability of these networks and related ASSR responses across participants. We predicted that individual sensitivity to roughness (versus stimulus frequency) might be determined by the relative magnitude of ASSR and ERP signals, respectively. Relating these EEG signals to individual aversion profiles across 33 participants, we found that aversion to rough (low) sounds correlated with the magnitude of ASSR signals (relative to high ERP magnitude). This supports our initial prediction, namely that sensitivity to rough sounds is linked to neural steady-state responsiveness, at the individual level. In light of the notion that ASSRs reflect the entrainment of widespread, non-canonical neural networks involved in salience processing, this finding suggests that the relative aversion to rough sounds, in particular at 40 Hz, might provide a valuable proxy of inter-individual excitability in those networks.

Anatomical and functional sources of ASSRs and related aversion

The phenomenological distinction between roughness and pitch is commonly attributed to a dual coding strategy of temporal and spectral information in the canonical auditory pathway [START_REF] Bendor | The role of inhibition in a computational model of an auditory cortical neuron during the encoding of temporal information[END_REF][START_REF] Brugge | Coding of repetitive transients by auditory cortex on Heschl's gyrus[END_REF][START_REF] Joris | Neural Processing of Amplitude-Modulated Sounds[END_REF]. However, recent results suggest that these two response modes may reflect the recruitment of distinct neural routes for sound processing. While pitch sounds are processed along the canonical auditory pathway, rough sounds may additionally target an alternative, non-canonical pathway ultimately targeting medial cortical and limbic regions 5,16,17 . Here, we hypothesize that the inter-individual variability in neural and emotional responses to 40 Hz sounds might be accounted for by increased excitability of these non-canonical networks relative to canonical auditory pathways.

The anatomical origin and trajectory of such pathways, and how they trigger negative emotional responses remain elusive. A recent stream of findings in mouse models has shown that 40 Hz sounds entrain neural activity in medial temporal and limbic -not primarily auditory-regions notably involved in arousal, emotional and attentional processing 19,24 . Interestingly, rough sounds transiently interrupt microsaccadic rhythms, suggesting that they impact oculomotor processes via non-primary auditory brainstem nuclei. That these sounds interact with non-auditory attentional processes concurs with other observations highlighting their impact on aversion (5), vigilance state and arousal [START_REF] Picton | Human Auditory Steady-State Responses: The Effects of Recording Technique and State of Arousal[END_REF][START_REF] Griskova | The amplitude and phase precision of 40 Hz auditory steady-state response depend on the level of arousal[END_REF]. This further highlights the propensity of rough sounds to affect early arousal and low-level emotional processes in a way that possibly prescinds higher-order auditory functions. Although more experimental work is needed to decipher the circuits involved and their influence on affective responses, the current results support the idea of an additional, non-canonical auditory circuitalthough possibly overlapping with classical ones-involved in the affective, non-linear processing of roughness.

Mapping between affective and neural responsiveness to rough sounds

Although brain responses to 40 Hz sounds are widely studied in neuroscience and frequently used in audiological assessment [START_REF] Picton | Estimating audiometric thresholds using auditory steady-state responses[END_REF], why and how these sounds affect arousal and induce negative emotional responses remains neglected. It has been known for decades that steady-state EEG responses to 40 Hz sounds are affected in various brain disorders, whether neuropsychiatric [START_REF] Spencer | Gamma-Band Auditory Steady-State Responses Are Impaired in First Episode Psychosis[END_REF][START_REF] Thun | The 40-Hz auditory steady-state response in patients with schizophrenia a meta-Analysis[END_REF], neurodevelopmental [START_REF] Lehongre | Altered Low-Gamma Sampling in Auditory Cortex Accounts for the Three Main Facets of Dyslexia[END_REF][START_REF] Wilson | Children and Adolescents with Autism Exhibit Reduced MEG Steady-State Gamma Responses[END_REF] or neurodegenerative [START_REF] Van Deursen | 40-Hz steady state response in Alzheimer's disease and mild cognitive impairment[END_REF][START_REF] Horvath | EEG and ERP biomarkers of Alzheimer's disease: a critical review[END_REF][START_REF] Nimmrich | Neuronal Network Oscillations in Neurodegenerative Diseases[END_REF]. Interestingly, many of these pathologies often feature (i) deficits in salience processing and arousal ( 28) and (ii) excitation/inhibition imbalance in medial temporal, limbic and salience-related networks [START_REF] Bookheimer | Patterns of Brain Activation in People at Risk for Alzheimer's Disease[END_REF][START_REF] Dickerson | Increased hippocampal activation in mild cognitive impairment compared to normal aging and AD[END_REF][START_REF] Palop | Aberrant Excitatory Neuronal Activity and Compensatory Remodeling of Inhibitory Hippocampal Circuits in Mouse Models of Alzheimer's Disease[END_REF][START_REF] Palop | Network abnormalities and interneuron dysfunction in Alzheimer disease[END_REF]. Many of these disorders further evidence atypical affective responses and are often associated with increased anxiety [START_REF] Ferretti | Anxiety and Alzheimer's disease[END_REF][START_REF] Pavlova | Lifetime prevalence of anxiety disorders in people with bipolar disorder: a systematic review and meta-analysis[END_REF][START_REF] Temmingh | Anxiety in Patients with Schizophrenia: Epidemiology and Management[END_REF][START_REF] Wigham | The interplay between sensory processing abnormalities, intolerance of uncertainty, anxiety and restricted and repetitive behaviours in autism spectrum disorder[END_REF]. However, the link between neural excitability and emotional reactivity to these sounds -whether in neuropsychiatric or neurotypical populations -remains unclear.

Anxious people often exhibit exaggerated responses to salient or unpleasant sensory events [START_REF] Tyler | A Review of Hyperacusis and Future Directions: Part I. Definitions and Manifestations[END_REF]. Here, we show that interindividual sensitivity to low sounds (relative to high ones) correlates both with neural responsiveness to these sounds and with individual anxiety trait. In other words, the most anxious participants are more averse and show stronger neural entrainment to rough, low sounds than participants whose sound aversion is mostly driven by the frequency (here confounded with intensity) of the sounds. Considering recent works suggesting that 40 Hz sounds entrain non-canonical networks involved in aversion processing [START_REF] Arnal | The rough sound of salience enhances aversion through neural synchronisation[END_REF], these results suggest that these networks are more excitable in anxious participants, who in turn perceive these sounds as more unpleasant. As these effects more tightly link neural and emotional responsiveness to anxiety trait rather than state, one could conjecture that such systems are overall more excitable in anxious participants rather than being over activated at the time of the experiment.

Altogether, these results support the prediction that emotional and electrophysiological responses to 40 Hz sounds concurrently probe responsiveness of aversion networks and participants' anxiety at the inter-individual level. The neuroanatomical and functional origin of such networks remains to be clarified e.g., in animal or pharmacological models focusing on emotional reactivity to these sounds. However, adding to previous works showing aberrant 40 Hz ASSRs in a variety of neuropsychiatric disorders, our results point to the possibility that anxiety implicates -and can be probed by-the overresponsiveness of brain networks involved in aversion to rough sounds.

Figure 3 :

 3 Figure 3: Time-frequency analysis of low steady-state responses (ASSR). a. Left: time course of ASSR responses at each stimulation frequency. Middle: topography of Low (30-60 Hz) ASSR power averaged over the late [.2-.6 s] time-window. Right: ASSR power per stimulation frequency averaged across the late time-window. Grey dots represent individual data points. b. Time-course of low ASSR power for low (red) and high (blue) responders. Upper right inset: averaged ASSR power difference between low and high responders in the late-time time window. c. Time-course of Pearson's correlation between low ASSR power and aversion ratings. Black lines represent significant correlation values at p<.05 level (cluster corrected). Upper right inset: Pearson's correlation r-values per stimulation frequency in the late-time time window. Star indicates significance at p<.05.

Figure 4 :FiguresFigure 1 :

 41 Figure 4: Relationship between auditory aversion and anxiety at the inter-individual level a. Pearson's correlation between aversion ratings and individual anxiety levels assessed via STAI questionnaires for STAI-Trait (red) and STAI-State (blue). Shaded areas represent SEM of the correlation. Left: pilot experiment. Middle: EEG experiment. Right: same as in left and middle panels grouping interpolated data from the pilot and the EEG experiment (see supplementary information). Time-course of Pearson's correlation between Low ASSR power and anxiety STAIT-T scores. Black lines represent significant correlation values at p<.05 level (cluster corrected). Right: Pearson's correlation r-values per stimulation frequency in the late-time time window. Star indicates significance at p<.05. Right: topography of correlation r-values between Low ASSR power and anxiety STAIT-T scores.

Figure 2 :

 2 Figure 2: ERP responses across frequencies. a. ERP responses across time and averaged across all EEG channels as a function of stimulus frequency. b. Top left panel: Topography of the N1 response. Top right panel: N1 magnitude as a function of stimulus frequency. Bottom left panel: Topography of the late [0.2 -0.6 s] responses. Bottom right panel: late response magnitude as a function of stimulus frequency. Selected channels (exhibiting maximal magnitude) are highlighted in white c. ERP responses across frequencies do not differ between Low (red) and high (blue) responders. d. Top panel: difference between low [30-60Hz] and high [130-250 Hz] frequency ERPs for low (red) and high (blue) responders. Stars indicate significant (p<.05 corrected) statistical difference. Bottom panel: Correlation between aversion ratings and low minus high ERPs. black dots indicate (cluster corrected) statistical significance at p<.05 level.

Figure 3 :

 3 Figure 3: Time-frequency analysis of low steady-state responses (ASSR). a. Left: time course of ASSR responses at each stimulation frequency. Middle: topography of Low (30-60 Hz) ASSR power averaged over the late [.2-.6 s] time-window. Right: ASSR power per stimulation frequency averaged across the late time-window. Grey dots represent individual data points. b. Time-course of low ASSR power for low (red) and high (blue) responders. Upper right inset: averaged ASSR power difference between low and high responders in the late-time time window. c. Time-course of Pearson's correlation between low ASSR power and aversion ratings. Black lines represent significant correlation values at p<.05 level (cluster corrected). Upper right inset: Pearson's correlation r-values per stimulation frequency in the late-time time window. Star indicates significance at p<.05.

Figure 4 :

 4 Figure 4: Relationship between auditory aversion and anxiety at the inter-individual level a. Pearson's correlation between aversion ratings and individual anxiety levels assessed via STAI questionnaires for STAI-Trait (red) and STAI-State (blue). Shaded areas represent SEM of the correlation. Left: pilot experiment. Middle: EEG experiment. Right: same as in left and middle panels grouping interpolated data from the pilot and the EEG experiment (see supplementary information). Time-course of Pearson's correlation between Low ASSR power and anxiety STAIT-T scores. Black lines represent significant correlation values at p<.05 level (cluster corrected). Right: Pearson's correlation r-values per stimulation frequency in the late-time time window. Star indicates significance at p<.05. Right: topography of correlation r-values between Low ASSR power and anxiety STAIT-T scores.

  To conclude, the current work shows that aversion to 40 Hz sounds and ensuing brain responses constitute viable biomarkers to probe emotional responsiveness and neural excitability in non-canonical auditory networks. Such neuromarkers open new perspectives to assess the effect of therapeutic interventions on neural and emotional responsiveness at the individual level, in a variety of brain disorders featuring excessive emotional reactivity or neural responsiveness to temporally salient, rough stimulation. Behavioral aversion ratings. a. Left panel: Pilot experiment: behavioral aversion ratings as a function of stimulus frequency. White circles and error bars represent mean across participants +/-SEM. Grey lines correspond to individual aversion curves. Middle left: correlation of aversion scores across frequencies and individuals. Positive correlation values (yellow) show that interindividual variability is preserved across neighboring frequencies. Negative correlation values (blue) show that participants who rate low frequencies (30-60 Hz, low) as more aversive rate high frequencies (130-250 Hz, high) as less aversive, and reciprocally, evidencing a distinction between low and high responders. Middle right: determination of low minus high responder groups: participants with positive low minus high aversion values are categorized as low responders and reciprocally for high responders. Right panel: individual aversion rating curves categorized as low (red) and high (blue) responders. b. Same as in (a.) for the EEG experiment. c. Same as (a.) and (b.) grouping interpolated data from the pilot and the EEG experiment (see supplementary information). ERP responses across frequencies. a. ERP responses across time and averaged across all EEG channels as a function of stimulus frequency. b. Top left panel: Topography of the N1 response. Top right panel: N1 magnitude as a function of stimulus frequency. Bottom left panel: Topography of the late [0.2 -0.6 s] responses. Bottom right panel: late response magnitude as a function of stimulus frequency. Selected channels (exhibiting maximal magnitude) are highlighted in white c. ERP responses across frequencies do not differ between Low (red) and high (blue) responders. d. Top panel: difference between low [30-60Hz] and high [130-250 Hz] frequency ERPs for low (red) and high (blue) responders. Stars indicate significant (p<.05 corrected) statistical difference. Bottom panel: Correlation between aversion ratings and low minus high ERPs. black dots indicate (cluster corrected) statistical significance at p<.05 level.

	Figure 1: Figure 2:

Acknowledgements:

We thank Arnaud Desvachez and Blanca Marin Bosch for their contribution to experimental work. This work was funded by Swiss National Science Foundation project grant 163040 (A.G.), Fondation Pour l'Audition FPA RD-2020-10 (L.A.). Author contributions: L.H.A designed the experiments. F.S and L.H.A. collected the data. K.D. and L.H.A. performed the analysis. L.H.A. drafted the manuscript. F.S., K.D., S.M., S.S., K.I., A.L.G and L.H.A. corrected and approved the manuscript. Competing interests: The authors declare no competing interests. Correspondence and requests for materials should be addressed to Luc H. Arnal.

Legends


Nicolas Mr

Delbart

Nicolas Jean

Paul Delbart

Jean-Paul Delbart

Laurent Kergoat

Thuy Le Toan

Julien Lhermitte

Ghislain Picard

Determination of phenological dates in boreal regions using Normalized Difference Water

Keywords: phenology, greening-up, senescence, NDWI, short-wave infrared, boreal, Siberia, SPOT-

Monitoring and understanding plant phenology is important in the context of studies of terrestrial productivity and global change. Vegetation phenology such as dates of onsets of greening up and leaf senescence have been determined by remote sensing using mainly the Normalized Difference Vegetation Index (NDVI). In boreal regions, the results suffer from significant uncertainties because of the effect of snow on NDVI. In this paper, SPOT VEGETATION S10 data over Siberia have been analysed to define a more appropriate method. The analysis of time series of NDVI, Normalized Difference Snow Index (NDSI) and Normalized Difference Water Index (NDWI), together with an analysis of in situ phenological records in Siberia, shows that the vegetation phenology can be detected using NDWI, with small effect of snow. In Spring, the date of onset of greening up is taken as the date at which NDWI starts increasing, since NDWI decreases with snowmelt and increases with greening up. In the Fall, the date of onset of leaf coloring is taken as the date at which NDWI starts decreasing, since NDWI decreases with senescence and increases with snow accumulation. The results are compared to the results obtained using NDVI-based methods, taking in situ phenological records as the reference. NDWI gives better estimations of the start of greening up than NDVI (reduced RMSE, bias and dispersions, and higher correlation), whereas it does not improve the determination of the start of Abstract leaf coloring. A map of greening up dates in central Siberia obtained from NDWI is shown for year 2002 and the reliability of the method is discussed.

Introduction

Monitoring and understanding plant phenology, which is the timing of recurrent biological events, is gaining importance in the context of global change. Because of the long-term interest by biologists in phenology, many long time series of in situ phenological records exist, for example from meteorological stations or from "phenological gardens". However, spatially and temporally continuous observations of phenology, required to study the effects of climatic changes on plant phenology, are difficult to generate from sparse ground stations. Consequently such observations have been derived from remote sensing data, especially from the Normalized Difference Vegetation Index (NDVI), an indicator of the density of chlorophyll and leaf tissue [START_REF] Tucker | African land-cover classification using satellite data[END_REF] calculated from the red and near-infrared reflectances from the Advanced Very High Resolution Radiometer (AVHRR). NDVI temporal variations are related to the seasonal changes in the amount of photosynthetic tissues: for typical mid-to-high latitude ecosystems, NDVI increases when the deciduous vegetation greens up in spring and decreases when the vegetation changes color in autumn. Based on these variations, various algorithms have been applied to derive phenological dates from NDVI in previous studies [START_REF] Justice | Analysis of the phenology of global vegetation using meteorological satellite data[END_REF][START_REF] Reed | Measuring phenological variability from satellite imagery[END_REF][START_REF] Lüdeke | The use of satellite NDVI data for the validation of global vegetation phenology models : application to the Frankfurt Biosphere Model[END_REF][START_REF] Moulin | Global-scale assessment of vegetation phenology using NOAA/AVHRR satellite measurements[END_REF][START_REF] White | A continental phenology model for monitoring vegetation responses to interannual climatic variability[END_REF][START_REF] Duchemin | Monitoring phenogical key stages and cycle duration of temperate deciduous forest ecosystems with NOAA/AVHRR data[END_REF].

Studies using AVHRR NDVI have shown that at northern latitudes the onset of greening up progressively earlier by a few days since 1982 [START_REF] Myneni | Increased plant growth in the northern high latitudes from 1981 to 1991[END_REF], Tucker et al. 2001[START_REF] Zhou | Variations in northern activity inferred from satellite data of vegetation index during 1981 to 1999[END_REF][START_REF] Shabanov | Analysis of interannual changes in northern Vegetation activity observed in AVHRR Data from 1981 to 1994[END_REF]. This important result is in agreement with the lengthening of the growing season predicted by phenology models [START_REF] Schwartz | Green-wave phenology[END_REF]) and observed in phenological garden records [START_REF] Menzel | Trends in phenological phases in Europe between 1951 and 1996[END_REF], and with the earlier onset of photosynthesis detected in the atmospheric CO 2 concentration measurements [START_REF] Keeling | Increased activity of northern vegetation inferred from atmospheric CO2 measurements[END_REF][START_REF] Randerson | Increases in early season ecosystem uptake explain recent changes in the seasonal cycle of atmospheric CO2 at high northern latitudes[END_REF].

However, detecting the onset of vegetation greening up from NDVI in boreal regions is difficult as the onset of NDVI increase corresponds to the beginning of snowmelt [START_REF] Moulin | Global-scale assessment of vegetation phenology using NOAA/AVHRR satellite measurements[END_REF]. [START_REF] Shabanov | Analysis of interannual changes in northern Vegetation activity observed in AVHRR Data from 1981 to 1994[END_REF] have analysed the time evolution of red and near-infrared reflectances in northern Europe during spring months from 1981 to 1991, and concluded that the trend observed from NDVI is not related only to an earlier vegetation onset of greening up but also to a reduction in the snow cover extent, which agrees with observations made by [START_REF] Dye | Seasonality and trends of snow-cover, vegetation index, and temperature in Northern Eurasia[END_REF].

Consequently, the effect of snow on the optical remote sensing data -in general NDVI of a given pixel increases when the snow cover fraction in the pixel decreases, and decreases when this fraction increases -must be taken into account when estimating the phenological dates in northern latitudes. To reduce the effects of snow on the NDVI signal, several solutions were proposed. [START_REF] Suzuki | West-east contrast of phenology and climate in northern Asia reavealed using a remotely sensed vegetation index[END_REF] define the onset of greening up as the first date of the year when NDVI exceeds a fixed threshold selected to be higher than NDVI of snow. [START_REF] Shabanov | Analysis of interannual changes in northern Vegetation activity observed in AVHRR Data from 1981 to 1994[END_REF] derive greening up dates only on snow free areas. [START_REF] Zhang | Monitoring vegetation phenology using MODIS[END_REF] replaces the vegetation index value by the last snow free record if snow is detected using a snow spectral index. Then, the date of onset of greening up is the date at which the modified vegetation index starts increasing.

Thanks to the spectral bands available on the recent optical sensors SPOT-VEGETATION (SPOT-VGT) and MODIS, other spectral indices can be calculated and should allow better detection of phenological dates in Northern latitudes. In particular, two interesting spectral indices are:

the Normalized Difference Snow Index (NDSI), a combination of blue and middle infrared bands (Hall et al., 1995) which is an indicator of the snow cover [START_REF] Salomonson | Estimating fractional snow cover from MODIS using the normalized difference snow index[END_REF], the Normalised Difference Water Index (NDWI), a combination of middle infrared and near infrared bands introduced by Gao to assess the water content of vegetation [START_REF] Gao | NDWI -a normalized difference water index for remote sensing of vegetation liquid water from space[END_REF]. NDWI has been used to monitor water stress in semi-arid environment [START_REF] Fensholt | Derivation of a shortwave infrared water stress index from MODIS near-and shortwave infrared data in a semiarid environment[END_REF], to map burnt areas in boreal forest (Fraser and Li, 2002), and to characterise land cover and vegetation type (Xiao et al., 2002 and[START_REF] Boles | Land cover characterization of Temperate East Asia, using multi-temporal VEGETATION sensor data[END_REF]. The latter two papers mentioned the potential of this index for phenological date assessment. [START_REF] Boles | Land cover characterization of Temperate East Asia, using multi-temporal VEGETATION sensor data[END_REF] suggested also that the NDWI decrease preceding the increase due to vegetation development was due to snowmelt. Similarly, the authors suggested that the NDWI increase after senescence was due to snow accumulation.

The aim of this paper is to develop a method using NDWI to derive greening up and leaf coloring dates in boreal regions and to compare the results with those obtained using existing methods based on NDVI. In section 1, we describe the study area located in Central Siberia, the in situ records (dates of beginning of snowmelt and snow cover formation, tree leaf appearance and tree leaf senescence), and the remote sensing data which includes indices derived from SPOT-VGT S10 data. In section 2, we analyse the temporal variations of NDVI, NDSI and NDWI from 1999 to 2002. Together with the in situ records, the variations of NDSI are used to identify the amount of NDVI and NDWI variations due to snow cover variations, and indirectly the part of these variations due to vegetation changes.

A method to estimate the dates of onsets of greening up and leaf coloring from NDWI is developed and described in section 3. Other existing methods using NDVI are described and applied to our datasets to derive the phenological dates. Section 4 presents the comparison between the dates obtained by remote sensing and the in situ observations. The reliability of each method and their ability to detect spatial and temporal phenological variations are assessed.

Finally, the greening up date using NDWI is mapped in Siberia for 2002.

1 Study area and data

Study area

The study area ranges from 80°E to 120°E, and from 50°N to 70°N in Central Siberia (Fig. 1).

The area is intensively studied in the SIBERIA-II project (http://www.siberia2.uni-jena.de/), which aims at fully accounting the carbon budget of vegetation and the greenhouse gas exchanges with the atmosphere, by means of multi-sensor remote sensing and vegetation modelling.

Climate in central Siberia is continental, characterised by very low temperatures in winter and a long snow period, typically from September to May. The main vegetation types include tundra in the North and in mountainous areas, taiga between latitudes 55°N and 65°N, and steppe, croplands and temperate forest in the South. Table 1 gives the area proportions for each cover type according to the Global Land Cover map GLC 2000 [START_REF] Fritz | The Global Land Cover for the Year 2000, Harmonisation, mosaicing and production of the Global Land Cover 2000 database (Beta Version)[END_REF]. The main tree species in the taiga belt belong to the following genus : Larch (Larix), including Siberian Larch (L. siberica), Spruce (Picea), Pine (Pinus) including Siberian Pine (P. sibirica) and Scot Pine (P. silvestris), Birch (Betula), and Poplar (Populus).

In situ data

Phenological data used in this paper have been compiled by the Komarov Botanical Institute, Saint Petersburg. Russian phenological data originate from various networks, such as hydrometeorological stations, agriculture stations and nature reserves [START_REF] Fedotova | Russian phenology: history and present day (poster)[END_REF].

Methodologies for phenological observations are described in [START_REF] Bulygin | Dendrology. Phenological observations on deciduous trees[END_REF] and [START_REF] Elagin | Methodology for collecting and processing data of phenological observations on trees and shrubs[END_REF].

The records include the dates of leaf appearance and autumn leaf coloring for birch, aspen and larch, the dates of snow events such as beginning of snowmelt and beginning of snow cover formation, all obtained through visual observations. One record set dates back to 1902, but many were from more recent years (table 2). All records, from sites 1 to 10 are used for the analysis, whereas only records from sites 1 to 5 are used for validation of remote sensing methods.

Remote Sensing data

SPOT-VGT 10-day-composite data (S10) at the resolution of 0.008928° (about 1km) in platecarree projection were downloaded (http://free.vgt.vito.be/). Data consist in reflectances in spectral bands B0, B2, B3 and SWIR, which correspond respectively to the blue (0.43-0.47µm), red (0.61-0.68 µm), near-infrared (0.78-0.89 µm) and short-wave infrared (1.58-1.75 µm) domains. These reflectances are combined to calculate NDVI, NDWI and NDSI as: NDVI= (B3-B2)/(B3+B2)

(1), NDWI=(B3-SWIR)/(B3+SWIR) (2), and NDSI=(B0-SWIR)/(B0+SWIR) (3).

For each pixel, one acquisition -i.e. one reflectance value for each spectral band-is selected for each 10-day period, based on a Maximum Value Composite (MVC) technique [START_REF] Holben | Characteristics of maximum-value composite images from temporal AVHRR data[END_REF].

The selected acquisition is expected to be the least affected by atmospheric disturbances. The exact acquisition date is provided in the dataset. Duration between two consecutive S10 data ranges from 1 to 20 days.

Despite MVC procedure, S10 data may still be affected by clouds and aerosols. Aerosols affect more strongly the shortest wavelengths (B0, B2) signals. Therefore, NDSI is more sensitive to aerosols than NDVI, which is more sensitive than NDWI.

2 Data analysis

In situ data

For each of the 10 sites, table 3 summarises the mean, standard deviation and temporal range of the key phenological and snow dates. Means and standard deviations of the time intervals between events are also reported. The standard deviation refers to interannual variation, and not to the uncertainties in the in situ observations which are difficult to estimate since they depend significantly on the observers.

Spring events

In all the records, first leaves appear after the beginning of snowmelt (column 6 in table 3). The time interval between the beginning of snowmelt and the appearance of first leaves is of the order of 34 days, and shows significant interannual variations (standard deviation of 17 days on average). Interannual variations of the snowmelt date (column 5) are larger than those of leaf appearance (column 2).

For every site, leaf appearance dates for the different deciduous species are close. In particular birch and larch show a very similar phenology in spring (column 3), with no difference on average (all sites and all years averaged). The average time difference between birch and aspen (column 4) is 3 days (std 4.7). Such a high uniformity of phenological events is also pointed out by [START_REF] Linkosalo | Regularities and patterns in the spring phenology of some boreal trees[END_REF] for Finnish forests, and ascribed to similar responses of trees to climate conditions. The spatial gradient is 34 days in average between the most Southern and Northern sites, whereas at each site mean interannual variabilitions are about 30 days.

Autumn events

The various species start their foliage senescence at various times. In general, aspen starts coloring earlier than birch, followed by larch (not shown in table 3). The time interval between the senescence of the various species is not constant. For example it ranges from 0 to 40 days between aspen and larch for the years under study. In addition, for a given species at a given site, the interannual variabilitions are large: for example, the beginning of birch senescence varies from -20 days to 20 days around the average date. In general, the senescence starts before the first snowfall. However, larch trees start their senescence quite often after the formation of the permanent snow cover.

Analysis of spectral indices from SPOT-VGT

Fig. 2 and3 show typical variations of the three spectral indices in relation with the phenological and snow events. They outline the logic of the algorithm, which is developed in section 3. Fig. 2 presents the 1999 time series of NDSI, NDVI and NDWI for deciduous forest pixels near Krasnoyarsk, a southern site of the study area, with the in situ dates of snowmelt, appearance of birch leaves, beginning of coloring of birch leaves and the beginning of the snow cover formation. Fig. 3 presents the 2000 time series of the three indices for larch pixels near Aksarka, a northern site, with the dates of beginning of snowmelt and first appearance of larch needles.

At the resolution of VGT data, a pixel contains a mixture of snow, green vegetation and non vegetated areas. When applied to optical data, the term snowmelt corresponds to the decay of the area covered by snow in the pixel, while the vegetation greening up corresponds to the increase of the area of green vegetation in the pixel. In autumn, leaf coloring corresponds to the increase of the area of yellow or brown vegetation in the pixel. Last, snow cover formation corresponds to the increase of the snow-covered area in the pixel.

Spring events

In Fig. 2 and 3 NDSI starts decreasing at the beginning of snowmelt. This decrease either continues until or after the appearance of leaves (Fig. 2), until complete snow disappearance, or ends before leaf appearance(Fig. 3) if snowmelt ends before leaf appearance. As expected, NDSI exhibits large variations during summer, due to atmospheric effects.

NDVI increases as the proportion of green leaves increases in the pixel, and is close to zero for a pixel covered by snow. Thus, NDVI increases during both snowmelt and leaf appearance, resulting in either one continuous increase (Fig. 2), or in two distinct increases if snowmelt ends before the leaf appearance (Fig. 3). In the first case, the two events can not be distinguished in the NDVI time series. In the second case, the distinction is difficult since there is no significant change such as a change in the sign of the slope.

NDWI is related to the quantity of water per unit area in the canopy [START_REF] Gao | NDWI -a normalized difference water index for remote sensing of vegetation liquid water from space[END_REF]. It therefore increases during leaf development. However, as opposed to NDVI, it decreases during snowmelt [START_REF] Boles | Land cover characterization of Temperate East Asia, using multi-temporal VEGETATION sensor data[END_REF]. NDSI time series show that greening up may start before or after complete snowmelt. If snowmelt did not melt totally well before leaf appearance (one week or more), NDWI first decreases and then increases (Fig. 2), displaying a "peak" at its minimum. If snowmelt is complete before leaf appearance, then NDWI remains stable during a period which may last between a few days to a few weeks before increasing (Fig. 3).

If greening up occurs during snowmelt, the NDWI decrease due to snowmelt may mask the NDWI increase due to greening up, and NDWI may start increasing later than the actual onset of greening up. However, in situ records show that snowmelt starts in average several weeks before greening up. Thus, in the majority of the cases, a large proportion of snow cover has already melt at the onset of greening up, and the moment at which NDWI starts increasing corresponds to the onset of greening up.

Autumn events

As illustrated in Fig. 2, NDSI remains stable during vegetation coloring, and increases during the formation of the snow cover. NDVI decreases both during vegetation coloring and snow cover formation, and shows only a single and continuous decrease. It is difficult to ascribe specific portions of such a decrease to one or the other event.

NDWI decreases during vegetation coloring. The slow increase in NDWI after the first snow (e.g. on day 270 in Fig. 2c) can be interpreted in two ways:

the snow cover formation is slow at the beginning, there is compensation between the NDWI decrease due to vegetation coloring and the NDWI increase due to snow. In some cases (not shown in Fig. 2 and 3), NDWI does not decrease at all during autumn, which indicates that snowfalls start before or at the beginning of vegetation coloring.

3 Algorithm development

Assessment of noise and data filtering

The noise due to atmosphere and to the sensor must be taken into account in our detection algorithm, and is quantified in this section. For this purpose, we study the temporal variations of the three indices at the end of July when changes due to phenology and snow are expected to be minimum. The distributions of the variations of the three indices between the last two 10-day periods of July are computed. For the three indices, the mean is close to zero, and the standard deviations, which are subsequently considered as noise level, are 0.06 for NDVI, 0.1 for NDSI, and 0.03 for NDWI. Additionally, we assume that the variations due to noise in late spring are the same as those in July.

The MVC processing retain data contaminated by clouds or high aerosol loading if all data within the 10-day period are affected. For the late spring and summer period, the data are filtered out using a threshold of 0.4 on NDSI positive variations between consecutive 10-day periods.

Autumn data were not filtered because successive snow falls and snowmelts have the same signature on NDSI than clouds and aerosols.

Estimation of the date of onset of greening-up

Following the above analysis, a method using NDWI is developed to extract the date of onset of vegetation greening up. For comparison purposes, the dates of onset of greening up are computed by applying three existing methods using NDVI.

Methodology for NDWI time series

During the spring, minimum NDWI corresponds to the vegetation state before greening up. It can be used to detect the onset of vegetation greening up, except if there is a plateau between snowmelt and vegetation greening up, as described in section 2.2.1. In this case, small fluctuations of the NDWI may arbitrarily position the minimum in the plateau period, causing a mismatch between minimum NDWI and the onset of vegetation greening up. Therefore, our algorithm retains the later NDWI record which is less than the minimum NDWI increased by an estimate of the noise level (Fig. 4a and 4b, and eq. 4). In addition for these boreal ecosystems, we assumed that onset of greening up must occur prior to day-of-year 200 (18 th -19 th July).

[ ] ( )

( ) ε + < ∈ = min ) ( 200 , 0 max NDWI t NDWI t t greening (4)
ε was chosen to fulfill the following constraints :

ε should be larger than the noise affecting NDWI time profile.

-ε must be smaller than the first NDWI increase due to vegetation growth.

To assess this increase, statistics of NDWI increase following the minimum NDWI are analysed.

For all the pixels of the study area, the daily rate of increase is calculated for the three 10-day periods following the minimum NDWI. In 80% of non-agricultural pixels, the maximum rate occurs during the first or the second 10-day period after the minimum NDWI. This means that the increase is rapid as soon as it starts. The distribution of the maximum rate for different forest and cover types is shown in Fig. 5. When expressed as a percentage of the overall spring NDWI increase (hereafter referred to as NDWI amplitude), the average rate is 4% per day, except for some agricultural areas. Finally, either the first or the second 10-day period of increase accounts on average for 40% of the amplitude. Such a value appears to be an upper limit for ε.

Based on these considerations, ε was chosen to be 20% of the spring NDWI amplitude. For pixels displaying an amplitude of 0.2 or more (71% of all pixels, see amplitude map in Fig. 6a and its distribution by land cover in Fig. 6b), this means that the detection is not affected by noise induced variations up to 0.04, which is above the estimated standard deviation of noise (estimated as 0.03 NDWI units, see section 3.1). Only 4% of non-agricultural pixels have a maximum increase rate lower than 2% of spring amplitude. For these pixels, there is a siginificant risk to overestimate the date of onset of greening up by one 10-day period or more.

Methodologies for NDVI time series

Several methods have been used to detect the onset of greening up from vegetation index data [START_REF] Lloyd | A phenological classification of terrestrial vegetation cover using shortwave vegetation index imagery[END_REF][START_REF] Reed | Measuring phenological variability from satellite imagery[END_REF][START_REF] Moulin | Global-scale assessment of vegetation phenology using NOAA/AVHRR satellite measurements[END_REF][START_REF] White | A continental phenology model for monitoring vegetation responses to interannual climatic variability[END_REF][START_REF] Schwartz | Assessing satellite-derived start-of-season measures in the conterminous USA[END_REF][START_REF] Zhang | Monitoring vegetation phenology using MODIS[END_REF][START_REF] Suzuki | West-east contrast of phenology and climate in northern Asia reavealed using a remotely sensed vegetation index[END_REF] among other). Among them, two methods were selected and applied to our data. The first method relies on the detection of the beginning of NDVI increase.

Two variants of this method are tested (methods 1a and 1b). In the second method, the onset of greening up is defined as the date NDVI exceeds a threshold. The threshold is either pixel dependant (method 2a), or spatially uniform (method 2b).

Method 1a was developed by [START_REF] Zhang | Monitoring vegetation phenology using MODIS[END_REF][START_REF] Zhang | Climate controls on vegetation patterns in northern mid-and high latitudes inferred from MODIS data[END_REF] to detect spatial variations of phenological dates from the Enhanced Vegetation Index [START_REF] Huete | Overview of the radimetric and biophysical performance of the MODIS vegetation indices[END_REF] calculated from the MODIS nadir BRDF-adjusted reflectances (NBAR) for latitudes higher than 35°N. We apply it on the NDVI time series. It consists in fitting logistic functions y(t) to the NDVI time series. y(t) is defined by:

d bt a e c t y + + + = ) ( 1 ) ( ( 5
)
where a and b are fitted coefficients linked to temporal variations of NDVI, c+d is the maximum NDVI value and d the minimum NDVI value. Two distinct functions are fitted for the greening up and the leaf coloring. Since the phenological cycle may be multimodal, local extrema are detected using a polynomial fit (9 th order) in order to select the main mode. The onset of the NDVI increase is then taken as the date of maximum positive curvature of the ascending function.

In [START_REF] Zhang | Climate controls on vegetation patterns in northern mid-and high latitudes inferred from MODIS data[END_REF], EVI data for which snow is detected (using MODIS quality assurance flags and MODIS Land Surface Temperature) are replaced by snow free data. In this paper, we use a variant of this method to reduce the effect of snow accumulated on the canopy. We assume that the maximum NDVI between January and March (when deciduous species are leaf-less) corresponds to the minimum accumulation of snow on the canopy. NDVIs lower than this maximum are replaced by this maximum. Method 1a consists in estimating the date at which this modified NDVI starts increasing as defined above. Method 1b is the same applied to the nonmodified NDVI.

In the threshold methods [START_REF] Lloyd | A phenological classification of terrestrial vegetation cover using shortwave vegetation index imagery[END_REF][START_REF] White | A continental phenology model for monitoring vegetation responses to interannual climatic variability[END_REF][START_REF] Schwartz | Assessing satellite-derived start-of-season measures in the conterminous USA[END_REF][START_REF] Suzuki | West-east contrast of phenology and climate in northern Asia reavealed using a remotely sensed vegetation index[END_REF], the onset of greening up is assumed to occur when NDVI exceeds a certain value. In [START_REF] White | A continental phenology model for monitoring vegetation responses to interannual climatic variability[END_REF], the date of onset of greening up is defined as the day at which NDVI has increased by half its total amplitude. We modify this method in order to keep the strong noise removal induced by the fit of the logistic function as used in method 1b. Therefore method 2a consists in estimating the date at which the fitted function has increased by half its amplitude

In method 2b, the same threshold value is used for all the pixels. [START_REF] Suzuki | West-east contrast of phenology and climate in northern Asia reavealed using a remotely sensed vegetation index[END_REF] choose a value of 0.2 (AVHRR dataset) considering that a snow covered surface would have a lower NDVI. Since AVHRR NDVI and SPOT-VGT NDVI are not directly comparable due to differences in sensor spectral response, orbits and data processing, we retain the threshold on VGT NDVI, which gives the smaller bias with the in situ data (see section 4.1). This retained NDVI threshold is 0.4.

Extraction of leaf coloring onset date

Method for NDWI

We define the beginning of NDWI decay, which corresponds to leaf coloring or senescence, as the date for which NDWI has decayed by 20% of its autumnal amplitude (other quantities were tried with poorer results).

Method 1 for NDVI Following [START_REF] Zhang | Monitoring vegetation phenology using MODIS[END_REF], a logistic function is fitted to the descending part of the NDVI time series. Senescence is defined as the start of NDVI decay, which is the date of minimum curvature (largest negative curvature value).

Method 2 for NDVI Alternatively, senescence is estimated by the date for which the fitted NDVI has decayed by 20% of its autumnal amplitude. This method follows the threshold approach of [START_REF] White | A continental phenology model for monitoring vegetation responses to interannual climatic variability[END_REF], although in their method no function is fitted to the data.

Methodology of validation

The greening up and senescence dates derived from remote sensing data are compared with in situ data. To avoid mis-registration of in situ site and SPOT-VGT pixel, it is preferable to compare the in situ dates to the dates averaged on pixels of the same vegetation type located in an area (window) surrounding the in situ site. This is justified by the fact that phenology is driven by climate, and thus does not display strong short-scale gradient. Different window sizes were tried, to obtain a good trade-off between a representative number of pixels corresponding to the in situ vegetation class, which favours larger windows, and the proximity to the in situ site, which favours small windows.

Finally the dates of onset of greening up dates using NDVI and NDWI for the deciduous forest pixels (GLC2000 classes 2, 5 and 6) in the window around each in situ measurement point, are compared to the average date of birch and larch leaf appearance (which are very close to each other, as shown in section 2.1) from 1999 to 2002.

Dates of onset of leaf coloring found on pixels of classes 2 and 6 (deciduous broadleaf and mixed forests) in the window are plotted against in situ dates of birch senescence. Dates found on pixels of class 5 (larch) are plotted against the in situ dates of larch senescence, as significant difference exist between the date of senescence of broadleaves and larch as shown in section 2.1.2.

The root mean square error (RMSE), bias, dispersion and correlation coefficient are calculated from all pixels and all years. Correlation coefficient is also calculated from the mean value of all pixels in one window. Formulations of RMSE, bias and dispersion as given in [START_REF] Willmott | Some comments on the evaluation of model performance[END_REF] are given in eq. 6 to 8 : 7), and

N y x RMSE i N i i 2 ) ( - = ∑ (6), N y x Bias N i i i ∑ - = ) ( (
∑ - - - = N i i i Bias y x N Dispersion 2 ) ( 1 1 (8),
where x i is the date obtained by remote sensing, y i is the in situ date for a given pixel in a given year. N is the total number of pixels multiplied by the number of years.

Results

Greening-up

The comparison between retrieved greening up dates and in situ dates are presented in Fig. 7.

Circles represent the mean value of the dates retrieved for the window around a given in situ site.

They are plotted against the average in situ date of appearance of birch and larch leaves. Vertical lines represents the standard deviation. The size of the circles is proportional to the number of the corresponding pixels found in each window. The window size is 0.14° for sites 1, 2, 4 and 5, and 0.34° for site 3. This allows to have a number of deciduous forest pixels ranging from 21 to 32 for all sites, ensuring a balanced importance of all sites in the calculations.

Results from methods 1a and 1b are presented in Fig. 7a and7b. The two methods display similar dispersion and correlation coefficients. However, the biases are different, being 19 days for method 1b, and 6 days for method 1a. This indicates that it is relevant to replace low NDVIs by the highest winter NDVI. For both methods, the bias is higher for the northernmost site, which shows leaf appearance dates later than day of year 150. In [START_REF] Zhang | Climate controls on vegetation patterns in northern mid-and high latitudes inferred from MODIS data[END_REF], snow detection is further completed using an additional filter based on MODIS Land Surface Temperature, which may further reduce the bias and the dispersion.

The results obtained with method 2a are similar to those obtained with method 1a. The bias is larger for the northernmost pixels, whereas in the southern taiga, the remote sensing estimates are relatively closer to the in situ measurements. As the method is based on the detection of the date of half-amplitude on the fitted logistic, the result seems to indicate that in these pixels, the first half of the NDVI increase is due to snowmelt whereas the second half is due to vegetation greening up. However the amplitude proportions depend on land cover or latitude: in the North the half amplitude occurs earlier than the in situ leaf appearance, which occurs between days 155 and 180. This means that more than half of the NDVI increase is due to snowmelt for the Northern sites.

The dates obtained with the threshold method, NDVI method 2b, are highly dispersed, which means that the dates of onset of greening up can not be determined precisely. Interannual variations seem difficult to detect, because of the dispersion. However the bias is small, as we determined the threshold in order to minimize it, and is similar in the Northern and in the Southern sites. This suggests that North-South spatial gradients in the dates of greening up can be determined, as shown in [START_REF] Suzuki | West-east contrast of phenology and climate in northern Asia reavealed using a remotely sensed vegetation index[END_REF].

The best results are obtained with the NDWI method, for South and North sites. There is almost no bias with in situ measurements. The dispersion is much lower than with NDVI-method 2b.

RMSE, bias and dispersion are also calculated year by year, providing values similar to their values calculated with all years. This shows that the spatial variability can be retrieved for any given year. Calculations are also made on each individual site (all years considered) to evaluate the ability of the method to retrieve interranual variations at a given site. At Krasnoyarsk, Taseevo and Aksarka (sites 3, 4 and 5), but not at Ilir (site 1), the correlation between the retrieved dates and the in situ dates is good, indicating that the interannual variations could be retrieved. Unfortunately the small amount of in situ data does not allow to firmly conclude on the robustness of the method to determine the interannual variations.

To assess the effect of the window size on the results, window sizes ranging from 0.06° to 0.36° were tested. Whatever the window size, the dates obtained with the NDWI method display the best overall statistics, with a RMSE ranging from 6.5 days to 9.5 days.

As a result, the determination of the date of onset of greening up using NDWI appears robust, due to the fact that the effects of snow and vegetation on NDWI are in opposite direction, combined with NDWI's relative lack of sensitivity to aerosols.

Leaf coloring

Fig. 8 shows the comparison of the date of onset of leaf coloring obtained using NDWI and NDVI methods 1 and 2 with in situ dates of beginning of senescence. Dates obtained for broadleaf pixels are compared to the in situ date of beginning of senescence of birch (circles), and dates obtained on larch pixels are compared to the in situ dates of beginning of senescence of larch (triangles). Window sizes are similar to those in 4.1.

The three methods using NDVI or NDWI give poor results, with low R 2 and high RMSE (especially for NDWI-based method). Moreover, the spatial variations are not properly observed, as the retrieved dates in the North are very close to the dates in the South, whereas in situ data show significant differences. We examined in the following possible causes of such discrepancies.

-Because of the low amplitude of NDWI variations in autumn, the detection of the beginning of index decay may not be robust. The 20% of the variation amplitude may be too small to be detectable, but other threshold values were tested with even poorer results.

-Lower sun angles in autumn increase atmospheric effects.

-Leaf senescence of different species and different trees start at different times. Thus, extracting a date of onset of leaf coloring for a pixel and comparing it to in situ senescence measurement for one species is not straightforward. Moreover, leaf senescence corresponds to various changes in chemical and structural properties, resulting in changes in spectral responses which may start long before leaf abscission. This is the opposite of spring phenology, where leaf appearance is much more well defined in time.

-The decrease in the water content observed by NDWI may not be linked to the decrease of the chlorophyll content detected by the in situ observations [START_REF] Ceccato | Detecting vegetation leaf water content using reflectance in the optical domain[END_REF]).

-The decrease in NDVI at the end of season "may be related to abiotic factors, such as extended period of cloudiness" [START_REF] Reed | Measuring phenological variability from satellite imagery[END_REF].

In summary, contrary to the dates of onset of greening up, the dates of onset of leaf coloring retrieved from NDVI or NDWI are not accurate. In the following, only the dates of onset of greening up dates are mapped.

Map of greening up dates in Siberia

Fig. 9 shows the dates of onset of greening up obtained using NDWI in year 2002 for the study area represented by the rectangle in Fig. 1. The general pattern is a South to North gradient, with earlier onset of greening up in the South. The exception to the South-North gradient is found in regions corresponding to mountainous areas located in the North, in the North East of the Lake Baikal, and in the South of the area, which display later onset of greening up compared to areas at the same latitude.

Discussions

The dates of onset of greening up of forest pixels with a large proportion of deciduous trees or understorey vegetation are expected to have uncertainties similar to uncertainties represented by RMSE at validation sites. It is likely that a large part of these uncertainties come from the temporal resolution of the VGT S10 data.

In addition to the overall validation based on comparison with independent in situ data, different sources of error in the retrieval algorithms can be assessed:

1) Phenological dates are also calculated for agricultural pixels, although the 20% threshold used in the NDWI algorithm is not optimised for cropland, as NDWI variations show a different behaviour compared to forests (Fig. 5). The corresponding uncertainties are unknown, since no in situ data for crops were available for validation.

2) The detection is dependent on the ε value, which is a proportion of the NDWI increase amplitude in spring. When this amplitude is small, NDWI variations we want to detect are close to the variations due to noise. Fig. 6a shows that the amplitude is low in the three mountainous regions and in the tundra as expected, because of the small amount of vegetation, and also possibly because of the long period of overlap between snowmelt and vegetation greening. The estimated dates at these regions are thus likely to be more affected by noise and likely to be less accurate. Low amplitude is also found in pixels where the proportion of deciduous vegetation is small compared to the proportion of evergreen vegetation. Onset of greening up detection is therefore more uncertain in these regions. A land cover map with deciduous/evergreen proportion may be used to exclude pixels where the deciduous proportion is small. However for all pixels where the amplitude is low, a flag corresponding to a higher level of uncertainty can be assigned directly based on the signal amplitude without any a priori information.

3) If snowmelt and greening up overlap during a long period, NDWI variations with snowmelt and greening up compensate, making NDWI start increasing later than the actual onset of greening up. However, no such case was found in our in situ records. Moreover, because of the compensation between the NDWI decrease and increase, the overall spring amplitude is low.

Such pixels can consequently be assigned a higher level of uncertainty as described above.

4) Last, NDWI time series may be sensitive to water intercepted by leaves and to abrupt increase in soil moisture [START_REF] Madeira | Effect of wetness on turfgrass canopy reflectance[END_REF], Xiao et al. 2002). Given the fact that our algorithm uses NDWI thresholds and does not retain isolated NDWI peaks, possible impacts of these effects on phenology are likely to be limited. The method may also be affected by short term land cover changes, such as floods or fires. Such effects can be removed by using flood and fire information, acquired independently.

Conclusion

In this paper, we have presented an alternative method to the use of NDVI to determine a key phenological date:the date of onset of greening up in the boreal regions. This method relies on the properties of NDWI, calculated from SPOT-VGT S10 data, which is similar to the NDVI but with the red band replaced by the short wave infrared band. NDWI decreases during snowmelt and increases during greening up, thus allowing distinction between snowmelt and vegetation greening up in spring, which is not possible with NDVI. The date of onset of greening up is taken as the date NDWI starts increasing in spring. The comparison with in situ phenological records shows that NDWI is more efficient to estimate the date of onset of greening up than other methods based on NDVI in the Siberian context.

Estimating the date of onset of leaf coloring in autumn proved to be less satisfactory, for both NDVI and NDWI, since poor agreement was found with the in situ dates of beginning of leaf senescence. This result could be due either to the heterogeneity in the dates of beginning of tree senescence within a pixel, or to the complexity of leaf senescence process and associated changes in spectral responses, or to abiotic factors such as atmospheric effects. It may thus not be possible to accurately measure the length of the active period of boreal vegetation through the use of these indices.

The date of onset of greening up was mapped for year 2002 for an area located between 80° and 120° in longitude, and between 50° and 70° in latitude. The accuracy of spring phenology retrieval differs among pixels, as the quality of the determination of the date of onset of greening up depends on properties of each time series, like the spring amplitude and the rate of increase, which change the sensitivity of the detection to the residual noise.

One limitation is the use of 10-day composite data, which prevents the detection of fine variations. In the future, the method could be refined using VGT daily data in order to reduce the error linked to the temporal resolution of the S10 data. The method could also be applied to MODIS data (bands 2 and 6) at the 500 metre resolution.

For long time series however, a method needs to be developed to interpret the NOAA/AVHRR archive, with a focus on the removal of the effect of snow. Our results may help analysing AVHRR Pathfinder data to retrieve leaf appearance dates over the whole Siberia from 1982 onwards. The results could be used to calibrate or validate phenological models incorporated in the dynamic vegetation models [START_REF] Botta | A global prognostic scheme of leaf onset using satellite data[END_REF], and to estimate the impact of phenological changes on terrestrial carbon budget [START_REF] White | Canopy duration has little influence on annual carbon storage in the deciduous broad leaf forest[END_REF].

Figure captions

Figure 1: Extract of the Global Land Cover 2000 map [START_REF] Fritz | The Global Land Cover for the Year 2000, Harmonisation, mosaicing and production of the Global Land Cover 2000 database (Beta Version)[END_REF]. The figure displays 10 classes formed by regrouping the 22 original GLC classes (in brackets) . Rectangle represents the study area. Numbers on the map represent the site numbers of in situ plots as in table 2.

Figure 2 :

 2 Figure 2: Time evolution of NDSI (a), NDVI (b) and NDWI (c) in 1999 of all broadleaf pixels located in the South-East of Krasnoyarsk, between 93°5'E and 93°12'E, and between 55°48'N and 55°54'N. Each of the individual curves corresponds to one pixel. The in situ dates of beginning of snowmelt, birch leaf appearance, beginning of birch senescence and first snow are represented by vertical lines.

Figure 3 :

 3 Figure 3: Time evolution of NDSI (a), NDVI (b) and NDWI (c) in 2000 of all deciduous needleleaf pixels located close to Aksarka, between 67°45'E and 67°51'E, and between 66°24'N and 66°30'N. The in situ dates of beginning of snowmelt and birch leaf appearance are represented by vertical lines.

Figure 4 :

 4 Figure 4: Detection principle: black line is the theoretical NDWI evolution , dots are the 10-daysynthesis acquisitions, A is the actual onset of greening up. t min is the date of minimum NDWI among acquisitions. Threshold level is indicated. Circled dot is the last acquisition below the threshold and is the acquisition retained as the onset of greening up. Two main cases are represented: a) NDWI decreases with snowmelt and increases with vegetation development. b) NDWI plateaus before increasing.

Figure 5 :

 5 Figure 5: Distribution of the maximum NDWI increase rate calculated for all pixels of the study area, by land cover, in 2002.

Figure 6 :

 6 Figure 6: a) Map of NDWI increase amplitude during spring in 2002. b) Distribution of NDWI increase amplitude during spring, by land cover, for the study area, in 2002.

Figure 7

 7 Figure 7, a, b,c, d, e: Comparison of in situ measurements of leaf appearance dates with estimates derived from the five remote sensing methods. a) start of NDVI increase after replacement of low NDVIs (NDVI method 1a), b) start of NDVI increase keeping low NDVI data (NDVI method 1b), c) NDVI half-peak (NDVI method 2a), d) date at which NDVI exceeds the 0.4 threshold (NDVI method 2b), e) start of NDWI increase. For each site and each year, a circle represent the mean value of deciduous forest pixels. Symbols size is proportional to the number of pixels in each window. Vertical bars are the standard deviations.

Figure 8 a

 8 Figure 8 a,b,c: Same as Fig. 7, but for the comparison of onset of leaf coloring dates derived from remote sensing data with the in situ senescence dates. For each site and each year, circles represent the mean value of pixels classified as broadleaf or mixed type forest. Triangles represent the values of deciduous needleleaf forest. a) start of NDVI decrease (NDVI method 1), b) start of NDVI decrease (NDVI method 2), c) start of NDWI decrease.

Figure 9 :

 9 Figure 9: Greening-up onset dates (day of year) estimated using the NDWI method for the Siberian study area in 2002.

Figure 3

 3 Figure 3

  Figure 3

Figure 5

 5 Figure 5

  Figure 5

  Figure 6 :

Figure 6 Figure 7

 67 Figure 6

Figure 7 Figure 8

 78 Figure 7

Figure 8 Figure 9

 89 Figure 8

Figure 9

 9 Figure 9

Acknowledgements

Vegetation phenology data were provided by Dr. Violetta Fedotova, of the Phenology research group at the Institute of Botany, Komarov Russian Academy of Sciences, St. Petersburg. This work was conducted in the framework of the Siberia II project ( Multi-sensor concepts for Greenhouse Gas Accounting of Northern Eurasia), EC Framework 4 Contract: EVG1-CT-2001-00048. Many thanks to Niall Hanan for his useful comments.

Table 1: Surface proportions of each land cover type in the study area. Original GLC class numbers are indicated (in brackets).

Land cover type

Surface proportion Deciduous needleleaf forest [START_REF] Boles | Land cover characterization of Temperate East Asia, using multi-temporal VEGETATION sensor data[END_REF] 34.3% Evergreen needleleaf forest (4) 16.5% Mixed leaf Forest (6) 10.7% Deciduous shrubs (12) 7.6% Herbaceous and sparse shrubs (13,14) 7.3% Agriculture (16 , 17, 18) 6.4% Deciduous broadleaf forest (2) 5.5% Bare soils or barrens (19,22) 2.1% Water bodies (20) 2.6% Evergreen shrubs (11) 0.4% Tables Table 2: Available in situ phenology records : location names, geographical coordinates, and temporal coverages.

Location Coordinates

Temporal coverages of the records 1 Ilir 53°30'N, 100°42 'E 1927-1930, 1939-1940, 1979-2002 2 Volchno-Burlingskoe 54°N, 81°E 198954°N, 81°E -200254°N, 81°E 3 Krasnoyarsk 56°N, 93°E 190254°N, 81°E -198254°N, 81°E , 199854°N, 81°E -200054°N, 81°E , 2002 4 4 Taseevo 57°11'N, 94°50 'E 1960-2002 5 Aksarka 66°30'N, 67°49 'E 1924-1936, 1951-2002 6 Tachtip 52°44 'N, 90°E 1935-1959, 1960-1970 (GLC 7,8,9,10) Water bodies (GLC 20)