Guillaume Lassalle
email: guillaume.lassalle@onera.fr

Sophie Fabre

Anthony Credoz

Remy Hedacq

Georges Bertoni

Dominique Dubucq

Arnaud Elger

Rémy Hédacq

Application of PROSPECT for estimating total petroleum hydrocarbons in contaminated soils from leaf optical properties

Keywords: soil contamination, leaf optical properties, PROSPECT, pigment, total petroleum hydrocarbons

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

Introduction

In the last decades, increasing efforts have been made to improve the monitoring of past and current industrial activities [START_REF] Datta | A review on different pipeline fault detection methods[END_REF][START_REF] Jha | Advances in Remote Sensing for Oil Spill Disaster Management: State-of-the-Art Sensors Technology for Oil Spill Surveillance[END_REF][START_REF] Shukla | Application of robotics in onshore oil and gas industry-A review Part I[END_REF]. The assessment of health and environmental risks remains a major issue as it requires efficient detection of contamination deriving from production facility failures and bad practices [START_REF] Barraza | Distribution, contents and health risk assessment of metal(loid)s in small-scale farms in the Ecuadorian Amazon: An insight into impacts of oil activities[END_REF][START_REF] Correa Pabón | Reflectance and imaging spectroscopy applied to detection of petroleum hydrocarbon pollution in bare soils[END_REF][START_REF] Perrodin | Ecological risk assessment of urban and industrial systems: A review[END_REF]. For this purpose, new approaches using hyperspectral spectroscopy have recently emerged in the oil and gas industry [START_REF] Khanna | Comparing the potential of multispectral and hyperspectral data for monitoring oil spill impact[END_REF][START_REF] Scafutto | Characterization of mineral substrates impregnated with crude oils using proximal infrared hyperspectral imaging[END_REF]. In exploiting the light absorption features of oil, these approaches have proved efficient for detecting and quantifying Total Petroleum Hydrocarbons (TPH) from the spectral signature of soils [START_REF] Correa Pabón | Spectroscopic characterization of red latosols contaminated by petroleum-hydrocarbon and empirical model to estimate pollutant content and type[END_REF][START_REF] Scafutto | Hyperspectral remote sensing detection of petroleum hydrocarbons in mixtures with mineral substrates: Implications for onshore exploration and monitoring[END_REF]. However, their applicability remains limited in vegetated areas, since only few radiations get through the foliage.

As an alternative, recent studies have proposed exploiting the spectral signature of vegetation to detect reflectance changes induced by the presence of oil in soils [START_REF] Arellano | Detecting the effects of hydrocarbon pollution in the Amazon forest using hyperspectral satellite images[END_REF][START_REF] Lassalle | Assessing Soil Contamination Due to Oil and Gas Production Using Vegetation Hyperspectral Reflectance[END_REF][START_REF] Noomen | Spectral and spatial indicators of botanical changes caused by long-term hydrocarbon seepage[END_REF][START_REF] Sanches | Assessing the impact of hydrocarbon leakages on vegetation using reflectance spectroscopy[END_REF].

The biophysical and biochemical parameters driving leaf optical properties are likely to be affected when vegetation is exposed to TPH [START_REF] H. Ur R. Athar | Influence of sub-lethal crude oil concentration on growth, water relations and photosynthetic capacity of maize (Zea mays L.) plants[END_REF][START_REF] Balliana | Development of Canavalia ensiformis in soil contaminated with diesel oil[END_REF][START_REF] Baruah | Effect of crude oil contamination on the chlorophyll content and morpho-anatomy of Cyperus brevifolius (Rottb.) Hassk[END_REF]. This includes alterations in leaf anatomy (e.g. cuticle and parenchyma thickness), pigment and water contents [START_REF] H. Ur R. Athar | Influence of sub-lethal crude oil concentration on growth, water relations and photosynthetic capacity of maize (Zea mays L.) plants[END_REF][START_REF] Balliana | Development of Canavalia ensiformis in soil contaminated with diesel oil[END_REF][START_REF] Han | The effects of petroleum-contaminated soil on photosynthesis or Amorpha fruticosa seedlings[END_REF]. These effects mainly result from the modification of soil properties and the reduction of root uptake capacities [START_REF] Balasubramaniyam | Scanning electron microscopic investigations of root structural modifications arising from growth in crude oil-contaminated sand[END_REF][START_REF] Kisic | The effects of drilling fluids and crude oil on some chemical characteristics of soil and crops[END_REF][START_REF] Nie | Plants' use of different nitrogen forms in response to crude oil contamination[END_REF].

They vary according to the species, the level of contamination and the duration of exposure.

Alterations in leaf biophysical and biochemical parameters lead to reflectance changes, which help discriminate between healthy and affected vegetation [START_REF] Lassalle | Assessing Soil Contamination Due to Oil and Gas Production Using Vegetation Hyperspectral Reflectance[END_REF][START_REF] Rosso | Reflectance properties and physiological responses of Salicornia virginica to heavy metal and petroleum contamination[END_REF]. Previous work focused on pigment alteration, since it is responsible for the increase in reflectance typically observed in the visible (VIS) when plants are exposed to TPH [START_REF] Arellano | Field spectroscopy and radiative transfer modelling to assess impacts of petroleum pollution on biophysical and biochemical parameters of the Amazon rainforest[END_REF][START_REF] Lassalle | Hyperspectral signature analysis of three plant species to long-term hydrocarbon and heavy metal exposure[END_REF][START_REF] Lassalle | Detection and discrimination of various oil-contaminated soils using vegetation reflectance[END_REF]. The amplitude of this increase is influenced by the level of contamination [START_REF] Sanches | Assessing the impact of hydrocarbon leakages on vegetation using reflectance spectroscopy[END_REF][START_REF] Zhu | Soil TPH concentration estimation using vegetation indices in an oil polluted area of eastern China[END_REF]. This suggests that leaf pigments could be used for detecting and quantifying oil contamination in vegetated areas. However, no method has been yet proposed in literature. Previous studies focused on quantifying heavy metal contamination deriving from mining or agricultural activities, a very different case of soil contamination [START_REF] Shi | Monitoring arsenic contamination in agricultural soils with reflectance spectroscopy of rice plants[END_REF][START_REF] Zhou | Evaluating metal effects on the reflectance spectra of plant leaves during different seasons in post-mining areas, China[END_REF].

Leaf pigments undergo significant changes throughout the year, depending on the season [START_REF] Archetti | Unravelling the evolution of autumn colours: an interdisciplinary approach[END_REF][START_REF] Demarez | Seasonal variation of leaf chlorophyll content of a temperate forest. Inversion of the PROSPECT model[END_REF][START_REF] Stylinski | Seasonal patterns of reflectance indices, carotenoid pigments and photosynthesis of evergreen chaparral species[END_REF]. Seasonal variations in pigment contents are therefore likely to hide the alterations induced by TPH. This can cause confusion and currently represents a major limitation for TPH detection under natural conditions. All the more so that not all species are good indicators of TPH, and their sensitivity to these contaminants is likely to affect the relative influence of TPH and seasons on leaf pigments and reflectance [START_REF] Lassalle | Assessing Soil Contamination Due to Oil and Gas Production Using Vegetation Hyperspectral Reflectance[END_REF][START_REF] Credoz | Experimental study of hyperspectral responses of plants grown on mud pit soil[END_REF][START_REF] Emengini | Comparative Analysis of Spectral Responses of Varied Plant Species to Oil Stress[END_REF].

The species established around industrial facilities are particularly tolerant to high levels of TPH [START_REF] Arellano | Field spectroscopy and radiative transfer modelling to assess impacts of petroleum pollution on biophysical and biochemical parameters of the Amazon rainforest[END_REF][START_REF] Arellano | Plant family-specific impacts of petroleum pollution on biodiversity and leaf chlorophyll content in the Amazon rainforest of Ecuador[END_REF]. This makes it difficult to detect and estimate TPH, as this implies being able to differentiate small alterations in pigment contents using leaf optical properties. Many approaches have been developed in hyperspectral spectroscopy for estimating leaf pigment contents [START_REF] Blackburn | Quantifying chlorophylls and carotenoids at leaf and canopy scales, Remote Sens[END_REF][START_REF] Gitelson | Detection of red edge position and chlorophyll content by reflectance measurements near 700 nm[END_REF][START_REF] Sims | Relationships between leaf pigment content and spectral reflectance across a wide range of species, leaf structures and developmental stages[END_REF][START_REF] Verrelst | Optical remote sensing and the retrieval of terrestrial vegetation bio-geophysical properties -A review[END_REF].

The physically-based model PROSPECT [START_REF] Jacquemoud | PROSPECT: A model of leaf optical properties spectra[END_REF] proved to be one of the most reliable approaches across a wide range of ecological contexts [START_REF] Jacquemoud | PROSPECT + SAIL models: A review of use for vegetation characterization[END_REF][START_REF] Jiang | Estimation of leaf traits from reflectance measurements: Comparison between methods based on vegetation indices and several versions of the PROSPECT model[END_REF]. The model has been recently improved and now enables pigment separation [START_REF] Féret | PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments[END_REF][START_REF] Féret | Towards modeling leaf optical properties through a complete lifecycle[END_REF]. This is of major importance for detecting and estimating TPH, since not all pigments are affected in the same way by contamination [START_REF] Lassalle | Hyperspectral signature analysis of three plant species to long-term hydrocarbon and heavy metal exposure[END_REF][START_REF] Lassalle | Detection and discrimination of various oil-contaminated soils using vegetation reflectance[END_REF]. PROSPECT could therefore be used for detecting and estimating TPH-induced changes in leaf pigment contents and consequently for determining the level of soil contamination [START_REF] Arellano | Field spectroscopy and radiative transfer modelling to assess impacts of petroleum pollution on biophysical and biochemical parameters of the Amazon rainforest[END_REF][START_REF] Arellano | Plant family-specific impacts of petroleum pollution on biodiversity and leaf chlorophyll content in the Amazon rainforest of Ecuador[END_REF].

This study proposes to assess the potential of leaf optical properties for predicting TPH concentrations in the soil, under natural temperate conditions. The proposed method is based on the retrieval of leaf pigment contents using the PROSPECT model, which are thereafter linked to TPH concentrations. The study was carried out on a range of plant species over a complete seasonal cycle in order to identify the most suitable conditions for predicting TPH. The proposed method was compared to other approaches typically used in similar studies.

Materials and Methods

Study site and species

The study was carried out on an industrial brownfield located in a temperate region. The site was subject to intensive oil and gas activities for 20 years, and production residues contaminated by TPH and heavy metals (HM) accumulated in the soil. Since then, it has been colonized by vegetation. The study focused on the predominant species, Rubus fruticosus L., which is generally well-established on industrial sites in temperate regions. Moreover, this species has already proved great potential in hyperspectral spectroscopy for detecting oil-induced pigment alterations [START_REF] Lassalle | Assessing Soil Contamination Due to Oil and Gas Production Using Vegetation Hyperspectral Reflectance[END_REF][START_REF] Lassalle | Detection and discrimination of various oil-contaminated soils using vegetation reflectance[END_REF][START_REF] Credoz | Experimental study of hyperspectral responses of plants grown on mud pit soil[END_REF]. Other species found on the brownfield, including Quercus pubescens Wild., Populus x canadensis Moench., Cornus sanguinea L. and Buddleja davidii Franch., were also tested for their sensitivity to TPH. The same species were also studied on an uncontaminated site (hereafter described as "control site") for comparison purposes. This site was located a few kilometers away from the brownfield and exhibited similar soil properties.

Field sampling and measurements

A mesh of 23 vegetated plots was defined, covering most of the brownfield. Plots consisted in areas of one meter in diameter and were spaced by 10 meters. The upper soil layer (0-20 cm) was analyzed for contaminants at the center of each plot. Analyses revealed C10-C40 TPH concentrations ranging from 16 to 77 g.kg -1 , the main contributors being in of the C21-C40 fraction (dense hydrocarbons). These concentrations correspond to a moderate to high level of exposure for plants, in comparison with other studies [START_REF] Balliana | Development of Canavalia ensiformis in soil contaminated with diesel oil[END_REF][START_REF] Zhu | Soil TPH concentration estimation using vegetation indices in an oil polluted area of eastern China[END_REF]. High concentrations of HM and Polycyclic Aromatic Hydrocarbons (PAH) were also found. The detailed soil analyses are given in the Supporting Information.

Leaf samples were collected on site. For each species found in each plot, three young leaves were randomly sampled on different individuals. Their reflectance was measured in the [400:2500] nm domain using an ASD FieldSpec 4 Hi-Res spectroradiometer (Malvern Panalytical, Malvern, UK) and a leaf-clip with an internal light source. These measurements avoid the influence of leaf orientation and density and the effects of soil and atmosphere. Leaf radiance was acquired on a small black background panel and converted to reflectance as described in [START_REF] Lassalle | Detection and discrimination of various oil-contaminated soils using vegetation reflectance[END_REF]. Measurements were carried out from June 2017 to October 2018, on a seasonal basis (Tab. 1). This included six different sampling dates, with samples taken in contrasted environmental conditions and at different phenological stages. On the same dates, additional spectral signatures were acquired on the control site. There were no TPH detected in the control soil and low HM concentrations corresponded to the geochemical background. Only B. davidii and R. fruticosus were found on the brownfield in March 2018. A total of 1620 spectral signatures had been collected on the brownfield by the end of the study. 152 leaves were used for determining water (Cw), dry matter (Cm) and pigment contents, following the procedures described in [START_REF] Yilmaz | Remote sensing of vegetation water content from equivalent water thickness using satellite imagery[END_REF][START_REF] Diepens | Toxicity effects of an environmental realistic herbicide mixture on the seagrass Zostera noltei[END_REF]. Pigments were analyzed by High Pressure Liquid Chromatography (HPLC). The chlorophyll content (Cab) was expressed as the sum of chlorophylls a and b, and the carotenoid content (Cxc) as the sum of β-carotene and xanthophylls.

These analyses are detailed in the Supporting Information. They served to validate PROSPECT inversion outputs, as described in section 2.3.

PROSPECT inversion and validation

The PROSPECT model was developed for simulating the hemispherical reflectance and transmittance of leaves based on their known structural and biochemical parameters (listed in Tab.

2). [START_REF] Jacquemoud | PROSPECT: A model of leaf optical properties spectra[END_REF]. Inversion of the model allows these parameters to be retrieved from the leaf spectral signatures [START_REF] Jiang | Estimation of leaf traits from reflectance measurements: Comparison between methods based on vegetation indices and several versions of the PROSPECT model[END_REF][START_REF] Féret | PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments[END_REF]. Although PROSPECT was initially calibrated for measurements carried out with an integrating sphere, those performed with a leaf-clip can also be used [START_REF] Li | Retrieval of leaf biochemical parameters using PROSPECT inversion: A new approach for alleviating ill-posed problems[END_REF][START_REF] Sonobe | Using spectral reflectance to estimate leaf chlorophyll content of tea with shading treatments[END_REF]. The inversion consists in finding the set of structural and biochemical parameters -denoted by the vector θ -that best simulates the leaf spectral signature measured [START_REF] Féret | PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments[END_REF]. This is achieved by successively testing sets of parameters and retaining the one that minimizes the cost function J, defined as:

𝐽(θ) = ∑(R meas (λ) -R sim (λ, θ)) 2 λ max λ min (1)
where Rsim and Rmeas stand for the simulated and measured reflectance, respectively [START_REF] Li | Retrieval of leaf biochemical parameters using PROSPECT inversion: A new approach for alleviating ill-posed problems[END_REF][START_REF] Sonobe | Using spectral reflectance to estimate leaf chlorophyll content of tea with shading treatments[END_REF][START_REF] Sun | Analyzing the performance of PROSPECT model inversion based on different spectral information for leaf biochemical properties retrieval[END_REF]. The cost function is computed in the [λmin:λmax] spectral range [START_REF] Féret | PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments[END_REF]. Since only leaf reflectance was measured in our study, we simplified the original function that included transmittance [START_REF] Jacquemoud | PROSPECT: A model of leaf optical properties spectra[END_REF].

PROSPECT inversions were first used to estimate the biochemical parameters of the 152 leaves sampled for analysis. Anthocyanins were recently added to the latest version of the model (PROSPECT-D), but their concentrations remained below the detection limit for all the species in our study. We therefore used the PROSPECT-5 version to perform the inversions as described in [START_REF] Qiu | Improving the PROSPECT model to consider anisotropic scattering of leaf internal materials and its use for retrieving leaf biomass in fresh leaves[END_REF][START_REF] Wang | Applicability of the PROSPECT model for estimating protein and cellulose + lignin in fresh leaves[END_REF][START_REF] Berger | Evaluation of the PROSAIL model capabilities for future hyperspectral model environments: A review study[END_REF]. The six leaf parameters (N, Cab, Cxc, Cbp, Cw and Cm) were simultaneously determined through minimization of the cost function (1) over the entire reflective domain (400-2500 nm). (N corresponds to the leaf structure parameter and mainly influences reflectance in the near-infrared domain (750 -1300 nm) [START_REF] Jacquemoud | PROSPECT: A model of leaf optical properties spectra[END_REF][START_REF] Ourcival | Exploring the relationships between reflectance and anatomical and biochemical properties in Quercus ilex leaves[END_REF][START_REF] Slaton | Estimating near-infrared leaf reflectance from leaf structural characteristics[END_REF].) Inversions were performed using an iterative optimization approach. Many methods, such as look-up tables and conventional optimization algorithms have been tested for this purpose in previous studies [START_REF] Darvishzadeh | Inversion of a radiative transfer model for estimating vegetation LAI and chlorophyll in a heterogeneous grassland[END_REF][START_REF] Zarco-Tejada | Needle chlorophyll content estimation through model inversion using hyperspectral data from boreal conifer forest canopies[END_REF]. However, these methods can sometimes pose problems. For example, local minima and algorithm nonconvergence lead to biased or inaccurate estimation of biochemical parameters [START_REF] Altmann | Supervised nonlinear spectral unmixing using a postnonlinear mixing model for hyperspectral imagery[END_REF][START_REF] Jacquemoud | Extraction of Vegetation Biophysical Parameters by Inversion of the PROSPECT + SAIL Models on Sugar Beet Canopy Reflectance Data. Application to TM and AVIRIS Sensors[END_REF]. To avoid this, we used the differential evolution (DE) algorithm [START_REF] Storn | Differential Evolution -A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces[END_REF]. The DE combines the advantages of evolutionary and genetic algorithms, which proved to be good candidates for solving inversion problems in remote sensing [START_REF] Das | Recent advances in differential evolution-An updated survey[END_REF][START_REF] Ghosh | Self-adaptive differential evolution for feature selection in hyperspectral image data[END_REF][START_REF] Zhong | An Adaptive Differential Evolution Endmember Extraction Algorithm for Hyperspectral Remote Sensing Imagery[END_REF]. To avoid extensive computational time, the parameters tested for optimization were constrained to ranges set out in Table 2. For each sample, the best set of leaf parameters was retained. Inversion outputs were validated in two ways. The Root Mean Squared Error (RMSE), the BIAS and the Standard Error of Prediction Corrected from the bias (SEPC) were used to assess the fit between the simulated and the measured spectral signatures on a wavelength-by-wavelength basis [START_REF] Féret | PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments[END_REF]. The estimated parameters were compared to biochemical analysis using the same metrics and the R 2 . After validation of the model's performances, we applied the inversion of PROSPECT to the 1620 leaves sampled in order to determine their pigment contents. The same procedure was applied to the leaves sampled on the control site.

Variability of leaf pigment contents and TPH estimation

The study then focused on the chlorophyll and carotenoid contents retrieved from PROSPECT inversions. Their seasonal variability was examined for each species using one-way ANOVA and Tukey post-hoc tests to look for significant differences among sampling dates. Pigment contents were then linked to soil TPH concentrations in order to assess their ability to predict the level of contamination. For this purpose, data were divided into 50% train and 50% test sets using the Kennard-Stone algorithm, which ensures that representative samples are selected [START_REF] Kennard | American Society for Quality[END_REF][START_REF] Shi | Estimation of arsenic in agricultural soils using hyperspectral vegetation indices of rice[END_REF].

Univariate regression models were adjusted between pigment contents and TPH concentrations on the train set [START_REF] Zhu | Soil TPH concentration estimation using vegetation indices in an oil polluted area of eastern China[END_REF][START_REF] Shi | Estimation of arsenic in agricultural soils using hyperspectral vegetation indices of rice[END_REF]. The models were then applied to the test set to predict TPH concentrations.

This procedure was carried out on the five species successively, on each sampling date. The accuracy of the predictions was assessed by calculating the R 2 and the RMSE on the train and test sets. In addition, we computed the BIAS, the SEPC and the Residual Predictive Deviation (RPD) on the test sets to complete the accuracy assessment. The RPD is defined as the ratio of the standard deviation of the measured TPH concentration to the RMSE computed between the predicted and measured concentrations [START_REF] Shi | Monitoring arsenic contamination in agricultural soils with reflectance spectroscopy of rice plants[END_REF]. The interpretation of RPD values strongly depends on the context [START_REF] Shi | Monitoring arsenic contamination in agricultural soils with reflectance spectroscopy of rice plants[END_REF][START_REF] Kawamura | Vis-NIR spectroscopy and PLS regression with waveband selection for estimating the total C and N of paddy soils in Madagascar[END_REF]. RPD values above 2 and 3 are often considered as good and excellent predictions of TPH, respectively, and values below 2 as poor ones. A comprehensive flowchart of the method is presented in Figure 1.

In order to assess the accuracy of our approach, the results were compared to those obtained with approaches directly linking TPH concentrations to leaf reflectance in the VIS, including pigmentrelated vegetation indices (listed in [START_REF] Erudel | Criteria Comparison for Classifying Peatland Vegetation Types Using In Situ Hyperspectral Measurements[END_REF]) and partial least square regression (PLSR) with reflectance transformation (derivatives, continuum removal, etc.) [START_REF] Lassalle | Assessing Soil Contamination Due to Oil and Gas Production Using Vegetation Hyperspectral Reflectance[END_REF][START_REF] Shi | Monitoring arsenic contamination in agricultural soils with reflectance spectroscopy of rice plants[END_REF]. These approaches usually perform well for tracking reflectance changes induced by contaminants [11-12, 25, 66].

The approaches were compared based on the RMSE and RPD values obtained on the test sets. The spectral signatures measured were used to retrieve leaf structural and biochemical parameters using PROSPECT-5 inversions. The estimated parameters were compared to chemical analyses, and the simulated spectral signatures were compared to the measured ones.

The estimated pigment contents (Cab and Cxc) obtained for each season were compared, using ANOVA and Tukey post-hoc tests, and were then used to predict TPH concentrations with the help of univariate models. The predicted TPH concentrations were then compared to those obtained from soil analyses.

Results and Discussion

Measured spectral signatures

Since the five species are well-established on the brownfield, they are naturally tolerant to high levels of TPH contamination [START_REF] Lassalle | Detection and discrimination of various oil-contaminated soils using vegetation reflectance[END_REF][START_REF] Credoz | Experimental study of hyperspectral responses of plants grown on mud pit soil[END_REF][START_REF] Benyó | Physiological and molecular responses to heavy metal stresses suggest different detoxification mechanism of Populus deltoides and P. x canadensis[END_REF]. Consequently, the leaves did not exhibit marked visible stress symptoms, and their spectral signatures strongly differed from those observed for species highly affected by oil (Fig. 2) [START_REF] Sanches | Assessing the impact of hydrocarbon leakages on vegetation using reflectance spectroscopy[END_REF][START_REF] Rosso | Reflectance properties and physiological responses of Salicornia virginica to heavy metal and petroleum contamination[END_REF][START_REF] Emengini | Early detection of oil-induced stress in crops using spectral and thermal responses[END_REF]. The spectral signatures did however vary temporally and among species during the study. The variability among the five species highlighted natural differences in leaf anatomy and biochemistry and in sensitivity to TPH exposure [START_REF] Arellano | Plant family-specific impacts of petroleum pollution on biodiversity and leaf chlorophyll content in the Amazon rainforest of Ecuador[END_REF][START_REF] Sims | Relationships between leaf pigment content and spectral reflectance across a wide range of species, leaf structures and developmental stages[END_REF][START_REF] Knapp | Variability in leaf optical properties among 26 species from a broad range of habitats[END_REF]. When considering a single species, the spectral signatures also varied temporally with season changes, especially in the VIS and the red-edge regions. The summer of 2017 was characterized by repeated heat waves (June), and reflectance was particularly low for all species, varying between 550 and 680 nm and affecting the red-edge region on the following sampling dates, an indication of changes in leaf pigment contents [START_REF] Zhu | Monitoring plant response to phenanthrene using the red edge of canopy hyperspectral reflectance[END_REF][START_REF] Sanches | Unravelling remote sensing signatures of plants contaminated with gasoline and diesel: An approach using the red edge spectral feature[END_REF]. Likewise, reflectance in the near-and short-wave infrared regions was also affected by the seasons, a sign of changes in leaf anatomy and water content [START_REF] Ourcival | Exploring the relationships between reflectance and anatomical and biochemical properties in Quercus ilex leaves[END_REF][START_REF] Clevers | Using spectral information from the NIR water absorption features for the retrieval of canopy water content[END_REF][START_REF] Foley | Foliar spectral properties following leaf clipping and implications for handling techniques[END_REF]. All species do not share the same sensitivity to environmental changes from season to season [START_REF] Demarez | Seasonal variation of leaf chlorophyll content of a temperate forest. Inversion of the PROSPECT model[END_REF][START_REF] Stylinski | Seasonal patterns of reflectance indices, carotenoid pigments and photosynthesis of evergreen chaparral species[END_REF]. Only small seasonal differences in the leaf reflectance of R. fruticosus were observed in relation to the other species. In comparison, C. sanguinea underwent substantial changes in October, exhibiting dry blackish leaves with very low VIS reflectance. This unknown symptom differed from the autumn leaf coloration described on Cornus species [START_REF] Feild | Why leaves turn red in autumn. The role of anthocyanins in senescing leaves of red-osier dogwood[END_REF].

The spectral signatures observed on the brownfield are thus the result of leaf anatomy and biochemistry, which are strongly influenced by the species' sensitivities to TPH and the seasons.

This explains some of the differences observed among the species and is a significant factor in the prediction of TPH, as described in sections 3.2. to 3.4. PROSPECT therefore proved efficient for tracking seasonal changes in leaf biochemistry, as suggested in previous studies [START_REF] Demarez | Seasonal variation of leaf chlorophyll content of a temperate forest. Inversion of the PROSPECT model[END_REF].

Evaluation of PROSPECT inversions outputs

Retrieved biochemical parameters

Chemical analyses performed on the 152 leaves indicated that PROSPECT accurately estimated biochemical parameters, regardless of the species. An illustration is provided for R. fruticosus in Figure 3. The RMSE between the simulated and measured signatures was particularly low in the VIS (< 1.3%), and resulted in accurate retrieval of chlorophyll and carotenoid contents.

Carotenoids were more difficult to estimate because of predominant absorption of chlorophylls in the VIS [START_REF] Féret | PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments[END_REF]. However, the RMSE remained below or equal to 1.82 µg.cm -2 for these pigments, which still indicated good accuracy. As suggested in previous work, separating carotenoids into carotenes and xanthophylls in PROSPECT might improve estimations, but this is not yet possible in the current version of the model [START_REF] Féret | PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments[END_REF][START_REF] Zhang | An extended PROSPECT: Advance in the leaf optical properties model separating total chlorophylls into chlorophyll a and b[END_REF]. Leaf water and dry matter contents were also correctly estimated (Fig. 3). All these results were consistent with those obtained on other datasets [START_REF] Li | Retrieval of leaf biochemical parameters using PROSPECT inversion: A new approach for alleviating ill-posed problems[END_REF][START_REF] Sonobe | Using spectral reflectance to estimate leaf chlorophyll content of tea with shading treatments[END_REF][START_REF] Sun | Analyzing the performance of PROSPECT model inversion based on different spectral information for leaf biochemical properties retrieval[END_REF][START_REF] Liu | A novel principal component analysis method for the reconstruction of leaf reflectance spectra and retrieval of leaf biochemical contents[END_REF][START_REF] Zhang | An extended PROSPECT: Advance in the leaf optical properties model separating total chlorophylls into chlorophyll a and b[END_REF].

Figure 3. Comparisons of the predicted and measured biochemical parameters of R. fruticosus leaves. Water and dry matter contents were expressed in mg.cm -2 for graphical convenience.

Seasonal and interspecific variability of leaf pigments

The chlorophyll contents estimated throughout study ranged from 9.3 to 90 µg.cm -2 on the brownfield, and were higher for all species on the control site with the exception of B. davidii.

This species was not affected by exposure to TPH and is therefore not a suitable indicator of oil contamination. Detailed results can be found in the Supporting Information section. The reduction in chlorophyll content observed for the other species is commonly described on oil-exposed vegetation [START_REF] Balliana | Development of Canavalia ensiformis in soil contaminated with diesel oil[END_REF][START_REF] Arellano | Plant family-specific impacts of petroleum pollution on biodiversity and leaf chlorophyll content in the Amazon rainforest of Ecuador[END_REF][START_REF] Sanches | Unravelling remote sensing signatures of plants contaminated with gasoline and diesel: An approach using the red edge spectral feature[END_REF][START_REF] Nakata | Growth and physiological responses of Triticum aestivum and Deschampsia caespitosa exposed to petroleum coke[END_REF]. Oil contamination reduces water and nutrient availability for plants as well as root uptake capacities [START_REF] H. Ur R. Athar | Influence of sub-lethal crude oil concentration on growth, water relations and photosynthetic capacity of maize (Zea mays L.) plants[END_REF][START_REF] Balasubramaniyam | Scanning electron microscopic investigations of root structural modifications arising from growth in crude oil-contaminated sand[END_REF][START_REF] Nie | Plants' use of different nitrogen forms in response to crude oil contamination[END_REF]. This results in anatomical and biochemical alterations in leaves, including a decrease in pigment contents [START_REF] Balliana | Development of Canavalia ensiformis in soil contaminated with diesel oil[END_REF][START_REF] Baruah | Effect of crude oil contamination on the chlorophyll content and morpho-anatomy of Cyperus brevifolius (Rottb.) Hassk[END_REF]. On the brownfield, chlorophyll content significantly varied among species (p < 0.05). The chlorophyll contents of R. fruticosus and B.

davidii systematically attained the highest values, whereas those of Q. pubescens and P.

canadensis never exceeded 40 µg.cm -2 on average. Chlorophyll content also varied significantly for each species according to sampling dates (p < 0.05), as suggested by the temporal evolution of spectral signatures (Fig. 2). These variations were however less marked than those of other temperate species [START_REF] Demarez | Seasonal variation of leaf chlorophyll content of a temperate forest. Inversion of the PROSPECT model[END_REF]. Some smaller differences among species were observed on the control site, indicating natural differences in pigment contents [START_REF] Li | Factors influencing leaf chlorophyll content in natural forests at the biome scale[END_REF][START_REF] Martin | Genetic variation in leaf pigment, optical and photosynthetic function among diverse phenotypes of Metrosideros polymorpha grown in a common garden[END_REF][START_REF] Sánchez-Azofeifa | Differences in leaf traits, leaf internal structure, and spectral reflectance between two communities of lianas and trees: Implications for remote sensing in tropical environments[END_REF]. Pigment contents on the brownfield were thus expressed as percentages of those of the control to factor in these differences and to better compare species' sensitivities to TPH exposure (Fig. 4).

The chlorophyll content of C. sanguinea leaves taken from the brownfield remained close to that of the control site on all dates, indicating this species' high tolerance to oil contamination (Fig. 4).

Q. pubescens was the most sensitive species, regardless of the season. The leaf chlorophyll content almost remained below 60% of that of the control for this species. R. fruticosus and P. canadensis were more tolerant to TPH but their pigment contents still exhibited severe alterations (~70% in June and October). This is consistent with the findings of previous studies [START_REF] Lassalle | Hyperspectral signature analysis of three plant species to long-term hydrocarbon and heavy metal exposure[END_REF][START_REF] Lassalle | Detection and discrimination of various oil-contaminated soils using vegetation reflectance[END_REF]. Such differences in sensitivity among species have already been reported in tropical regions [START_REF] Arellano | Plant family-specific impacts of petroleum pollution on biodiversity and leaf chlorophyll content in the Amazon rainforest of Ecuador[END_REF], but the factors influencing this sensitivity remain poorly understood. Oil contamination, as found on our study site, is made of a complex cocktail of contaminants -mainly TPH and HM. All these contaminants are likely to affect leaf biochemistry, depending on each species' sensitivity to each of them. R. fruticosus, for example, is particularly tolerant to HM [START_REF] Dorrington | Some aspects of tissue accumulation and tolerance to available heavy metal ions by Rubus Fruticosus L. A colonizer of spoil tips in S.W. England[END_REF][START_REF] Nujkić | Impact of metallurgical activities on the content of trace elements in the spatial soil and plant parts of Rubus fruticosus L[END_REF]. Several species, including P. canadensis, developed detoxification mechanisms to protect themselves against the effects of HM and PAH [START_REF] Benyó | Physiological and molecular responses to heavy metal stresses suggest different detoxification mechanism of Populus deltoides and P. x canadensis[END_REF][START_REF] Kvesitadze | Biochemical mechanisms of detoxification in higher plants[END_REF][START_REF] Shahid | Foliar heavy metal uptake, toxicity comparison and detoxification in plants: A comparison of foliar and root metal uptake[END_REF].

Q. pubescens, P. canadensis and R. fruticosus exhibited a decrease in relative chlorophyll content in April (spring) (p < 0.05), after the appearance of new leaves (Fig. 4). Their sensitivity to TPH exposure was therefore more pronounced during this period, but remained unchanged between June and October in 2017 and 2018, despite important differences in environmental Leaf carotenoid content also varied among species, but did not always vary significantly from one season to the next. These results are shown in the Supporting Information. According to Archetti et al. [START_REF] Archetti | Unravelling the evolution of autumn colours: an interdisciplinary approach[END_REF], seasonal variations in leaf color -and thus in leaf reflectance -tend to result from changes in chlorophyll content rather than from carotenoid de novo synthesis. Carotenoid content was up to 60% lower on the brownfield than on the control site, indicating alterations induced by exposure to TPH. The observations made regarding relative chlorophyll content also apply to carotenoids. Previous studies mentioned that carotenoid pigments (i.e. carotenes and xanthophylls) are not affected in the same way by oil contaminants [START_REF] Lassalle | Hyperspectral signature analysis of three plant species to long-term hydrocarbon and heavy metal exposure[END_REF][START_REF] Lassalle | Detection and discrimination of various oil-contaminated soils using vegetation reflectance[END_REF]. For example, the lutein content might increase when plants are exposed to TPH, but the zeaxanthin content might be altered. Consequently, separating carotenoid pigments would benefit the detection of TPHaffected vegetation, but remains difficult when using hyperspectral spectroscopy [START_REF] Féret | PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments[END_REF][START_REF] Zhang | An extended PROSPECT: Advance in the leaf optical properties model separating total chlorophylls into chlorophyll a and b[END_REF].

The degree of pigment alteration therefore reflected species' sensitivity to TPH, as suggested by the spectral signatures observed on the brownfield (Fig. 2). Results indicate that this sensitivity varies by season, which is likely to affect TPH predictions, as described in section 3.4.

Prediction of TPH from leaf pigments

Leaf pigment contents from the five species were tested to estimate TPH concentrations in the soil. The results are summarized in Table 3. A strong relationship was observed on the train sets between leaf chlorophyll content and TPH concentrations for R. fruticosus, Q. pubescens and P.

canadensis in June and October (R 2 train ≥ 0.74). The more contaminated the soil, the lower the chlorophyll content (Fig. 5). This type of plant response has frequently been reported under experimental conditions [START_REF] H. Ur R. Athar | Influence of sub-lethal crude oil concentration on growth, water relations and photosynthetic capacity of maize (Zea mays L.) plants[END_REF][START_REF] Baruah | Effect of crude oil contamination on the chlorophyll content and morpho-anatomy of Cyperus brevifolius (Rottb.) Hassk[END_REF], but never in a natural context and using hyperspectral spectroscopy. Leaf chlorophyll content helped predict TPH concentrations with good (RPD > 2)

or excellent (RPD > 3) accuracy, as indicated by the results obtained on the test sets (Tab. 3).

Conversely, poor predictions were observed for B. davidii, confirming the absence of effects induced by TPH on this species. Over the entire study, the best predictions were obtained with R.

fruticosus and Q. pubescens, for which the lowest RMSE -3.63 and 3.20 g.kg -1 respectively -and the highest RPD -4.10 and 5.17 respectively -were observed (Fig. 5). Successful predictions were also achieved with P. canadensis, sometimes even exceeding those of the other species, such as in October 2017 (RMSE = 5.86 g.kg -1 and RPD = 3.04). Inaccuracies in pigment estimation also affected TPH predictions. For example, the RMSE observed for chlorophylls on R. fruticosus (2.62 µg.cm -2 , Fig. 3) was greater than the pigment variations for TPH concentrations between 40 and 60 g.kg -1 (Fig. 5). In that sense, making improvements to PROSPECT, including the separation of chlorophylls a and b, might help improve TPH predictions [START_REF] Zhang | An extended PROSPECT: Advance in the leaf optical properties model separating total chlorophylls into chlorophyll a and b[END_REF].

Predictions were more accurate for TPH concentrations of less than 30 g.kg -1 , because small variations rapidly induced changes in chlorophyll content (Fig. 5). Above 30 g.kg -1 , TPH were more difficult to predict as the chlorophyll content tended to stabilize. This type of difficulty might also arise outside of the range studied. Lethal concentrations of TPH would become impossible to predict for example. Likewise, predictions might be difficult with R. fruticosus, Q. pubescens and P. canadensis at very low concentrations (µg to mg.kg -1 ), because little or no pigment alterations are expected for these species. At such low concentrations, the use of a species more sensitive to TPH might however make predictions possible [START_REF] Zhu | Soil TPH concentration estimation using vegetation indices in an oil polluted area of eastern China[END_REF]. There was no observable relationship between leaf chlorophyll content and TPH concentrations in March and April for R. fruticosus, Q. pubescens and P. canadensis (Tab. 3). Consequently, the predictions of TPH on these dates were poor. The leaf chlorophyll content of R. fruticosus was the same as that of the control site in March (Fig. 4), meaning that the species was thus not affected by TPH, which explains the inaccurate predictions. April was characterized by the appearance of new leaves in the three species, along with the most pronounced alteration in leaf chlorophyll content. The plants' sensitivity to TPH exposure is particularly pronounced in the growing stages [START_REF] Lassalle | Assessing Soil Contamination Due to Oil and Gas Production Using Vegetation Hyperspectral Reflectance[END_REF][START_REF] Merkl | Phytoremediation in the tropics-the effect of crude oil on the growth of tropical plants[END_REF]. As a result, chlorophyll content in April 2017 displayed similar alterations regardless of the contamination level, which had an impact on TPH predictions. A similar response was observed for C. sanguinea leaves over the entire study (Tab. 3). In this case, chlorophyll content made it possible to detect but not to quantify TPH (RMSE > 15 g.kg -1 , RPD < 1.10). Consequently, the relationship between the leaf chlorophyll content and TPH concentrations observed for three out of the five species appeared after the growing period (April). This temporal pattern is of critical importance for the long-term monitoring of oil activities, as it shows that a single campaign carried out at the right time can be sufficient to assess the level of soil contamination. Leaf chlorophyll and carotenoid contents are usually well-correlated [START_REF] Jiang | Estimation of leaf traits from reflectance measurements: Comparison between methods based on vegetation indices and several versions of the PROSPECT model[END_REF][START_REF] Féret | PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments[END_REF]. However, the exposure of plants to TPH altered this correlation. The relationship between carotenoid content and TPH was weaker than for chlorophyll content (R 2 ≤ 0.60) and made predictions more difficult.

The results obtained with R. fruticosus, Q. pubescens and P. canadensis were very similar. The best RMSE was observed with R. fruticosus (8.31 g.kg -1 ), but represented poor predictions nonetheless (RPD < 2) (see Supporting Information). As discussed above, the contrasting effects of TPH on carotenoid pigments might have altered the relationship between leaf carotenoid content and TPH concentrations and therefore contributed to inaccuracies. Predictions might be greatly improved by separating carotenoid pigments and focusing on particular ones.

When compared to the proposed approach, vegetation indices and PLSR did not give any more accurate TPH predictions for B. davidii and C. sanguinea, confirming that these species are not suitable for predicting TPH. These results are presented in the Supporting Information. Predictions were better than those obtained with carotenoid content in the other species, but never reached those achieved with chlorophyll content. RMSE never dropped below 5.86 g.kg -1 and RPD never exceeded 2.57. Consequently, our approach outperformed vegetation indices and PLSR. One possible explanation is that in the VIS, leaf optical properties are mainly influenced by chlorophylls, and to a lesser extent by carotenoids -and other pigments -that were less correlated to TPH concentrations in our study [START_REF] Blackburn | Quantifying chlorophylls and carotenoids at leaf and canopy scales, Remote Sens[END_REF][START_REF] Sims | Relationships between leaf pigment content and spectral reflectance across a wide range of species, leaf structures and developmental stages[END_REF]. Therefore, TPH predictions using vegetation indices and PLSR suffered from the influence of carotenoids, which introduced inaccuracies. Conversely, the proposed approach using PROSPECT made it possible to isolate chlorophyll content, thereby improving TPH predictions.

Conclusion

This study demonstrates that TPH concentrations in oil-contaminated soils can be estimated from leaf optical properties using PROSPECT. The model performed well for simulating the spectral signatures and for retrieving pigment contents of leaves from exposed vegetation. Leaf chlorophyll content was strongly linked to TPH concentrations and allowed accurate predictions of these contaminants. The response of carotenoid pigments was complex and made predictions more difficult. Our approach was more effective than vegetation indices and PLSR, thanks to pigment separation. This study reveals that the accuracy of predictions strongly varies among species, because of their different sensitivities to TPH exposure. Three out of the five species studied (R. fruticosus, Q. pubescens and P. canadensis) were successful in predicting TPH.

Moreover, the sensitivity of a single species also varies throughout the year, according to season changes. Changes in the species' sensitivity, during the growing period (April) for example, affect TPH predictions.

Considering these factors, the proposed approach can be used as a reliable tool for assessing soil contamination on oil production facilities. In the future, the approach could be applied to airborne images with high spatial and spectral resolutions for monitoring large industrial areas with dense vegetation. For this purpose, the combination of the PROSPECT and SAIL models (called PROSAIL [START_REF] Jacquemoud | PROSPECT + SAIL models: A review of use for vegetation characterization[END_REF]), which takes into account canopy architecture, Leaf Area Index (LAI) and viewing geometry, will be addressed in further research.

Figure 1 .

 1 Figure 1. Flowchart of the method presented in this study (N: structure parameter, Cab:

Figure 2 .

 2 Figure 2. Seasonal evolutions of the mean spectral signatures for the five species on the

3. 2 . 1 .

 21 Simulated spectral signatures PROSPECT performed well for simulating the spectral signatures of leaves measured during the study. The DE algorithm always converged during inversions and the RMSE computed between the simulated and measured signatures ranged from 0.2 to 2% reflectance across the entire spectrum for the five species. RMSE, BIAS and SEPC values were consistent with those observed for other datasets[START_REF] Jiang | Estimation of leaf traits from reflectance measurements: Comparison between methods based on vegetation indices and several versions of the PROSPECT model[END_REF][START_REF] Féret | PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments[END_REF][START_REF] Liu | A novel principal component analysis method for the reconstruction of leaf reflectance spectra and retrieval of leaf biochemical contents[END_REF]. These results are presented in the Supporting Information, along with examples of simulated and measured signatures. Although spectral signatures varied according to the sampling date, the RMSE remained almost unchanged throughout the study.

  conditions among seasons and years. Only R. fruticosus and B. davidii were found on the brownfield in March 2017, and displayed chlorophyll content very similar to that of the control site.

Figure 4 .

 4 Figure 4. Temporal evolutions of the relative chlorophyll contents of leaves retrieved from

Figure 5 .

 5 Figure 5. Relationship between predicted leaf chlorophyll content and Total Petroleum

Table 1 .

 1 Description of the sampling campaigns carried out on the brownfield. Only R. fruticosus and C. sanguinea were found on all 23 plots, except in March 2018.

	Date	Species found (n)	Spectral signatures (n)	Season
	June 2017	5	309	Summer (heat waves)
	October 2017	5	309	Autumn
	March 2018	2	75	Late winter
	April 2018	5	309	Spring
	June 2018	5	309	Summer
	October 2018	5	309	Autumn

Table 2 .

 2 Parameter bounds used to perform PROSPECT inversions. Cab and Cxc bounds were extended from those of literature to factor in seasonal variations of leaf pigment contents.

	Parameter	Unit	Range
	Leaf structure (N)	-	1 -5
	Chlorophylls (Cab)	µg.cm -2	1 -100
	Carotenoids (Cxc)	µg.cm -2	1 -50
	Brown pigments (Cbp)	µg.cm -2	0.01-1
	Water (Cw)	g.cm -2	0.001 -0.1
	Dry matter (Cm)	g.cm -2	0.001 -0.1

Table 3 .

 3 Assessment of Total Petroleum Hydrocarbons (TPH) in the soil predicted from leaf chlorophyll content on the train and test datasets. RMSE, BIAS and SEPC are given in g.kg -1 (abs: absent from the site). The best prediction results are marked in bold.

	Species	Season		train			test
			R 2	RMSE	R 2	RMSE BIAS SEPC	RPD
		June 2017	0.06	16.78	0.04	16.24	3.04	15.95	1.00
		October 2017	0.05	19.27	0.07	18.14	-4.40	17.60	0.90
	B. davidii	March 2018 April 2018	0.01 0.01	17.01 21.13	0.03 0.02	16.47 -10.35 12.81 20.01 -5.50 19.24	0.99 0.81
		June 2018	0.04	16.18	0.07	15.98	0.12	15.98	1.02
		October 2018	0.01	19.04	0.02	17.76	-1.31	17.71	0.92
		June 2017	0.89	5.54	0.94	3.63	-0.21	3.62	4.10
		October 2017	0.61	10.44	0.74	7.64	1.56	7.48	1.95
	R. fruticosus	March 2018 April 2018	0.02 0.05	15.47 13.31	0.19 0.12	17.05 18.31	-0.93 -4.68	17.02 17.70	0.87 0.81
		June 2018	0.75	8.56	0.78	6.64	-0.58	6.62	2.24
		October 2018	0.79	6.50	0.77	8.28	-0.35	8.27	1.80
		June 2017	0.01	15.71	0.01	15.75	1.93	15.63	1.02
		October 2017	0.03	17.81	0.04	17.20	0.32	17.19	0.93
	C. sanguinea	March 2018 April 2018	abs. 0.01	abs. 14.24	abs. 0.02	abs. 17.33	abs. -1.34	abs. 17.28	abs. 0.92
		June 2018	0.01	15.88	0.01	15.98	-0.03	15.98	1.00
		October 2018	0.04	18.39	0.01	17.45	0.86	17.43	0.92
		June 2017	0.78	8.41	0.85	5.70	2.24	5.24	2.89
		October 2017	0.66	10.09	0.86	5.82	0.65	5.79	2.83
	Q. pubescens	March 2018 April 2018	abs. 0.01	abs. 15.80	abs. 0.03	abs. 17.94	abs. 0.84	abs. 17.92	abs. 0.92
		June 2018	0.72	8.93	0.96	3.20	-0.20	3.19	5.17
		October 2018	0.57	11.82	0.83	5.88	0.66	5.85	2.80
		June 2017	0.61	11.03	0.80	6.85	0.22	6.85	2.60
		October 2017	0.76	7.62	0.89	5.86	0.67	5.82	3.04
	P. canadensis	March 2018 April 2018	abs. 0.05	abs. 16.91	abs. 0.13	abs. 16.26	abs. 5.64	abs. 15.25	abs. 1.09
		June 2018	0.62	10.70	0.75	7.83	-1.16	7.74	2.27
		October 2018	0.74	8.93	0.81	6.62	-3.18	5.81	2.69

Acknowledgements

This work was performed in the frame of the NAOMI project between TOTAL and the ONERA, with the support of the Ecolab and Dynafor research units of Toulouse. The authors gratefully acknowledge E. Buffan-Dubau, D. Lambrigot and O. Berseille for their assistance in pigment analysis, and E. Hannan for correcting the manuscript. Declarations of interest: None. Funding: Financial support of this work was provided by TOTAL.


Yacine Bengana

Laurette S Tuckerman

Spirals and ribbons in counter-rotating Taylor-Couette flow: frequencies from mean flows and heteroclinic orbits

Keywords: numbers: 47, 20, Ky, 47, 20, Qr

A number of time-periodic flows have been found to have a property called RZIF: when a linear stability analysis is carried out about the temporal mean (rather than the usual steady state), an eigenvalue is obtained whose Real part is Zero and whose Imaginary part is the nonlinear Frequency. For two-dimensional thermosolutal convection, a Hopf bifurcation leads to traveling waves which satisfy the RZIF property and standing waves which do not. We have investigated this property numerically for counter-rotating Couette-Taylor flow, in which a Hopf bifurcation gives rise to branches of upwards and downwards traveling spirals and ribbons which are an equal superposition of the two. In the regime that we have studied, we find that both spirals and ribbons satisfy the RZIF property. As the outer Reynolds number is increased, the ribbon branch is succeeded by two types of heteroclinic orbits, both of which connect saddle states containing two axially stacked pairs of axisymmetric vortices. One heteroclinic orbit is non-axisymmetric, with excursions that resemble the ribbons, while the other remains axisymmetric.

I. INTRODUCTION

A number of time-periodic flows have been found to have the following property: when a linear stability analysis is carried out about the temporal mean (rather than the usual steady state from which they bifurcate), an eigenvalue emerges whose imaginary part reproduces the frequency of the nonlinear oscillation and whose real part is near zero. This property was first discovered by Barkley [START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF] and extensively studied [START_REF] Pier | On the frequency selection of finite-amplitude vortex shedding in the cylinder wake[END_REF][START_REF] Noack | A hierarchy of low-dimensional models for the transient and post-transient cylinder wake[END_REF][START_REF]Global linear stability analysis of time-averaged flows[END_REF][START_REF] Mantič-Lugo | Self-consistent mean flow description of the nonlinear saturation of the vortex shedding in the cylinder wake[END_REF][START_REF] Mantič-Lugo | A self-consistent model for the saturation dynamics of the vortex shedding around the mean flow in the unstable cylinder wake[END_REF][START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF] for the cylinder wake (for which it is strongly verified) and then for an open cavity [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF][START_REF] Meliga | Harmonics generation and the mechanics of saturation in flow over an open cavity: a second-order self-consistent description[END_REF][START_REF] Bengana | Bifurcation analysis and frequency prediction for open cavity flow[END_REF] (for which it holds quite well, though not exactly). In these two-dimensional flows, a horizontal velocity imposed at infinity is deviated as it encounters objects (a cylinder) or openings (a cavity) and, for sufficiently high Reynolds number, undergoes a Hopf bifurcation.

In Turton, Tuckerman & Barkley [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF], this property was named RZIF -a mnemonic for the fact that the eigenvalue obtained by linearizing about the mean flow has a Real part of Zero and an Imaginary part which is the nonlinear Frequency -and was studied for a quite different kind of flow. Thermosolutal convection is driven by a vertical density gradient in a fluid mixture between two horizontal plates, which in turn is caused by imposed temperature and concentration differences at the two bounding plates. When the temperature and concentration have opposing effects on the density, instability is manifested as a Hopf bifurcation. In a two-dimensional domain with horizontally periodic boundary conditions, this is a classic case of the breaking of O(2) symmetry (translations and reflection) and leads to the simultaneous formation of traveling and standing wave branches, at most one of which can be stable [START_REF] Knobloch | Oscillatory convection in binary mixtures[END_REF]. Turton et al. [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF] discovered that the traveling waves in thermosolutal convection provide another clear realisation of the RZIF property, while the standing waves provide an equally clear counterexample: the real part of the eigenvalue obtained by linearizing about the mean part of the standing waves is far from zero and its imaginary part is far from the nonlinear frequency.

We have sought to continue this investigation on another pair of hydrodynamic traveling and standing waves. In Taylor-Couette flow, when only the inner cylinder rotates, a steady bifurcation leads to the classic Taylor vortices. In counter-rotating Taylor-Couette flow, however, these vortices are inclined to the rotation axis, and hence nonaxisymmetric. A Hopf bifurcation then gives rise to branches of the well-known spirals [START_REF] Krueger | On the relative importance of Taylor-vortex and non-axisymmetric modes in flow between rotating cylinders[END_REF][START_REF] Langford | Primary instabilities and bicriticality in flow between counter-rotating cylinders[END_REF][START_REF] Antonijoan | Non-linear spirals in the Taylor-Couette problem[END_REF] and the somewhat lesser-known ribbons [START_REF] Demay | Calcul des solutions bifurquées pour le problème de Couette-Taylor avec les 2 cylindres en rotation[END_REF][START_REF] Tagg | Nonlinear standing waves in Couette-Taylor flow[END_REF][START_REF] Chossat | The Couette-Taylor Problem[END_REF][START_REF] Pinter | Competition between traveling fluid waves of left and right spiral vortices and their different amplitude combinations[END_REF][START_REF] Pinter | Bifurcation of standing waves into a pair of oppositely traveling waves with oscillating amplitudes caused by a three-mode interaction[END_REF][START_REF] Pinter | Wave-number dependence of the transitions between traveling and standing vortex waves and their mixed states in the Taylor-Couette system[END_REF][START_REF] Hoffmann | Transitions between Taylor vortices and spirals via wavy Taylor vortices and wavy spirals[END_REF][START_REF] Altmeyer | Secondary bifurcation of mixed-cross-spirals connecting travelling wave solutions[END_REF][START_REF] Deguchi | Fully nonlinear mode competitions of nearly bicritical spiral or Taylor vortices in Taylor-Couette flow[END_REF], which, like standing waves, are an equal superposition of spirals moving axially upwards and downwards. As in [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF][START_REF] Knobloch | Oscillatory convection in binary mixtures[END_REF], this is a consequence of the breaking of O(2) symmetry, here in the axial direction, when end effects are disregarded. Figure 1 shows the appearance of these two types of flows.

Counter-rotating Taylor-Couette flow has an additional symmetry when compared to the two-dimensional thermosolutal cases studied in [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF][START_REF] Knobloch | Oscillatory convection in binary mixtures[END_REF]. Along with the O(2) axial z symmetry, the Taylor-Couette configuration has SO(2) symmetry, consisting of rotations (not reflection) in the azimuthal θ direction. The tilted spirals or azimuthally tiled ribbons are not axisymmetric, and hence they break SO(2) symmetry. In such cases, bifurcations always lead to states which rotate in the azimuthal direction; see e.g. [START_REF] Ecke | Hopf bifurcation with broken reflection symmetry in rotating Rayleigh-Bénard convection[END_REF]. Thus, although ribbons are standing waves in the axial direction, i.e. equal superpositions of axially-upwards and axially-downwards traveling waves, they rotate in the azimuthal direction, as do the spirals.

As the rotation velocity of the outer cylinder is varied, the ribbon branch is succeeded by two types of heteroclinic orbits, both of which connect saddle states consisting of two pairs of axisymmetric vortices. One heteroclinic orbit is non-axisymmetric, with excursions that resemble the ribbons, while the other remains axisymmetric.

II. METHODS

To investigate this problem numerically, we have used the spectral finite-difference code of Willis [START_REF] Dessup | The self-sustaining process in Taylor-Couette flow[END_REF] to solve the Navier-Stokes equations in a cylindrical annulus,

∂ t U = -(U • ∇) U -∇P + ∇ 2 U (1a) ∇ • U = 0 (1b) U = e θ Re in , e θ Re out at r in , r out (1c)
where the length scale is d ≡ r out -r in and the time scale is d 2 /ν, where ν is the kinematic viscosity. The radius ratio is η ≡ r in /r out , the two Reynolds numbers are defined as Re j ≡ r j Ω j d/ν where j = in, out, and the numerical resolution is (N r , N θ , N z ) = (33, [START_REF] Tagg | Nonlinear standing waves in Couette-Taylor flow[END_REF][START_REF] Tagg | Nonlinear standing waves in Couette-Taylor flow[END_REF].

We use the parameters in Pinter et al. [START_REF] Pinter | Competition between traveling fluid waves of left and right spiral vortices and their different amplitude combinations[END_REF][START_REF] Pinter | Bifurcation of standing waves into a pair of oppositely traveling waves with oscillating amplitudes caused by a three-mode interaction[END_REF][START_REF] Pinter | Wave-number dependence of the transitions between traveling and standing vortex waves and their mixed states in the Taylor-Couette system[END_REF], namely radius ratio η ≡ r in /r out = 0.5, inner Reynolds number Re in = 240 and azimuthal wavenumber M 0 = 2 (azimuthal wavelength 2π/M 0 = π). We have chosen the axial wavelength λ z = 1.2, slightly less than the value studied intensively by [START_REF] Pinter | Competition between traveling fluid waves of left and right spiral vortices and their different amplitude combinations[END_REF][START_REF] Pinter | Bifurcation of standing waves into a pair of oppositely traveling waves with oscillating amplitudes caused by a three-mode interaction[END_REF][START_REF] Pinter | Wave-number dependence of the transitions between traveling and standing vortex waves and their mixed states in the Taylor-Couette system[END_REF] to avoid a secondary bifurcation they describe near the onset of spirals and ribbons for λ z = 1.3. Because counter-rotating Taylor-Couette flow is centrifugally unstable only near the inner cylinder, here for r ∈ [1, 1.24], the preferred axial wavelength is correspondingly less than that which would be expected if the vortices filled the gap. Because Re in remains fixed and we will vary Re out , we will simplify the notation by using Re instead of Re out to designate the outer Reynolds number.

We have modified the Taylor-Couette code [START_REF] Dessup | The self-sustaining process in Taylor-Couette flow[END_REF] in order to timestep the equations linearized about a velocity field U:

∂ t u = -(U • ∇) u -(u • ∇) U -∇p + ∇ 2 u (2a) ∇ • u = 0 (2b) u = 0 at r in , r out (2c)
Time-integration of the linearized equations (2) about a steady solution U converges to the leading eigenvector, i.e. that whose corresponding eigenvalue has the most positive or least negative real part or growth rate. It is a means of carrying out the power method on the exponential of the right-hand-side of (2). Integration of (2) does not generally lead to a steady or periodic state, since the nonlinear terms which would otherwise saturate the amplitude is absent. Instead, it leads to states whose amplitude increases or decreases, but with a constant growth rate and frequency, in the case of a complex eigenvalue. We first use the linearized code to compute the Hopf bifurcation threshold from laminar Couette flow. For various timesteps ∆t, we calculate the growth rate σ as a function of Re and then use linear extrapolation or interpolation to find Re c such that σ(Re c ) = 0. (This procedure circumvents critical slowing down and so finds bifurcation thresholds faster and more accurately than could be done with nonlinear simulations.) The bifurcation thresholds calculated using ∆t = 10 -4 , 10 -5 , and 10 The branches of spirals and ribbons illustrated in figure 1 bifurcate at Re ≈ -585 towards increasing Re. Their functional form is:

u SPI± (r, θ, z, t) = m ûm (r) e im(z/λz±M0θ-ωSPIt) u RIB (r, θ, z, t) = k,m    ûkm,r (r)e r sin(kz/λ z ) ûkm,θ (r)e θ sin(kz/λ z ) ûkm,z (r)e z cos(kz/λ z )    e im(M0θ-ωRIBt) (3)
The spirals are traveling waves: they depend on z, θ, and t only via the single combination z/λ z + M 0 θ -ω SPI t. The ribbons have fixed axial nodal lines and are reflection symmetric; the axis of symmetry is z = λ z /2 in equation ( 3) and in figure 1. The ribbons are traveling waves in θ: they depend on θ and t via M 0 θ -ω RIB t. Many other complex and esthetic symmetry-breaking flows have been observed in numerical simulations of this configuration [START_REF] Pinter | Competition between traveling fluid waves of left and right spiral vortices and their different amplitude combinations[END_REF][START_REF] Pinter | Bifurcation of standing waves into a pair of oppositely traveling waves with oscillating amplitudes caused by a three-mode interaction[END_REF][START_REF] Pinter | Wave-number dependence of the transitions between traveling and standing vortex waves and their mixed states in the Taylor-Couette system[END_REF][START_REF] Hoffmann | Transitions between Taylor vortices and spirals via wavy Taylor vortices and wavy spirals[END_REF][START_REF] Altmeyer | Secondary bifurcation of mixed-cross-spirals connecting travelling wave solutions[END_REF][START_REF] Deguchi | Fully nonlinear mode competitions of nearly bicritical spiral or Taylor vortices in Taylor-Couette flow[END_REF], including cross spirals, oscillating cross spirals, mixed cross spirals and mixed ribbons. We will use numerical simulations of the nonlinear and linearized equations ( 1) and ( 2) to investigate the RZIF phenomenon mentioned in the introduction, in which linearization about the mean flow of a periodic solution yields an eigenvalue whose Real part is Zero and whose Imaginary part is the nonlinear Frequency. To do so, we will calculate the spirals and ribbons by time integrating the nonlinear equations (1) until an asymptotic state is reached. Because the spirals and ribbons are both azimuthally traveling waves, their temporal means are identical with their azimuthal means which are, in turn, the axisymmetric (m = 0) components of (3). We will then find the eigenvalues and eigenvectors about the resulting mean flows by integrating equations (2).

III. RZIF ANALYSIS

We begin by showing the bifurcation diagram for spirals and ribbons. Their common threshold is Re ≈ -585. On figure 2a, we represent the square root of kinetic energy of the two branches. In this parameter range, the spirals are unstable to ribbons, but can be computed by including only k = m modes or, like [START_REF] Pinter | Competition between traveling fluid waves of left and right spiral vortices and their different amplitude combinations[END_REF], by setting the m = 0 component of u r to zero. The ribbons themselves become unstable for Re -530. Figure 2b shows the frequency of the nonlinear spirals and ribbons, as well as those obtained by linearizing about laminar Couette flow (dashed line). Linearization about the laminar flow leads to the same set of eigenvalues for spirals or ribbons, since the eigenvectors corresponding to upgoing and downgoing spirals are related by symmetry, with their superposition leading to ribbons. For a supercritical Hopf bifurcation, these necessarily match the frequency of the nonlinear branches at the threshold. As Re is increased, the frequencies of the nonlinear spirals deviate greatly from those derived from the laminar flow, but the frequencies of the nonlinear ribbons remain close. This was also true of the standing waves of thermosolutal The dark green circles (upper set) are obtained by linearizing around a mean spiral flow whose orientation is opposite to that of the eigenvector while the blue ones (lower set) are obtained by linearizing around the mean flow of a ribbon. convection [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF], which inspired this study: the frequencies of the nonlinear standing waves were much closer to those obtained by linearization about the laminar (conductive) flow than were the frequencies of the nonlinear traveling waves.

Figure 3 shows the temporal means of the nonlinear spiral and ribbon states shown in figure 1 at Re = -550. Since spirals are traveling waves in both the θ and the z direction, their temporal means are equivalent to their spatial means and have no dependence on either direction. Ribbons are traveling waves in the θ direction and so their means have no θ dependence. We set U in the linearized equations (2) to the mean flows of spirals or ribbons shown in figure 3 and initialize u with fields of the spiral or ribbon form, respectively. We then integrate until the growth rate and the frequency are constant. The resulting field u(t) then cycles between the real and imaginary parts of the eigenvector which here are related by a spatial shift. Figure 4 shows samples of these eigenvectors. These greatly resemble the deviations of the corresponding nonlinear flows from laminar Couette flow shown in figure 1 and resemble even more the usual eigenvectors of spiral and ribbon form resulting from linearization about laminar Couette flow. Figure 5 presents the main result of this section, namely that both the spirals and the ribbons satisfy the RZIF property. The growth rates or real parts of eigenvalues are shown in figure 5a and the frequencies or imaginary parts are shown in figure 5b. The solid curves show the frequencies obtained from nonlinear simulations. Linearization about the mean flows of spirals (red crosses) or of ribbons (blue stars) lead to eigenvalues whose imaginary parts agree with the nonlinear frequencies. The corresponding growth rates are seen to be near zero, i.e. the mean flows are nearly marginally stable. (This has no relation to the marginal stability of the full spirals or ribbons to shifted versions of themselves, which arises from their θ or z dependence in a homogeneous domain.) Thus, both the spirals and the ribbons satisfy the RZIF property. This contrasts with the thermosolutal case, in which the traveling waves satisfied the RZIF property, while the standing waves did not.

Additional eigenvalues are shown in figure 6, which are not related to the basic RZIF property. Both correspond to spiral eigenvectors and have positive growth rates. These show the frequency and growth rate of spiral eigenvectors superposed on the mean flow of spirals of opposite helicity (dark green circles) or of ribbons (blue circle). The positive growth rate for the spirals of opposite chirality could be related to the linear instability of the spirals to ribbons; such a transition would take place via the addition of spirals of opposite chirality until equal amplitudes are reached. However, for these parameter values, ribbons show no tendency to make a transition to spirals, and yet the growth rates from the mean flow of ribbons to spirals are also positive (although less so than the growth rates from the mean flow of spirals). More understanding of the meaning of linearization about mean flows would be necessary to interpret these eigenvalues.

The RZIF property is implied by a near-monochromatic spectrum, as shown in [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF] as follows. Consider the evolution equation

∂ t U = LU + N (U, U) (4)
where L is linear and

N (•, •) is a quadratic nonlinearity. Let U = U + n =0 u n e inωt (5)
(with u -n = u * n ) be the temporal Fourier decomposition of a periodic solution to (4) with mean U and frequency ω. The n = 1 component of ( 4) is

iωu 1 = Lu 1 + N (U, u 1 ) + N (u 1 , U) L Ū u 1 + N (u 2 , u -1 ) + N (u -1 , u 2 ) + N (u 3 , u -2 ) + N (u -2 , u 3 ) + . . . N 1 (6)
If, as is often the case, ||u n || ∼ |n| , then N 1 = O( 3 ) may be neglected and RZIF is satisfied: the linear operator L Ū in (6) has the pure imaginary eigenvalue iω, corresponding to the frequency of the periodic solution. Hence the RZIF property is satisfied for near-monochromatic oscillations in any system with a quadratic nonlinearity. It remains to verify whether these periodic solutions are near-monochromatic. Figure 7 shows the temporal Fourier spectra of the azimuthal velocity for the spirals and the ribbons. For the thermosolutal case of [START_REF] Turton | Prediction of frequencies in thermosolutal convection from mean flows[END_REF], the ratio ||u 2 ||/||u 1 || was found to be much higher for the standing waves (greater than about 10 -1 ) than for the traveling wave branch (between 10 -3 and about 10 -1 ), explaining the satisfaction of RZIF for the traveling waves and not for the standing waves. In contrast, for this counter-rotating Taylor-Couette case, figure 7 shows that the temporal spectra for spirals and ribbons are quite similar. The ratio ||u 2 ||/||u 1 || is consistently less than 10 -2 for both flows over the range of our investigation. This is consistent with the fact that RZIF is satisfied for both flows.

IV. HETEROCLINIC ORBITS et al. [START_REF] Pinter | Competition between traveling fluid waves of left and right spiral vortices and their different amplitude combinations[END_REF][START_REF] Pinter | Bifurcation of standing waves into a pair of oppositely traveling waves with oscillating amplitudes caused by a three-mode interaction[END_REF][START_REF] Pinter | Wave-number dependence of the transitions between traveling and standing vortex waves and their mixed states in the Taylor-Couette system[END_REF] observed that the ribbon branch became unstable and was succeeded by oscillating cross-spirals as Re was increased. We observe this as well, near Re ≈ -525. As Re is further increased, however, the oscillating cross-spirals are themselves succeeded by a near-heteroclinic orbit, in which the system spends the overwhelming majority of its time in one of two saddles containing two pairs of axisymmetric vortices. The two saddles differ by a quarter of an axial period (one vortex). The system passes from one to the other via rapid non-axisymmetric excursions resembling the ribbon states. Snapshots of the phases of an excursion taken from a heteroclinic cycle at Re = -502 are shown in figure 8.

Figure 9(a,b,c) shows timeseries from this heteroclinic orbit. Time spent at the two saddles are seen as two plateaus. The approaches to and departures from the saddles are oscillatory; this appears clearly in the enlarged timeseries of figure 9b, as well as in the spirals of the phase portrait of figure 10a. In the phase portrait, the two saddles are distinguished, and the spiraling approaches to them are approximately perpendicular to one another. In the timeseries of figure 9c, the energy in the non-axisymmetric modes is seen to increase and to decrease logarithmically.

As Re is further increased, this regime is succeeded by another type of heteroclinic orbit, visualized in figure 11. Timeseries for this cycle are shown in figure 9(d,e,f) for Re = 490. The two plateaus are the same as for the previous heteroclinic cycle, but the approaches to and departures from the saddles are not oscillatory. These approaches and departures differ from excursion to excursion, as highlighted in the enlargements in figure 9e. The difference between the excursions is also seen in the phase portraits in figure 10b, which begins with a transient oscillatory phase. The black dots in the phase portrait are equally spaced in time and so accumulate on the saddles. Trajectories are seen approaching and leaving the saddles along different paths at different angles. Turning to the logarithmic energy timeseries of figure 9(f), after an initial transient (corresponding to the oscillatory transient in the phase portrait), the energy in the non-axisymmetric modes becomes and remains very small. Instead, it is the energy in the odd axial wavenumber modes which tracks the heteroclinic cycle. Since the saddles contain two stacked pairs of vortices, their (a,d) Complete timeseries alternating between two plateaus. (b,e) Enlargements of timeseries. b) Enlargement of a single period from a non-axisymmetric cycle, with further enlargement around the excursion. The differences between the growth/decay rates and frequencies during the spiralling approach and departure are clear. e) Enlargements near the peaks of the axisymmetric non-oscillatory excursions, highlighting their variation. (c,f) Logarithmic plot of energy, which grows and decays exponentially. The non-axisymmetric energy is shown in black and, for the axisymmetric heteroclinic cycle in (f), the energy corresponding to axisymmetric but axially antisymmetric modes is shown in red. Both types of heteroclinic cycles have the same plateaus, up to numerical roundoff.

u (3) θ - 10
FIG. 10. a) Phase portraits for a) oscillatory heteroclinic cycle at Re = -502 and for b) non-oscillatory heteroclinic cycle at Re = -490 using ur(r = 1.2, θ = 0, z = 0), ur(r = 1.2, θ = π/3, z = λz/3), and u θ (r = 1.2, θ = 2π/3, z = 2λz/3) corresponding to timeseries in figure 9. The trajectory in the oscillatory heteroclinic cycle spirals inwards along one curved surface and outwards along a perpendicular curved surface. Non-oscillatory phase portrait begins with a spiralling transient, then follows different paths into and out of the saddles. The color coding (blue and red) of these paths matches that in the timeseries, where can be seen subtle differences between the excursions. The saddles for the two cycles are the same; note the different scale for u axial Fourier spectrum contain only even multiples of the fundamental wavenumber 2π/λ z . Excursions are manifested by the logarithmic increase and decrease of the energy contained in odd multiples of 2π/λ z in the axial spectrum. The variation from cycle to cycle is also seen in the logarithmic timeseries.

The axisymmetric heteroclinic cycle is very similar to that observed in the 1:2 mode-interaction [START_REF] Armbruster | Heteroclinic cycles and modulated travelling waves in systems with O(2) symmetry[END_REF][START_REF] Nore | The 1:2 mode interaction in exactly counter-rotating von Kármán swirling flow[END_REF]. An erratic alternation between a few different types of excursions is also seen in the cycles of [START_REF] Nore | The 1:2 mode interaction in exactly counter-rotating von Kármán swirling flow[END_REF]. As explained in these references, this cycle is a consequence of the interaction of axial modes with wavenumbers 2π/λ z and 4π/λ z . Furthermore, the heteroclinic orbits can be characterized quantitatively by calculating the rates of approach and departure from the saddles which anchor the heteroclinic orbits. These are the two leading eigenvalues (in the axisymmetric subspace) of the saddles and are seen as the slopes of the black curve in the logarithmic timeseries in figure 9f. We plot these as a function of Re in figure 12b. The existence of the heteroclinic cycle requires eigenvalues of opposite signs; thus, the axisymmetric cycle exists for Re > -502, where one of the leading eigenvalues becomes positive while the other remains negative. Stability of the cycle requires that the rate of departure be no greater than the rate of approach, i.e. the absolute value of the negative eigenvalue must exceed the positive eigenvalue. This occurs for Re < -482. Hence the axisymmetric cycle exists and is stable within this axisymmetric subspace in the range Re ∈ [-502, -482].

The theory in [START_REF] Armbruster | Heteroclinic cycles and modulated travelling waves in systems with O(2) symmetry[END_REF][START_REF] Nore | The 1:2 mode interaction in exactly counter-rotating von Kármán swirling flow[END_REF] does not concern oscillatory heteroclinic cycles or complex eigenmodes. However, assuming that the analysis of the oscillatory cycles is similar to that of the non-oscillatory cycles, we plot the slopes (or real parts of the eigenvalues) seen in the logarithmic timeseries of figure 9c as a function of Re. We also plot the frequencies associated with these eigenvalues in figure 13a, associated with the oscillations in the timeseries in figure 9a,b and the spiraling in and out in the phase portraits of figure 10a. Using the same criteria as for the non-oscillatory cycles, figure 12a shows that the oscillatory heteroclinic cycle exists and is stable within the subspace of these eigendirections for Re ∈ [-522, -476].

The theory in [START_REF] Armbruster | Heteroclinic cycles and modulated travelling waves in systems with O(2) symmetry[END_REF][START_REF] Nore | The 1:2 mode interaction in exactly counter-rotating von Kármán swirling flow[END_REF] also does not address the relative stability of two kinds of heteroclinic cycles. The transition that we observe when Re is reduced from -502 to -490 implies that the non-axisymmetric heteroclinic cycles is unstable to the axisymmetric one at this value of Re. In fact, we observe non-axisymmetric heteroclinic cycles for Re < -496 and axisymmetric cycles for Re > -496, separated by the point at which axisymmetric and nonaxisymmetric departure rates cross, as seen in figure 13b. This supports the idea that the cycle whose departure rate from the saddles is smaller is unstable to the cycle whose departure rate is greater [START_REF] Ashwin | Attractors for robust heteroclinic cycles with continua of connections[END_REF]. Combining all of the criteria, the stability range of the non-axisymmetric cycle is Re ∈ [-522, -496], while that of the axisymmetric cycle is Re ∈ [-496, -482]. 12. Rate of approach to and departure from saddles as a function of Re. These are the slopes of the logarithmic plots of figure 9 and also the real parts of the leading eigenvalues. a) Non-axisymmetric oscillatory heteroclinic cycle. The cycle exists when the eigenvalues are of opposite sign, here for Re < -476. The cycle is stable within the space of these eigenmodes when the rate of approach exceeds the rate of departure, here for Re > -522. b) Axisymmetric non-ocillatory heteroclinic cycle. The cycle exists when the eigenvalues are of opposite sign, here for Re > -502. The cycle is stable within the space of these eigenmodes when the rate of approach exceeds the rate of departure, here for Re < -482. which matches our observations.

V. DISCUSSION

The problem of characterizing the nonlinear frequency of a periodic flow away from the threshold has led to one proposed solution: when a linear stability analysis is carried out about the temporal mean, an eigenvalue is obtained whose Real part is Zero and whose Imaginary part is the nonlinear Frequency. This RZIF property has heretofore been studied for the cylinder wake, an open cavity, and thermosolutal convection. It is strongly satisfied for the cylinder wake and the traveling waves of thermosolutal convection with oppositely directed density gradients, weakly satisfied for the open cavity, and not at all satisfied for the standing waves of thermosolutal convection. Although the RZIF property is a natural outcome of a near-monochromatic temporal spectrum, it is not entirely clear when this occurs.

We have investigated RZIF in another configuration, counter-rotating Taylor-Couette flow, which, like the thermosolutal convective case, has a Hopf bifurcation leading to branches of traveling and standing waves, here manifested as spirals and ribbons. In the thermosolutal case, the traveling waves display the RZIF property and have a temporal Fourier spectrum which is highly peaked, while the standing waves do not. Here, both the spirals and the ribbons display the RZIF property and both have temporal spectra which are equally peaked. However, the ribbons are standing waves only in the axial direction and are traveling waves in the azimuthal direction. This may be the cause of the difference between the Taylor-Couette and cases. The search for another counter-example to RZIF, especially in a purely hydrodynamic flow without additional fields, remains open.

During the course of our investigation, we discovered two heteroclinic cycles, one non-axisymmetric and the other axisymmetric. Both cycles are anchored by the same saddles: axisymmetric states containing two axially stacked pairs of vortices. The axisymmetric heteroclinic cycle is a manifestation of the classic 1:2 mode interaction, here the interaction between axisymmetric states with one and two pairs of vortices. The non-axisymmetric heteroclinic cycle is not of this type, since the approach to and from the saddles is oscillatory via transients resembling ribbons. However, it is anchored by the same axisymmetric saddles. The existence of two qualitatively different heteroclinic cycles connecting the same saddles is an unusual and intriguing feature. These cycles will be studied in more detail in a later investigation.

FIG. 1 .

 1 FIG. 1. Nonlinear spirals (a,b) and ribbons (c,d) at Re = -550. The deviation from laminar Couette flow is shown. Spirals of the opposite helicity exist and have identical properties but are not shown. a,c) (θ, z) velocity (arrows) with ur (colors) at r = 1.26, near the maximum of the deviation from laminar Couette flow. b,d) (r, z) velocity (arrows) with u θ (colors).

FIG. 2 .

 2 FIG. 2. a) Bifurcation diagram for spirals and ribbons. Square root of kinetic energy of the deviation from laminar Couette flow is shown as a function of Re. The ribbons are stable where shown, while the spirals are calculated by using a spiral initial condition. b) Frequencies of spiral and ribbon states as a function of Re. The dashed line shows the eigenvalues obtained by linearizing about laminar Couette flow.

FIG. 3 .FIG. 4 .FIG. 5 .FIG. 6 .

 3456 FIG. 3. Mean flow of spirals (a,b) and ribbons (c,d) at Re = -550. The deviation from laminar Couette flow is shown. a,c) (θ, z) velocity (arrows) with ur (colors) at r = 1.26, near the maximum of the deviation from laminar Couette flow. b,d) (r, z) velocity (arrows) with u θ (colors). Since the spirals are traveling waves in the θ and z directions, their means have no dependence on either. The ribbons are traveling waves in θ and hence their mean is axisymmetric.

FIG. 7 .

 7 FIG. 7. Spectra of spirals and ribbon states. a) Full spectrum for states at Re = -550, normalized by amplitude of fundamental frequency. b) Ratio of second harmonic to fundamental frequency.

FIG. 8 .

 8 FIG. 8. Non-axisymmetric heteroclinic cycle at Re = -502. Colors show the radial velocity in a (θ, z) slice at r = 1.26. The first (top left) and last (bottom right) visualizations show the saddles which anchor the heteroclinic cycles. These contain two pairs of axisymmetric vortices and differ by an axial phase shift of a single vortex. The visualizations between the two illustrate the excursion or rapid transition between the two saddles, via states which resemble ribbons.

FIG. 9 .

 9 FIG.9. Timeseries for (a,b,c) oscillatory and non-axisymmetric heteroclinic cycle at Re = -502 and for (d,e,f) non-oscillatory and axisymmetric heteroclinic cycle at Re = -490. (a,d) Complete timeseries alternating between two plateaus. (b,e) Enlargements of timeseries. b) Enlargement of a single period from a non-axisymmetric cycle, with further enlargement around the excursion. The differences between the growth/decay rates and frequencies during the spiralling approach and departure are clear. e) Enlargements near the peaks of the axisymmetric non-oscillatory excursions, highlighting their variation. (c,f) Logarithmic plot of energy, which grows and decays exponentially. The non-axisymmetric energy is shown in black and, for the axisymmetric heteroclinic cycle in (f), the energy corresponding to axisymmetric but axially antisymmetric modes is shown in red. Both types of heteroclinic cycles have the same plateaus, up to numerical roundoff.

FIG. 11 .

 11 FIG. 11. Axisymmetric heteroclinic cycle at Re = -490. Colors show the radial velocity in a (θ, z) slice at r = 1.26. The first (top left) and last (bottom right) visualizations show the saddles which anchor the heteroclinic cycles. These contain two pairs of axisymmetric vortices and differ by an axial phase shift of a single vortex. For this heteroclinic cycle, the states remain axisymmetric throughout the rapid transition between the two saddles.

  FIG.12. Rate of approach to and departure from saddles as a function of Re. These are the slopes of the logarithmic plots of figure9and also the real parts of the leading eigenvalues. a) Non-axisymmetric oscillatory heteroclinic cycle. The cycle exists when the eigenvalues are of opposite sign, here for Re < -476. The cycle is stable within the space of these eigenmodes when the rate of approach exceeds the rate of departure, here for Re > -522. b) Axisymmetric non-ocillatory heteroclinic cycle. The cycle exists when the eigenvalues are of opposite sign, here for Re > -502. The cycle is stable within the space of these eigenmodes when the rate of approach exceeds the rate of departure, here for Re < -482.

FIG. 13

 13 FIG.13. a) Imaginary parts of eigenvalues involved in non-axisymmetric oscillatory heteroclinic cycle. b) Real parts of eigenvalues involved in both heteroclinic cycles, axisymmetric and non-axisymmetric. Non-axisymmetric cycles are observed for Re < -496, when the largest eigenvalue (departure rate, dark green circles) corresponds to a non-axisymmetric eigenpair. Axisymmetric cycles are observed for Re > -496, when the largest eigenvalue (blue stars) corresponds to an axisymmetric eigenvector. The other endpoints of the stability ranges, shown as the black rightwards and red leftwards arrows, involve the next-to-largest eigenvalues (rates of approach, red diamonds and black squares); see figure12and its caption. Summarizing, the stability range of the non-axisymmetric cycle is Re ∈ [-522, -496], while that of the axisymmetric cycle is Re ∈ [-496, -482], which matches our observations.

ACKNOWLEDGMENTS

This work was performed in part using high performance computing resources provided by the Institut du Developpement et des Ressources en Informatique Scientifique (IDRIS) of the Centre National de la Recherche Scientifique (CNRS), coordinated by GENCI (Grand Equipement National de Calcul Intensif) through grant A0042A01119.