Christel Vidaller
email: christel.vidaller@gmail.com

Thierry Dutoit

Yosra Ibrahim

Hans Martin Hanslin

Armin Bischoff

Adaptive differentiation among populations of the Mediterranean dry grassland

Keywords: Common garden, local adaptation, plant origin, Poaceae, population differentiation, seedling recruitment

come

moisture often results in higher competition for other resources that represent a counteracting selective force [START_REF] Kadmon | Plant Competition along Soil Moisture Gradients: A Field Experiment with the Desert Annual Stipa Capensis[END_REF].

Ecosystem management is a third potential driver of adaptive differentiation in semi-natural systems such as grasslands. Grazing changes site conditions by biomass removal, plant damage, reduced competition, and creates heterogeneity in grassland vegetation [START_REF] Bullock | Restoration of biodiversity enhances agricultural production[END_REF][START_REF] Hufford | Local Adaptation and the Effects of Grazing on the Performance of Nassella pulchra: Implications for Seed Sourcing in Restoration[END_REF]. It seems to be likely that long-term grazing favors genotypes adapted to grazing [START_REF] Hufford | Local Adaptation and the Effects of Grazing on the Performance of Nassella pulchra: Implications for Seed Sourcing in Restoration[END_REF]. However, grazing may also obscure adaptation to other environmental factors if this selection pressure has not been consistent in the past [START_REF] Hufford | Local Adaptation and Effects of Grazing among Seedlings of Two Native California Bunchgrass Species: Implications for Restoration[END_REF].

Genetic differentiation is not only the result of adaptation to different environments. Stochastic processes such as genetic drift considerably contribute to population differentiation and may even overwhelm local adaptation [START_REF] Hereford | Limits to local adaptation in six populations of the annual plant Diodia teres[END_REF][START_REF] Hereford | A Quantitative Survey of Local Adaptation and Fitness Trade-Offs[END_REF]. Genetic drift and limited dispersal favor isolation by distance effects that are not related to environmental gradients [START_REF] Sexton | Genetic isolation by environment or distance: Which pattern of gene flow is most common?[END_REF]. [START_REF] Durka | Genetic differentiation within multiple common grassland plants supports seed transfer zones for ecological restoration[END_REF] found both isolation by distance caused by stochastic processes, and isolation by environment in seven tested plant species. The contribution of adaptive differentiation to genetic structure may even be low compared with genetic drift (Michalski andDurka, 2012, Sexton et al. 2013).

Genetic differentiation in quantitative, and thus potentially adaptive traits, is typically tested in common garden experiments [START_REF] Petit | Gene flow and local adaptation in two endemic plant species[END_REF][START_REF] Leinonen | Local adaptation in European populations of Arabidopsis lyrata (Brassicaceae)[END_REF][START_REF] Bischoff | Testing population differentiation in plant species-how important are environmental maternal effects[END_REF]. Tests for local adaptation require reciprocal transplantation to original collection sites [START_REF] Blanquart | A practical guide to measuring local adaptation[END_REF]. However, common garden experiments comparing the response to the manipulation of key environmental factors allow conclusions on adaptive differentiation and may reveal more specific information on drivers of local adaptation than reciprocal transplantations in the field [START_REF] Lenssen | Local adaptation of the clonal plant Ranunculus reptans to flooding along a smallscale gradient[END_REF][START_REF] Pahl | No evidence for local adaptation in an invasive alien plant: field and greenhouse experiments tracing a colonization sequence[END_REF].

We used the perennial grass Brachypodium retusum (Pers.) P. Beauv. as a model species to test for neutral and adaptive population differentiation at a regional scale in Southern France. The species is dominant in dry grasslands of the Western Mediterranean basin, but also occurs in matorral and open woodlands [START_REF] Caturla | Early post-fire regeneration dynamics of Brachypodium retusum Pers. (Beauv.) in old fields of the Valencia region (eastern Spain)[END_REF][START_REF] Catalan | Phylogeny and Evolution of the Genus Brachypodium[END_REF]. In contrast to studies showing a high colonization capacity of the species [START_REF] Caturla | Early post-fire regeneration dynamics of Brachypodium retusum Pers. (Beauv.) in old fields of the Valencia region (eastern Spain)[END_REF][START_REF] Luis | Fire and torrential rainfall: effects on the perennial grass Brachypodium retusum[END_REF][START_REF] Cassagne | Using a fire propagation model to assess the efficiency of prescribed burning in reducing the fire hazard[END_REF], spontaneous natural regeneration of B. retusum was poor in the La Crau nature reserve (Southern France) 30 years after conversion from arable use to traditional sheep grazing, complicating ecological restoration operations [START_REF] Coiffait-Gombault | Are old Mediterranean grasslands resilient to human disturbances?[END_REF]. We hypothesize that genetic differentiation contributes to the different colonization patterns, and that this differentiation results from adaptation to major environmental factors distinguishing habitats. In the study region, B. retusum occurs in lower calcareous mountain ranges and on more or less decalcified red Mediterranean soils of former riverbeds such as the Crau area. The latter habitat is also characterized by stones (pebbles and cobbles) shaped by the former rivers. These stones were found to favor growth and survival of B. retusum because rhizomes and roots can grow underneath and thus benefit from higher soil moisture during the hot and dry Mediterranean summer [START_REF] Buisson | Limiting processes for perennial plant reintroduction to restore dry grasslands: Perennial plant reintroduction in dry grasslands[END_REF]. Apart from soil characteristics and stone cover, these two habitats show also large differences in soil moisture and grazing intensity. The information on population differentiation and adaptation to environmental factors is crucial to predict colonization processes, but also to improve seed sourcing in ecological restoration [START_REF] Mckay | How local is local?"-a review of practical and conceptual issues in the genetics of restoration[END_REF][START_REF] Vander Mijnsbrugge | A question of origin: Where and how to collect seed for ecological restoration[END_REF].

In our study, we tested B. retusum populations from calcareous and closest red Mediterranean soils replicated at different geographical distances to evaluate whether local-scale differentiation along strong environmental gradients can be as important as large-scale differentiation driven by macroclimate or genetic drift [START_REF] Bischoff | Detecting local adaptation in widespread grassland species-the importance of scale and local plant community[END_REF][START_REF] Hereford | Limits to local adaptation in six populations of the annual plant Diodia teres[END_REF][START_REF] Durka | Genetic differentiation within multiple common grassland plants supports seed transfer zones for ecological restoration[END_REF]. We further analyzed differential responses of these populations to grazing simulated by clipping, to stone cover and to soil moisture. We specifically wanted to answer the following research questions:

(1) At which scales does population differentiation occur: what is the role of local-scale soil type compared with large-scale variation? (2) Is this differentiation adaptive in terms of soil conditions (stone cover), management (grazing simulated by clipping) and humidity (soil moisture, germination)?

MATERIAL AND METHODS

Study species

Ramose false brome (Poaceae, Brachypodium retusum P. Beauv.) is a rhizomatous perennial C3 grass species. Its rhizomes form a dense network under the soil surface, resulting in a high resilience to aboveground disturbance such as wildfires or grazing [START_REF] Caturla | Early post-fire regeneration dynamics of Brachypodium retusum Pers. (Beauv.) in old fields of the Valencia region (eastern Spain)[END_REF]. Vegetative growth is characterized by tillering at the stem base, but also by lateral branching at upper stem parts. Brachypodium retusum is an outcrossing, wind-pollinated species. It flowers from April to July and shows high clonal growth and reproduction when rhizome connections disintegrate. However, rhizomes are short and just allow a slow clonal spread (< 5 cm × year -1 in our experiments).

Collection sites and sampling

The study and sampling sites were situated in the central part of the French Mediterranean climate zone from Montpellier to the lower Durance valley. Seeds were collected from six populations on baserich calcareous soils (Calcaric cambisol) of lower mountain ranges with pH values of 7.8-8.9 (Table 1, Fig. 1). Close to each of these six populations on calcareous grasslands, six populations were also collected from usually decalcified red Mediterranean soils (Haplic cambisols) with pH values of 6.7-8.0. Three additional populations were collected in the Crau nature reserve as the focal site of Mediterranean dry grassland restoration [START_REF] Buisson | Creation of the natural reserve of La Crau: Implications for the creation and management of protected areas[END_REF] in order to obtain a better representation of this largest red Mediterranean soil area of the study region. For further analysis, the four Crau populations (distance: 4-11 km) were pooled to the red Mediterranean soil population of Saint-Martin-de-Crau (SM, Table 1).

Red Mediterranean soils show a higher clay content than calcareous soils, resulting in a better nutrient and water retention. Because of the higher soil fertility, grasslands on red Mediterranean soils were more intensively used (in particular grazing, less frequently mowing) than those on calcareous soils (Table 1). They are also rich in hydrated oxides of iron and aluminum, and in silica whereas calcareous soils are poorer in these elements. Up to 50% of red Mediterranean soil surface is covered by cobbles and pebbles deposited by former rivers. Calcareous soils are less deep on calcareous bedrock and contain numerous angular brittle stones.

The study area is characterized by a typical Mediterranean climate with hot and dry summers, and mild and humid winters. The six sampling sites differ in climatic conditions (Table 1) whereas differences between soil types within sites are small. Annual precipitation depends on the proximity to mountain ranges and is generally higher in the western than in the eastern parts of the study area. Close to the Mediterranean Sea, winter temperatures are higher resulting in a lower frost frequency.

Seed collection and preparation

In each population, all inflorescences of each of 50 1-m² patches were collected. The minimum distance between patches was 10 m and sampling was done in one or several parallel lines depending on population shape (linear versus polygon). Because of clonal growth of the species, the real patch size of individuals is unknown. A distance of 10 m was considered sufficient to obtain samples from different individuals [START_REF] Montalvo | Transplantation of the Subshrub Lotus scoparius: Testing the Home-Site Advantage Hypothesis[END_REF]. On average, one patch comprised 40 inflorescences resulting in a total of 2000 collected inflorescences per population. In germination experiments, the number of seed families may be higher than the number of patches because one patch may comprise more than one individual. This pooling at patch level reduces intrapopulation variation compared to half-sib designs based on individual seed families. In all other experiments, only one seedling per patch was kept for planting, corresponding to classical half-sib designs.

Population differentiation under common conditions

Population differentiation was tested in a growth chamber germination experiment and in a common garden study. Inflorescences were shaken in tubes with glass beads using a vortex mixer to separate seeds from glumes. Transparency under transmitted overhead projector light was used as a criterion to remove empty seeds without embryos. In December 2016, seeds were germinated in peat plugs (4 cm × 4 cm). Four seeds of each sampled patch were sown per plug, and altogether, 40 patches of each population were randomly selected for a total of 160 seeds and 40 plugs per population. The plug trays were placed in a growth chamber (temperature 15/20°C, 12/12 h, night/day) and watered regularly. Coleoptile emergence used as a proxy of germination was recorded every 2-3 days for seven weeks. In March 2016, 24 plugs per population (corresponding to 24 patches or seed families) were randomly selected and transplanted to a previously ploughed experimental site in Avignon-Montfavet. Thus, 48% of the initially collected patches were represented in the field. Seedling number was reduced to one (average size) per plug prior to transplantation. Two plugs per population were planted to random positions within each of 12 blocks (3 m × 3 m each). During the first four months, seedlings were moderately watered. Plots were weeded by hand to limit competition. Adaptation to differences in soil conditions, management, and soil moisture

In May 2016, stones of 10-20 cm length and 7-15 cm width were placed in a circle around half of the seedlings of the common garden experiment to test their effect on plant growth. The stones were taken from the red Mediterranean soil of the Caumont population (Ca, Table 1, Fig. 1). Measurements were the same as in the population differentiation study under common conditions detailed above. A subset of four population pairs (red Mediterranean and calcareous soils × Me, Ca, SR, SM sites) was grown one year earlier in a similar common garden experiment to simulate grazing by clipping established plants. Clipping of second year plants prevents high mortality compared to clipping of young seedlings. The design was similar to the differentiation study. In November 2014, seeds were sown in peat plugs and placed in a greenhouse. In February 2015, two seedlings of each population were planted in each of 12 blocks (24 seedlings in total). Because of low germination, seedling number was only one per block for the Me population from calcareous soil (12 seedlings in total). Plants were moderately watered during the first year, but not in the second year. In April 2016, half of the seedlings were clipped 2 cm above the soil. We measured the same traits as in the differentiation study. Additionally, the number of inflorescences per plant was counted at the end of the season. This was not possible in the population differentiation study under common conditions because plants rarely (<10%) flower in the first year.

The response to soil moisture during germination was analyzed in a controlled experiment manipulating moisture levels. The major objective was to test adaptive differentiation in germination response under dry conditions. Five moisture levels from 0.3 to 20% were adjusted by adding water to plastic boxes filled with 1.15 kg of sand and using a Theta Probe ML2x sensor (Delta-T Devices Ltd, Cambridge, United Kingdom). A subset of 10 populations was chosen for this experiment. Six seeds of each population (four for calcareous population of SM due to limited seed number) were sown in each of 12 boxes and each of the five soil moisture levels; this resulted in a total of 360 seeds per population assigned to 60 boxes. Boxes were closed with a translucent lid to keep initial moisture levels. The boxes were placed in a growth chamber with a 12-h cycle of light and temperature (day: 25°C; night: 12°C). Germination and moisture were measured every three days for a period of six weeks. Germinated seeds were removed. At the end of the experiment, the water potential of the sand was measured using WP4C Dewpoint Water Potentiometer (Decagon Devices, Pullman, Washington, USA). Water potentials were close to 0 (saturation) for the highest and -0.18 MPa for the lowest moisture level (Appendix S1, see Supplemental Data with this article).

Statistical analysis

Linear models were used to analyze population and treatment effects on germination and growth traits. Population was divided into a soil and site component. The two adjacent populations from red Mediterranean and calcareous soils were considered as originating from the same site resulting in a total of six sites and two soil types per site. Thus, each population represented a particular combination of site and soil. Site, soil, and their interaction were included as fixed effects in the models. The effect size of site of origin, soil of origin, and the site × soil interaction was calculated as 100 × factor sum of square/total sum of squares.

The common garden studies were designed as completely randomized block experiments; thus, block was included as an additional random factor. The adaptation treatments (stone cover, clipping) were included as a fixed factor for measurements taken after treatment start. In some cases log10(x + 1) transformation was required to comply with assumptions of linear models (normality, homoscedasticity).

In the germination experiments, final germination percentage in plugs (population differentiation study) and in boxes (soil moisture test) was calculated as response variable. Germination percentage was arcsin (sqrt(x/100)) transformed prior to analysis [START_REF] Quinn | Experimental Design and Data Analysis for Biologists[END_REF]. Because the position of plugs was completely randomized before transplantation to the field, no block effect was considered (germination in the population differentiation study). The soil moisture test on germination was analyzed in a split-plot model. The box effect was included as a random factor nested in moisture level. Accordingly, moisture was tested against the box effect whereas soil, site, and all interactions were tested against the model error.

A Tukey honest significance test (HSD) posthoc test was calculated to analyze differences between treatment levels more in detail if the treatment main effect or interactions were significant. Posthoc tests were preferred over a priori linear contrasts because we had no clear a priori hypothesis on specific adaptations to tested environmental factors. All statistical analyses were run in R (R, version 3.3.1, R Development Core Team (2013)).

RESULTS

Population differentiation under common conditions

A high variation in germination capacity was found between populations ranging from 10% (Ca, red Mediterranean soil) to 75% (Ni, calcareous soil) (Fig. 2). Site and soil of origin as well as their interaction, had a significant influence on germination (Table 2). Populations from calcareous soils and from the westernmost sites germinated significantly better than those from red Mediterranean soils and eastern parts of the study area, respectively. Differences between soil types depended on site of origin, resulting in a significant interaction.

Population differentiation in plant growth was smaller than in germination, but still significant (Table 2, Fig. 2). Differentiation in growth traits did not decrease during the growing season (Appendix S2). The site of origin contributed more to the differentiation than soil of origin. Effect size of site of origin was between 3% of total variation for leaf length and 14% for leaf number whereas effect size of soil of origin was always lower than 1% (Fig. 3). The site × soil interaction explained 1% and 7.7%, respectively, of the model variation. In agreement with the germination results, the westernmost populations from Ni and Mon showed a higher performance in early (leaf length, cumulative leaf length) and late growth traits (tiller length and number, cumulative tiller length, tussock diameter). The soil type effect was only significant for leaf length being higher in populations from red Mediterranean soils than in those from calcareous soils. In several traits, in particular in leaf length, tiller length and plant diameter, a significant site × soil of origin interaction was detected. For the westernmost sites of origin, growth was higher in populations from red Mediterranean soils whereas no differences or even a higher performance on calcareous soils were found for populations from the eastern and central parts of the study area.

Adaptation to differences in soil conditions, management, and soil moisture

Stone cover around the seedlings significantly increased plant height, tiller number, and thus cumulative tiller length, but not cumulative leaf length (Table 2, Fig. 4). The effect of stone cover was the same for plants from the two soil types, while stone cover had contrasting effects on leaf number and cumulative leaf length as well as on tiller number and cumulative tiller length of plants from the different sites. The populations from SR including the SM populations showed a consistent positive response to stone cover, whereas the response of the Me populations was negative or neutral.

Simulated grazing (clipping) significantly reduced plant height, leaf number, tiller number, cumulative tiller length, cumulative leaf length, and inflorescence number (Table 3, Fig. 5). The negative effect on tiller number was the result of cutting underneath lateral branches that were also counted as tillers. While site × treatment interactions were not significant, several significant soil × treatment interactions were detected. Populations from red Mediterranean soils showed significantly higher compensatory growth resulting in smaller differences in plant height, cumulative tiller length, and inflorescences number between clipped and nonclipped plants.

The test of soil moisture effects on germination in different B. retusum populations revealed high germination rates even under dry and a significant decrease of germination under maximum soil moisture (at saturation level: 19.5%, Fig. 6). The highest germination was obtained for intermediate soil moisture levels. With 50% to 80% on average, germination rates were higher than in peat plugs used for the germination test on population differentiation. In agreement with the peat plug germination test, site and soil type effects were significant with a higher germination rate in populations from calcareous soils and from the westernmost sites, respectively (Fig. 6).

The interaction between soil type of origin and soil moisture was not significant and the one between site and soil moisture was only marginally significant suggesting similar adaptation to dry conditions. The marginally significant interaction was explained by a stronger negative response of the Ni populations to drought stress compared with populations from the three other sites.

DISCUSSION

Population differentiation under common conditions

Our common garden study demonstrated genetic differentiation among populations at a regional scale. We found significant differences between populations for germination capacity and several morphological traits such as leaf length, tiller number, and height whereas variation in leaf number per tiller was low. In particular, germination and tiller number are closely linked to seedling recruitment and establishment, thus supporting our hypothesis that different colonization patterns observed at different Mediterranean sites (high: [START_REF] Caturla | Early post-fire regeneration dynamics of Brachypodium retusum Pers. (Beauv.) in old fields of the Valencia region (eastern Spain)[END_REF][START_REF] Luis | Fire and torrential rainfall: effects on the perennial grass Brachypodium retusum[END_REF]low: Coiffait-Gombault et al., 2012) may result from genetic differentiation.

Site of origin better explained this differentiation than soil type. Independent of soil type, the two eastern populations showed significantly lower germination rates and growth than western populations. The significant effect of site is likely the result of climatic differences between eastern and western parts of the study area although our design does not allow us to distinguish random genetic drift from climatic differences as drivers of differentiation. The sites of the western populations (Ni, Mon) are closer to the Cévennes Mountains and receive more rainfall than those of the Eastern populations [START_REF] Ceresetti | Space-time structure of heavy rainfall events: application to the Cévennes-Vivarais region[END_REF]. Climate is known to be a major driver of population differentiation and adaptation [START_REF] Linhart | Evolutionary Significance of Local Genetic Differentiation in Plants[END_REF][START_REF] Etterson | Evolutionary potential of Chamaecrista fasciculata in relation to climate change. I. clinal patterns of selection along an environmental gradient in the great plains[END_REF][START_REF] Ravenscroft | Intraspecific functional differentiation suggests local adaptation to long-term climate change in a calcareous grassland[END_REF]. Water availability and temperature have a particularly strong effect on plant fitness leading to population differentiation in the case of spatial variation [START_REF] Shah | Interaction of drought and high temperature on photosynthesis and grain-filling of wheat[END_REF][START_REF] Xu | Combined effects of water stress and high temperature on photosynthesis, nitrogen metabolism and lipid peroxidation of a perennial grass Leymus chinensis[END_REF]. [START_REF] Yuan | Environmental conditions and genetic differentiation: what drives the divergence of coexisting Leymus chinensis ecotypes in a large-scale longitudinal gradient?[END_REF] found for example that population differentiation in the grass species Leymus chinensis was driven by precipitation, temperature, and soil water content. A mainly soil moisture driven differentiation was also shown for Boechera stricta [START_REF] Lee | Quantifying Effects of Environmental and Geographical Factors on Patterns of Genetic Differentiation[END_REF]. However, the authors also found that isolation by distance due to stochastic processes such as genetic drift considerably contributes to population differentiation. Adaptive differentiation and isolation-bydistance effects often act at the same time on the genetic structure of populations and their contribution to differentiation may be difficult to distinguish [START_REF] Michalski | Assessment of provenance delineation by genetic differentiation patterns and estimates of gene flow in the common grassland plant Geranium pratense[END_REF][START_REF] Durka | Genetic differentiation within multiple common grassland plants supports seed transfer zones for ecological restoration[END_REF]. Nevertheless, shoot and leaf growth (number and length) were lower in populations from dryer Eastern sites suggesting a selection for growth forms with smaller and less numerous leaves and shoots that are better protected against drought stress [START_REF] Liancourt | Competition and a short growing season lead to ecotypic differentiation at the two extremes of the ecological range[END_REF][START_REF] Pérez-Ramos | Tradeoffs between functional strategies for resource-use and droughtsurvival in Mediterranean rangeland species[END_REF]. The significant site effects on the cited growth traits were stable throughout the whole growing season and did not decrease with time suggesting a minor influence of environmental maternal effects on our results. Environmental maternal effects may contribute to population differentiation, but particularly affect early traits and vanish during plant development [START_REF] Roach | Maternal Effects in Plants[END_REF][START_REF] Donohue | Completing the cycle: maternal effects as the missing link in plant life histories[END_REF].

Soil type is also known to select for population differentiation in plants [START_REF] Linhart | Evolutionary Significance of Local Genetic Differentiation in Plants[END_REF]Raabova et al., 2011). In our study, we detected a significant soil type of origin effect on differentiation at a small scale of 10-20 km, but only for germination capacity and leaf length. Germination is under stronger control of environmental maternal effects than growth traits affecting population differentiation in such early traits [START_REF] Roach | Maternal Effects in Plants[END_REF][START_REF] Bischoff | Testing population differentiation in plant species-how important are environmental maternal effects[END_REF]. However, it is also possible that selection for clonal growth under higher grazing pressure of red Mediterranean soils has reduced allocation of resources to seeds compared with less grazed calcareous sites. Although providing some evidence for soil of origin effects, our results do not confirm earlier studies showing a high small-scale differentiation along strong environmental gradients [START_REF] Lenssen | Local adaptation of the clonal plant Ranunculus reptans to flooding along a smallscale gradient[END_REF][START_REF] Bischoff | Detecting local adaptation in widespread grassland species-the importance of scale and local plant community[END_REF][START_REF] Raabová | Ecological rather than geographic or genetic distance affects local adaptation of the rare perennial herb, Aster amellus[END_REF][START_REF] Bischoff | Scales and drivers of local adaptation in Brassica nigra (Brassicaceae) populations[END_REF]. Because no specific barriers prevented pollen dispersal at this scale, gene flow may have counteracted local adaptation to soil type [START_REF] Kawecki | Conceptual issues in local adaptation[END_REF].

Adaptation to differences in soil conditions, management, and soil moisture

Brachypodium retusum populations showed a differential response to stone cover and clipping, and to a lesser degree to soil moisture in the germination phase.

Our study demonstrated a clearly positive effect of stone cover on the growth of B. retusum, in particular on plant height, tiller number, and plant diameter. Stone cover increases soil moisture by shading and water condensation, and protects lower stem parts and roots from strong solar radiation and heat [START_REF] Peters | Nurse rocks are more important than nurse plants in determining the distribution and establishment of globose cacti (Mammillaria) in the Tehuacán Valley, Mexico[END_REF][START_REF] Buisson | Limiting processes for perennial plant reintroduction to restore dry grasslands: Perennial plant reintroduction in dry grasslands[END_REF]. A differential response of Aizoaceae species demonstrated that plants might be locally adapted to type and density of stone cover [START_REF] Ellis | Coexistence and differentiation of 'flowering stones': the role of local adaptation to soil microenvironment[END_REF]. In our study, however, the differential response was not related to the dominant soil types, but to site of origin. Although the texture and depth of calcareous soils and red Mediterranean soils are different, limestone layers and pieces close to the soil surface may have similar effects on microclimate. Thus, the differential response to stone cover cannot be related to differences in local environmental conditions at collection sites.

As expected, clipping largely reduced plant height, tiller number, and leaf number throughout the whole season resulting in a highly reduced reproduction. The effect was less strong in populations from red Mediterranean soils than in populations from calcareous sites suggesting an adaptation to higher grazing intensity. Grazing exerts strong directional selection and may result in genetic differentiation and local adaptation to long-term differences in grazing intensity [START_REF] Linhart | Evolutionary Significance of Local Genetic Differentiation in Plants[END_REF][START_REF] Hufford | Local Adaptation and the Effects of Grazing on the Performance of Nassella pulchra: Implications for Seed Sourcing in Restoration[END_REF].

Higher compensatory growth may be one important adaptation of grasses to higher grazing intensity [START_REF] Sarmiento | Adaptive strategies of perennial grasses in South American savannas[END_REF]. Although we do not exactly know whether our collection sites show such long-term differences, red Mediterranean soils were more grazed in the past than calcareous soils because slopes are less steep and rocky, and soils are more fertile-both of these attributes favor sheep grazing [START_REF] Verheye | Mediterranean soils. Land Use, Land Cover and Soil Sciences-Volume VII: Soils and Soil Sciences-296[END_REF].

Drought stress reduces germination, and usually plant species show a positive correlation between soil moisture and germination rates [START_REF] Bochet | Soil water availability effects on seed germination account for species segregation in semiarid roadslopes[END_REF][START_REF] Thomas | The interaction of temperature, water availability and fire cues regulates seed germination in a fire-prone landscape[END_REF]. The ability to germinate fast under moderate-to-low soil moisture may be crucial for species colonizing dry environments [START_REF] Schütz | Germination requirements and seedling responses to water availability and soil type in four eucalypt species[END_REF][START_REF] Bochet | Soil water availability effects on seed germination account for species segregation in semiarid roadslopes[END_REF]. All tested B. retusum populations showed a high germination capacity even under the lowest soil moisture level (-0.18 MPa = 0.3% volumetric water content), but a considerable drop under saturated conditions. The observed germination pattern differs from the results of [START_REF] Bochet | Soil water availability effects on seed germination account for species segregation in semiarid roadslopes[END_REF] who found a strong decrease in B. retusum germination already at water potentials of -0.01 and -0.05 MPa. Although the results may depend on age of seeds and specific experimental conditions this difference suggests genetic variability in germination response to soil moisture. Environmental maternal effects may contribute to population differentiation in germination, but there is little evidence for a contribution to adaptive responses through transgenerational plasticity [START_REF] Galloway | Transgenerational Plasticity Is Adaptive in the Wild[END_REF][START_REF] Bischoff | Testing population differentiation in plant species-how important are environmental maternal effects[END_REF]. In our study, the response to soil water potential was not much different among populations. Only the population of one of the two sites with the most humid climate (Ni) showed a stronger decline under dry conditions suggesting weak adaptation to climate-mediated soil moisture.

CONCLUSIONS

Significant population differentiation was found at a regional scale of 30-200 km. Differentiation between major soil types representing strong environmental differences were found at even smaller scales (10-20 km). However, differentiation between soil types was weaker than differentiation among sites, suggesting a major influence of climate and/or genetic drift on population differentiation.

Stone cover increased B. retusum growth and the differential response of populations suggests adaptation to this environmental factor. However, we could not relate this response to the home environment of tested populations. Further experiments are required to identify potential drivers of adaptation to stone cover (density at the surface or in different soil layers). The response to clipping was also clearly different among populations, and the differential response could be linked to the major soil type suggesting that different grazing intensity in the past has selected for different ecotypes. Brachypodium retusum showed high germination under dry soil conditions, but there was only weak evidence for a differential response of tested populations. The study demonstrated that population differentiation may well explain different colonization patterns observed in different Mediterranean regions [START_REF] Caturla | Early post-fire regeneration dynamics of Brachypodium retusum Pers. (Beauv.) in old fields of the Valencia region (eastern Spain)[END_REF][START_REF] Luis | Fire and torrential rainfall: effects on the perennial grass Brachypodium retusum[END_REF][START_REF] Coiffait-Gombault | Are old Mediterranean grasslands resilient to human disturbances?[END_REF]. Because observed differences may largely influence seedling establishment in Mediterranean dry grasslands, plant origin needs to be carefully considered in ecological restoration involving planting or sowing of B. retusum.

FIGURE 2 .

 2 FIGURE 2. Population differentiation decomposed in collection site and soil type within collection site (mean values ±SE, see Table 2 for statistics). A, Final germination. B, Cumulative leaf length. C, Cumulative tiller length. D, Plant diameter. Different lower case letters indicate significant differences (P < 0.05).

FIGURE 3 .

 3 FIGURE 3. Percentage of variation explained by collection site, soil type within collection site and site × soil interaction for growth traits (May).

FIGURE 4 .

 4 FIGURE 4. Differential response to stone cover (mean values +SE) between Brachypodium retusum populations in southern France. A, Plant height. B, Cumulative leaf length. C, Cumulative tiller length. D, Plant diameter.

FIGURE 5 .

 5 FIGURE 5. Differential response of a subset of populations to simulated grazing (mean values ±SE). Since site × clipping interaction was never significant only soil type × clipping combinations are presented. A, Plant height. B, Cumulative leaf length. C, Cumulative tiller length. D, Inflorescence number.

FIGURE 6 .

 6 FIGURE 6. Differential germination response of a subset of populations to soil moisture (mean values ±SE). Since the soil of origin x soil moisture interaction was never significant only site of origin effects are presented. Soil moisture (F 4,55 = 53.71, P < 0.001), site of origin (F 3,505 = 27.96, P < 0.001), soil of origin (F 1,505 = 33.86, P < 0.001), site of origin × soil of origin × moisture interaction (F 12,505 = 1.726).

TABLE 1 .

 1 Soil characteristics, grazing intensity, annual precipitation, mean annual air temperature (including January and July mean in brackets), and the number of frost days (T<0) at collection sites of the twelve populations. Sites are ordered from east to west. Climate data are based on daily averages from the nearest meteorological stations. *, Irregular mowing.Plant height, diameter, tiller number, leaf number, and the length of the longest leaf were measured once a month. Young B. retusum plants allocate resources to tillering and leaf growth. Later, tillers elongate and older leaves die. Therefore, we calculated cumulative leaf length as an aboveground biomass estimate for young seedlings (leaf length × number of leaves) and cumulative tiller length for older plants (plant height × number of tillers).

							Precipitation	Temperature (°C)
	Site code	Site name	Soil	Coordinates	pH Grazing	(mm*a -1 )	Mean [Jan/Jul]	Days T<0
	Me	Mérindol	Calcareous	43.756456 N, 5.199526 E	8.0	No*	583	14.5 [6.0/23.7]	49
			Red Mediterranean 43.750425 N, 5.210495 E	8.0	Low
	Ca	Caumont	Calcareous	43.898827 N, 4.941070 E	8.2	No	676	14.6 [5.9/24.2]	33
			Red Mediterranean 43.898319 N, 4.932395 E	6.9	No
	SR	Saint-Rémy-de-Provence Calcareous	43.768047 N, 4.831898 E	7.8	No	625	15.1 [6.6/24.3]	30
			Red Mediterranean 43.798952 N, 4.841659 E	6.8	High
	SM	Saint-Martin-de Crau	Calcareous	43.705354 N, 4.797561 E	8.1	No	517	15.3 [7.4/24.4]	24
			Red Mediterranean 43.567465 N, 4.834645 E	6.8	High
	Ni	Nîmes	Calcareous	43.784700 N, 4.259389 E	8.0	No	701	14.9 [6.7/24.1]	34
			Red Mediterranean 43.722868 N, 4.339264 E	7.8	No*
	Mon	Montpellier	Calcareous	43.645416 N, 3.691638 E	8.9	No*	870	14.1[6.4/23.2]	51
			Red Mediterranean 43.607580 N, 3.449269 E	6.7	No*

FIGURE 1. Position of six sampling sites each comprising two populations from different soil types.

TABLE 2 .

 2 Differentiation in Brachypodium retusum populations from different sites and soil types, and differential population response to stone cover treatment. Fvalues and significance levels from two-way ANOVA before stone treatment (growth traits: May) and three-way ANOVA afterwards (July) including the stone effect.

	* P

conditions (0.3% volumetric water content)

TABLE 3 .

 3 Differential response in a subset of populations to simulated grazing. F-values and significance levels from three-way ANOVA on site, soil (within site) and clipping effects. Vegetative traits measured in spring (May), inflorescence number in June of the second growing season. * P < 0.1; ** P < 0.05; *** P < 0.01; **** P < 0.001; NS, not significant.

df Leaf number Leaf length Plant height Tiller number Cumulative leaf length Cumulative tiller length Plant diameter Inflorescence number

	Site	3,132	1.20 NS	1.13 NS	0.27 NS	0.63 NS	1.47 NS	0.64 NS	1.30 NS	2.48 *
	Soil	1,132	0.02 NS	9.43 ***	3.16 *	0.05 NS	4.13 **	2.87 *	1.54 NS	3.91 *
	Clipping	1,132	12.31 ****	3.19 *	220.98 ****	10.48 ***	11.46 ****	211.72 ****	7.23 **	58.39 ***
	Soil × Site	3,132	0.67 NS	2.12 NS	0.81 NS	0.16 NS	2.05 NS	1.18 NS	0.67 NS	1.98 NS
	Clipping × Site	3,132	0.66 NS	0.49 NS	0.81 NS	0.02 NS	0.99 NS	1.49 NS	0.94 NS	0.23 NS
	Clipping × Soil	1,132	0.08 NS	0.09 NS	3.93 **	2.17 NS	0.12 NS	3.07 *	0.48 NS	3.37 *
	Clipping × Soil × Site	3,132	1.12 NS	0.61 NS	1.65 NS	0.76 NS	0.43 NS	1.08 NS	0.40 NS	0.29 NS

ACKNOWLEDGEMENTS

We thank Stéphane Ruy for water potential measurements and Vivien Carrere, Lara Amorsi, Hervé Ramone, Sébastien Dailly, Julien Formento, Hugo Lezarme, Kevin Maurin, Cannelle Moinardeau, Julie Chenot, Anne Aurière, Lorenzo Ottaviani, Jean-François Alignan, Olivier Blight, and Elise Buisson for field assistance. We are grateful to the Conservatoire des Espaces Naturels (CEN-PACA) for the permit to collect seeds in the Nature Reserve "La Crau". This work was supported by the Région PACA (Ph.D. grant) and the University of Avignon (Excellence program). We also thank two anonymous reviewers for their helpful and constructive comments.

SUPPORTING INFORMATION

Additional Supporting Information may be found online in the supporting information tab for this article.

LITERATURE CITED


M Khalid

L Pénard

E Mémin

Optical Flow For Image-Based River Velocity Estimation

Keywords: Optical Flow, River Velocimetry, PIV, LSPIV

 for a review). Our optical flow is based on the scalar transport equation and is augmented with a weighted diffusion term to compensate for small scale (non-captured) contributions. Additionally, since there is no ground truth data for such type of image sequences, we present a new evaluation method to assess the results. It is based on trajectory reconstruction of few Lagrangian particles of interest and a direct comparison against their manually-reconstructed trajectories. The new motion estimation technique outperformed traditional optical flow and PIV-based methods.

Introduction

To estimate the velocity of the free surface of a river using image-based techniques, there are two main issues to be addressed. The first is the 2D motion estimation in image space and the second issue is about how to transform 2D image estimations back to 3D. It should also handle the scaling of these estimations to their true physical scale. The metric surface velocity could be then obtained using the time between images. It is a valuable piece of information to compute rivers discharge using the well-known velocity-area method. The work presented here focuses on the first issue. We specifically target rivers with noticeable discharge or urban inundations. In these settings, the free surface manifests a visible and probably turbulent translatory motion.

For many years, PIV (Particle Image Velocimetry) (Adrian, 1991) was the only technique available for fluid motion estimation in image sequences. PIV computes the displacement between windows of a certain size in two consecutive images depicting spatio-temporal evolution of a fluid. These windows are tracked assuming a translational motion between images using some measure of similarity that is computed on them. To compute better similarity scores, the flow has to be well-seeded with particles. As a consequence, PIV faces strong difficulties on scalar images in which many areas have low intensity gradient. In many cases, there is a necessary postprocessing step to correct or remove spurious vectors, inevitably generated in regions with low intensity gradient. The post-processing step is also used to interpolate the sparse vector field to generate a denser one. In addition, choosing a suitable window size is a tricky task, one needs bigger window to accumulate more information for better score computation, but then the bigger the window, the more likely it contains complex motions far from the window's translational motion assumption. Despite its wide-spread usage specially in controlled lab environments, PIV lacks sound physics foundation, it treats all image sequences equally, regardless of the nature of the object in motion. Furthermore, it only computes sparse estimations, many of them might be discarded in the post-processing phase. In their seminal paper, [START_REF] Fujita | Large-Scale Particle Image Velocimetry for Flow Analysis in Hydraulic Engineering Applications[END_REF] took PIV outside the lab to measure the velocity of rivers in ortho-rectified image sequences. The ortho-rectification applies geometric transformation to produce new set of images in which the perspective effect is removed. The ortho-rectified sequence is produced in such a way that their real spatial resolution is now known. The image motion estimation technique (PIV in this case) is then applied to the newly created sequence. Hence, the name of their method LSPIV (Large Scale Particle Image Velocimetry) which became the benchmark for image-based river velocimetry [START_REF] Jodeau | Application and Evaluation of ls-piv Technique for the Monitoring of River Surface Velocities in High Flow Conditions[END_REF][START_REF] Dramais | Advantages of a Mobile LSPIV Method for Measuring Flood Discharges and Improving Stage-Discharge Curves[END_REF][START_REF] Muste | Capabilities of Large-Scale Particle Image Velocimetry to Characterize Shallow Free-Surface Flows[END_REF].

If we closely examine the original problem as depicted in Figure (1.1), the motion in 3D space is captured by a camera and as a result, a 2D motion is observed in images. The image-based velocimetry problem is exactly solved by taking the reverse route. We can distinguish the two aforementioned issues or steps that need to be solved along this route. Most LSPIV works however mix the two steps together which makes it difficult to determine the exact sources of uncertainty and the effects of each step error to the overall error. It is a valid idea to take the route shown by the blue arrows in Figure (1.1) and treat each issue individually. Recently, [START_REF] Patalano | Rectification of image velocity results (RIVeR): A simple and user-friendly toolbox for large scale water surface particle image velocimetry (PIV) and particle tracking velocimetry (PTV)[END_REF] applied the ortho-rectification only to the results but the motion estimation was separately applied to the original obliquely-viewed images. This paper focuses only on the first step, that is, motion estimation of river surface in image space.

Due to PIV limitations cited above, we found the so-called Optical Flow [START_REF] Horn | Determining Optical Flow[END_REF] to be well adapted to fluids. A nice review paper on the application of optical flow to fluids could be find in [START_REF] Heitz | Variational Fluid Flow Measurements from Image Sequences: Synopsis and Perspectives[END_REF].

Previous work

Stemmed from computer vision community in the early 80s due to [START_REF] Horn | Determining Optical Flow[END_REF] work. Optical flow consists of an energy functional to be minimized. It combines the information that could be extracted from images (called data term) and the expected solutions set (called regularization term) in a concise and unambiguous way:

Ω ∂I ∂t + ∇I• ω 2 + α ∇ω 2 dX (1.1)
where ∇ represents the spatial gradient for both image intensity function I and the 2D velocity field ω = [u, v]. A parameter α is introduced to define the trade-off between the data term and the regularization term. The first quadratic term (data term) is called the brightness constancy (BC). It assumes that points conserve their brightness value captured by the camera while moving in a small time interval. The second term is trying to minimize velocity gradient globally, assuming that neighbouring points belong normally to the same surface and tend to move together. This assumption resembles the one made in PIV that all points in the window have the same velocity. However, its formulation in optical flow is more relaxed, velocities of points on the same neighbourhood could vary, but this variation should be limited. This model witnessed many improvements over the years to overcome its limitations. First of all, the BC assumption used in the data term in (1.1) is derived using Taylor approximation to linearize the original non-linear pattern matching BC, I(X(t), t) = I(X(t + 1), t + 1). The linearized model assumes very small motion which is indeed not the case in general. To deal with large displacements, the use of different image resolutions in a hierarchical, coarse-to-fine way (called image pyramid) has been introduced [START_REF] Bergen | Hierarchical Model-Based Motion Estimation[END_REF]. This strategy is improved by using warping steps in which the second image is brought closer to the first image using the estimated velocity field of the previous pyramid level (Mémin andPérez, 1998, 2002). Consequently, only a small increment is sought at each level.

Secondly, the quadratic penalization is known to be sensitive to outliers. This is particularly seen in image regions with discontinuities where we might have two different objects moving in different directions at different speeds. The quadratic penalization tends to average these two motions which results in a solution clearly not in agreement with the real motion at the discontinuity area and around it. [START_REF] Black | The Robust Estimation of Multiple Motions: Parametric and Piecewise-Smooth Flow Fields[END_REF] proposed to use robust penalization methods to cope with outliers. This procedure gave better results i.e. sharper discontinuities.

Lastly, the BC assumption might not hold in many cases due to occlusions, imaging acquisition quality and/or illumination changes caused by surface reflections for instance. Therefore, it is a good practice to use a robust function for the data term as well [START_REF] Mémin | A Multigrid Approach for Hierarchical Motion Estimation[END_REF]. Other notable improvements were proposed by [START_REF] Mémin | A Multigrid Approach for Hierarchical Motion Estimation[END_REF], [START_REF] Brox | High Accuracy Optical Flow Estimation Based on a Theory for Warping[END_REF] in which authors used the original non-linear version of BC. The resulting non-linear equation is solved iteratively using a multi-grid optimization technique [START_REF] Mémin | A Multigrid Approach for Hierarchical Motion Estimation[END_REF] or a Graduated-Non-Convexity (GNC) [START_REF] Blake | Visual Reconstruction[END_REF] with two fixed point iterations optimization technique. Median filtering of intermediate flow (between pyramid levels) showed very good performance regarding outliers rejection. [START_REF] Sun | A Quantitative Analysis of Current Practices in Optical Flow Estimation and the Principles Behind Them[END_REF] elaborated on this and showed that it is a different energy function when median filtering is utilized.

The functional (1.1) was designed for rigid or quasi-rigid motions as they are the most encountered in natural scenes. Fluids in general have more complex motion patterns. [START_REF] Liu | Fluid Flow And Optical Flow[END_REF] formally established the relationship between optical flow and fluid flow based on the perspective projection of the transport equation on the 2D plane and found that optical flow is proportional to the path-averaged velocity field weighted with a relevant field quantity. This gave the physical foundation needed since the BC assumption is not based on any physical principle. Higher order regularization has been proposed before by [START_REF] Suter | Motion Estimation and Vector Splines[END_REF] based on the divergence and the curl of the vector field. However, minimizing such high order terms is known to be problematic. [START_REF] Corpetti | Fluid Experimental Flow Estimation Based on an Optical-Flow Scheme[END_REF] proposed a relaxed and easier term to minimize the gradient of divergence and curl. They also suggested a data term based on the Integrated Continuity Equation. They also showed that, from an optimization point of view, the L 2 regularization in equation (1.1) is equivalent to penalizing divergence and curl equally, which might be a problem depending on the fluid at hand. [START_REF] Heitz | Dynamic Consistent Correlation-Variational Approach for Robust Optical Flow Estimation[END_REF] suggested to refrain from multi-resolution strategy for PIV applications because the smoothing and sub-sampling of images will suppress the particles in the coarsest resolution levels. They introduced a model that combines local information (PIV) with a global dense optical flow formulation. A very promising line of research consists of deriving new data term based on LES (Large Eddy Simulation) of the transport equation [START_REF] Cassisa | Subgrid Scale Formulation of Optical Flow for The Study of Turbulent Flow[END_REF][START_REF] Chen | Optical Flow for Incompressible Turbulence Motion Estimation[END_REF]. The latter used a divergence norm penalization term as regularization to conserve the curl of the vector field. This derivation from transport equation separates the effects of information captured on the image pixels' grid and the contribution of the small non-captured scales (subgrid). This paper belongs to this line of research and proposes a simple way to model the small scales contribution in river flows.

In section 2 we present the necessary derivation for the optical flow data term followed by the regularization term in section 3. In section 4, we detail how the energy functional is discretized and minimized. An evaluation of the model is done on image sequences of rivers with comparisons against conventional optical flow model and PIV/LSPIV is presented in section 5. We then conclude in section 6. This paper extends and elaborate in greater details a previous work published at a conference [START_REF] Khalid | Application of Optical Flow for River Velocimetry[END_REF]. Main differences amongst others consist of an implementation section that makes the paper more self-contained, direct comparison with LSPIV on ortho-rectified images, rigorous parameters fine-tuning and more results and image sequences.

Data term

For the specific case of rivers sequences, we make the assumption that the image intensity function I is related to some passive scalar field concentration C. They are related by a perspective transformation that projects a real world quantity into the observed image plane forming the intensity function I. Let us examine the scalar transport equation:

∂C ∂t + ∇ • (Cω) - 1 Re Sc ∆C = 0 (2.1)
where ∇• is the divergence operator, ∆ is the Laplacian operator, Re and Sc are the Reynolds and Schmidt numbers, respectively. The transport equation above links the scalar quantity C to the sought velocity vector field ω. However, in order to be able extract the velocity, the function C has to contain information till Kolmogorov scale. Since C is captured by a camera, the pictured image corresponds indeed to a smooth filtered version of the scalar quantity. It is hence a large-scale representation of the scalar in which the small-scale contributions are omitted. In the case of fluid flows the small-scale effects can not be neglected. Their action on the large-scale drift component must be modeled [START_REF] Cui | A New Subgrid Eddy-Viscosity Model for Large-Eddy Simulation of Anisotropic Turbulence[END_REF]. Following [START_REF] Cassisa | Subgrid Scale Formulation of Optical Flow for The Study of Turbulent Flow[END_REF][START_REF] Chen | Optical Flow for Incompressible Turbulence Motion Estimation[END_REF], we propose a LES (Large Eddy Simulation) decomposition of the transport equation to model the small scales contributions in river sequences. The expression ∇ • (Cω) could be divided into an observed (∇ • (Iω)) and non-observed (∇ • (τ )) parts, leading to:

∂I ∂t + ∇ • (Iω) + ∇ • (τ ) - 1 Re Sc ∆I = 0 (2.2)
in image space. The molecular diffusion term 1 Re Sc ∆I is normally neglected. We also assume that the incompressibility condition for water still holds on the 2D plane (that the divergence of the 2D velocity field is zero). When applying the incompressibility condition and ignoring the molecular diffusion term, equation 2.2 becomes:

∂I ∂t + ∇I • ω + ∇ • (τ ) = 0 (2.3)
We see that the BC assumption appears again in addition to the new subgrid term, which means conventional optical flow is consistent with this physics-based derivation. As suggested by [START_REF] Cassisa | Subgrid Scale Formulation of Optical Flow for The Study of Turbulent Flow[END_REF], the non-observed term is considered related to turbulent viscosity τ = -D t ∇I where D t is a turbulent diffusion coefficient. We opted for a simple model to estimate D t using previous velocity estimations within the sequence and/or the pyramid levels of image pair at hand. We feed these velocity estimations to Prandtl mixing-length model (Prandtl, 1925). This model uses the streamwise velocity u and the Mixing Length l to estimate the turbulent viscosity ν t = l 2 du dy . This quantity is relevant for the computation of D t = νt sct where sc t is the turbulent Schmidt number which has an empirical value usually determined experimentally. It has been reported that sc t has widespread values between 0.2 and 3 in general [START_REF] Tominaga | Turbulent Schmidt Numbers for CFD Analysis with Various Types of Flowfield[END_REF] while other authors very recently suggested values around 1 to be optimal on a tracer transport study [START_REF] Gualtieri | On the Values for the Turbulent Schmidt Number in Environmental Flows[END_REF]. Accordingly, D t is considered equal to ν t . The Mixing Length l is defined as the distance traversed by a fluid parcel before it becomes blended in with neighbouring masses. To the best of our knowledge, there is no clear way to predict this value from images. It is taken here to be 1 because the differential optical flow formulation assumes infinitesimal displacements that don't exceed one pixel.

Due to sources or sinks in the fluid, the imaged surface is prone to intermittent changes because of out of plane (depth) motions. Some regions may rise up or sink down during the temporal evolution of the fluid causing changes in the intensity function. As a result, the brightness consistency assumption might not hold in many locations. Thus, a robust function is chosen for the data term. In practice we mainly use the Lorentzian function [START_REF] Black | The Robust Estimation of Multiple Motions: Parametric and Piecewise-Smooth Flow Fields[END_REF] but there are many others that could be used. The data term reads:

Ω Ψ ∂I ∂t + ∇I • ω + ∇ • (-D t ∇I) 2
where Ψ is the robust function. This could be further simplified as:

Ω Ψ ∂I ∂t + ∇I • ω -D t ∆I 2 (2.4)
which is the proposed data term.

Regularization

Regarding the regularization, two main choices were made. Firstly, it can be noticed that usually the river flows of interest do not exhibit strong large-scale eddies motion. Rivers with such eddies have indeed less interest in river velocimetry applications. There might be however many small scale vortices but due to river velocity and/or image acquisition speed, most of them are transported by the largescale quasi-translational motion. As a result, we opted for a first order gradient penalization like in functional (1.1). As stated earlier, this regularization penalizes the divergence and the curl of the flow field equally. On the one hand, penalizing the divergence is beneficial to enforce the 2D incompressibility assumed. On the other hand, penalizing a non existent or a weak property like the curl will reject all solutions with strong vorticity. For the case of floods image sequences, their strong translatory nature does not promote chaotic motions. In other words, the turbulence does not occur at individual pixels. There is always this notion of groups of pixels that move together and in which case this regularization is optimal. In addition, in uniform intensity areas, the image intensity gradient diminishes and it is the regularization that takes the lead over the data term since the latter depends on image gradient. These uniform areas should then be moving together for this regularization to be relevant. We argue that in uniform intensity image areas, it is unlikely to have very different motions, otherwise, because of water and flow properties including velocity, there would be a mixing phenomenon that would eventually disrupt the uniformity of the intensity function in that specific area. The reverse logic should apply, if the area is uniform in intensity, it is more likely that it contains similar velocity vectors i.e. points move together, in which case this regularization is definitely optimal. The second choice made is to take a quadratic function instead of robust one for the regularization. The reasoning behind this is: previously robust functions are utilized in the regularization of rigid motion sequences to improve the results on discontinuities, these discontinuities are due to different objects motions and/or different depth differences, etc... However, there is only one object of interest here, which is the river water surface. Though complex, its underlying motion is not discontinuous and there is no reason to introduce penalty terms enforcing the apparition of discontinuities between two adjacent points. The quadratic functions behavior is exactly suitable in this prospect as it tends to average the two different motions in the neighbourhood of interest. The regularization reads:

Ω α ∇ω 2 dX
The final model equation after combining the two terms is:

E(ω) = Ω Ψ ∂I ∂t + ∇I • ω -D t ∆I 2 + α ∇ω 2 dX (3.1)
We name the model SGSD (SubGrid Scale Diffusion).

Implementation

In practice, the original non-linear displacement model i.e. I(X(t + 1), t + 1) -I(X(t), t) is utilized in a coarse-to-fine fashion with warping as in (Brox et al., 2004, Mémin and[START_REF] Mémin | A Multigrid Approach for Hierarchical Motion Estimation[END_REF]. Not only to allow the estimation of large displacements, but also because the diffusion coefficient D t depends on the velocity computed on previous coarser pyramid level. The linear BC used earlier is only a Taylor approximation to this model while assuming small motions. The velocity vector field is initialized using conventional optical flow at the coarsest level (which in its turn initialized by a zerovalued vector field). The estimated vector field is propagated to the next finer level after applying the necessary scaling. The vector field is then used to warp the second image of the finer level towards the first image of the same level. Only a small increment of the vector field is sought at this level ω k+1 = ω k + dω k+1 , implying that we already know ω k and we look only for the increment dω k+1 . Equation (3.1) could be minimized using Euler-Lagrange equations. However, an Iterative Reweighted Least Squares IRLS approach, previously used by [START_REF] Mémin | A Multigrid Approach for Hierarchical Motion Estimation[END_REF], is equivalent to the variational Euler-Lagrange as shown by [START_REF] Liu | Beyond Pixels: Exploring New Representations and Applications for Motion Analysis[END_REF], yet simpler to derive while working on the discrete equivalent of equation ( 6). Let

I t = I(X + ω) -I(X), I x = ∂ ∂x I(X + ω), I y = ∂ ∂y I(X + ω)
Since we have introduced dω earlier, an expression of an incremental BC I(X + ω + dω) -I(X) could be linearized using Taylor expansion as:

I t + I x du + I y dv (4.1)
where dω = (du, dv). Equation (4.1) is equivalent to the linear BC in equation (1.1) but it is derived on a pyramid level after warping. The difference of this approach with original Horn and Schunck (HS) optical flow is that while both rely on a linearized expression, for HS this linearization is applied only once. While in equation (4.1) the linearization is performed successively on every pyramid level at least once. The discretized version of equation (3.1) applied at a given pyramid level reads:

E(du, dv) = X Ψ δ T X (I t + I x du + I y dv -D t ( I)) 2 + α[ δ T X F x (u + du) 2 + δ T X F y (u + du) 2 + δ T X F x (v + dv) 2 + δ T X F y (v + dv) 2 ] (4.2)
where δ T X is a column vector with all zeros except on the position X and F * ; * {x, y} are derivative filter matrices in the direction of the subscript. u, v, du and dv are all vectorized and I * ; * {x, y, t} are all diagonalized, for instance diag[I x ]. The goal is to find du and dv that minimizes the gradient dE ∂du ; dE ∂dv = 0. Using matrix calculus we have:

∂E ∂du = 2 X Ψ I x δ X δ T X I x du + I x δ X δ T X (I t + I y dv -D t,u ∆I) + F T x δ X δ T X F x + F T y δ X δ T X F y (u + du) (4.3)
where Ψ is the derivative of the robust function with respect to its parameters. Since δ X δ T X will add up to identity, some arrangement of the terms will yield:

∂E ∂du = 2 Ψ I 2 x + αL du + Ψ I x I y dv + Ψ (I x I t -I x (D t,u ∆I)) + αLu (4.4)
where L is a Laplacian filter defined as F T x IDF x + F T y IDF y , ID being the identity matrix. The identity here comes from the quadratic function, if we use a robust function, we would have a diagonal matrix composed of the robust weights instead. An equation for dv could be obtained in the same way.

∂E ∂dv = 2 Ψ I x I y du + Ψ I 2 y + αL dv + Ψ (I y I t -I y (D t,v ∆I)) + αLv (4.5)
The only unknowns beside dv and du are the two turbulent diffusion coefficients, these could now be simply estimated using previous velocity estimates

D t,u = ∂ ∂y u D t,v = ∂ ∂x v (4.6)
IRLS procedure considers the derivative of the robust function as a weight to an ordinary least squares problem. The convergence is achieved when no significant changes are observed in the weights or when a maximum number of iteration is reached. Two fixed point iterations are finally performed. In the inner loop, we solve for du and dv while continuously linearizing the non-linear model and updating the weight Ψ as shown in Algorithm 1. On convergence, the outer loop updates u and v and warp again the second image towards the first using the updated vector field. The final linear operator reads:

Ψ I 2 x + αL Ψ I x I y Ψ I x I y Ψ I 2 y + αL du dv = - Ψ (I x I t -I x (D t,u ∆I)) + αLu Ψ (I y I t -I y (D t,v ∆I)) + αLv (4.7)
As in [START_REF] Brox | High Accuracy Optical Flow Estimation Based on a Theory for Warping[END_REF][START_REF] Sun | A Quantitative Analysis of Current Practices in Optical Flow Estimation and the Principles Behind Them[END_REF], we use a continuation method called Graduated Non Convexity (GNC) [START_REF] Blake | Visual Reconstruction[END_REF] to minimize the energy function. The general idea behind it is that the solution of a smoothed version of Algorithm 1 Computation of optical flow on a pyramid level 1: for i = 1 to the max number of warping steps do 2:

Compute D t,u and D t,v as in equation (4.6) using current u and v

3:

Warp the second image towards the first using current u and v.

4:

Initialize du , dv to zero 5:

for j = 1 to max number of linearization steps do 6:

Linearize the problem using equation (4.1)

7:

Compute the weight Ψ .

8:

Solve for du and dv using equation (4.7)

9:

end for 10:

update u and v using du and dv 11: end for the problem represents a good initialization to solve the original problem. In practice, since we are dealing with a complicated optimization problem, the problem is firstly convexified using the quadratic formulation and then it is linearly and gradually changed until the original problem formulation. The energy function to be minimized is obtained using:

cE Q + (1 -c) E R (4.8)
where c [0, 1], E Q is the convex quadratic approximation of the problem, E R is the robust and possibly non convex original problem. The whole process is depicted in

Evaluation

One of the key factors of the maturity of optical flow approaches ensues from the availability of ground truth data. Authors compete in benchmarks such as Middlebury (Baker et al., 2011) and MP Sintel [START_REF] Butler | A Naturalistic Open Source Movie for Optical Flow Evaluation[END_REF] to evaluate their models. When optical flow was borrowed by experimental fluid dynamics community, similar datasets were required while considering the governing physics equations. A reference simulation based on Navier-Stokes equations is used to generate such images using the simulated (and hence, known) velocity fields [START_REF] Carlier | Report 1 on Production and Diffusion of Fluid Mechanics Images and Data[END_REF]. For river sequences however, a 3D simulation that respects the physics of rivers and then generate realistic 2D image sequences from it is a tedious job to do. Instead a simple idea has been utilized in this work to generate ground truth data. The idea is to track particles on Lagrangian basis throughout the sequence. First, a map of reference trajectories is obtained by tracking these particles manually. The manual tracking is performed at pixel-level (after zooming) using CellTracker free software [START_REF] Piccinini | CellTracker (not only) for dummies[END_REF]. The error related to manual selection depends on the quality of the images. The error in the sequences we used in this study is about 1-2 pixel maximum in the transverse direction. On the longitudinal direction, the error is 1-2 pixels maximum for high quality images and 2-3 for degraded ones. Generally, we select the brightest pixel within a group of pixels that forms the particle. Figure (5.1) shows the location of a chosen particle in the original image and the same particle after zooming in two different images within the sequence. Then, integrating in time (through a 4th order Runge-Kutta), we reconstruct trajectories starting from the same initial points from the estimated displacement fields. The reconstructed trajectories are then compared to their reference trajectories. We test the approach on different image sequences. The sequences differ mainly in their displacement magnitude, image quality, density of tracers (if any) and surface perturbations. The reference trajectories were reconstructed up to few pixels uncertainty. This is because particles change their orientations while moving and their form change as well when they are partially (or sometimes fully) occluded by water. We compare our method SGSD to a modified [START_REF] Horn | Determining Optical Flow[END_REF] variation (HS). In fact, (HS) method used here is exactly the same as SGSD but does not include the turbulent diffusion term. A PIV-based method [START_REF] Ray | Computation of Fluid and Particle Motion From a Time-Sequenced Image Pair: A Global Outlier Identification Approach[END_REF] that is capable of producing dense vector fields (necessary for trajectory reconstruction) is also compared to the two previous methods.

We use the same parameters for SGSD and for (HS) to highlight the difference of results with and without the diffusion term. Some of the aforementioned parameters are general and related to optical flow itself. For instance we use only two-stages GNC iterations in general. Additional GNC iterations improved the results for sequences with good seeding. Median filtering is another parameter to be considered. In general it is well-suited for sequences with less turbulence. If applied to turbulent sequences it might reject turbulent (but maybe correct) vectors and considers them as outliers. However, we found using it only between GNC iterations (and not between image pyramids) to be beneficial in all cases. The weighting factor between the data term and the regularization term is chosen as 0.6 for all sequences. We found that in general, the Lorentzian robust function works best for "good sequences". For noisy sequences or for those for which the 2D divergence assumption does not hold, the "Charbonnier" robust function [START_REF] Charbonnier | Two Deterministic Half-Quadratic Regularization Algorithms for Computed Imaging[END_REF] gives slightly better results than the Lorentzian. We compute the normalized distance error of trajectories of different methods and also for SGSD with different Sc t values:

Error(i) = (y(i) ref -y(i)) 2 + (x(i) ref -x(i)) 2 (y(1) ref -y(end) ref ) 2 + (x(1) ref -x(end) ref ) 2
where ref subscript refers to a reference trajectory component. We chose to normalize the error because we have different displacement values for different sequences. Note that smaller value for Sc t means larger diffusion coefficient and vice versa. We also exploit the dense nature of the estimation to compute relevant differential operators in a straightforward way.

Results on conventional images

We collect videos of rivers in motion from various sources like YouTube or from research institutes data repositories. We require a fixed camera to make sure the motion is only due to the river. We also require the existence of at least few particles for the purpose of evaluation. The chosen particles are selected in such a way that they stay visible (at least partially) throughout the sequence.

The Arc river first sequence

The Arc river, located in the french Alps, is known for its dark color water which enhances the contrast (the gradient) of the image intensity when coupled with white tracers. In addition, the image quality itself is good with no changes in scene lighting. The sequence has 41 images. The average displacement in the streamwise direction is approximately 11 pixels between two images. In Figure (5.2), we plot the mean SGSD dense vector field (with an offset of few pixels for visualization purposes). This gives a quick spatio-temporal idea on the nature of the river at hand. Next, we plot the trajectories of different methods along the reference trajectory. In Figure top), it is easy to observe that all methods were generally in agreement with the reference trajectory. This is mainly because of the good gradient signal in this sequence as mentioned earlier, but also to the quality of the images. Two zoomed areas around the middle and the end of the trajectory are shown. SGSD is almost identical to the reference trajectory, followed by (HS) and then PIV which appears to underestimate the motion magnitude. The zoomed window at the end of the trajectory confirms the results on the zoomed window in the middle (any false positive match in the middle will eventually result in a visible underestimation or overestimation in the window at the end). The normalized distance error in Figure (5.3, middle) confirms the results. In Figure (5.3,bottom) we plot the normalized distance error of SGSD using different values for Sc t . This did not change the results much, there is a maximum error of 0.015 for the case with no diffusion and around 0.05 for other cases with different diffusion coefficients. Still SGSD with Sc t = 1 gives the best results. In Figure (5.4, left), we take advantage of the dense estimation to compute a fine-detailed divergence field. In the same figure on the right, we could compute trajectories on an image cross-section. Observe that these trajectories do not necessarily start at a particle location like the one in Figure (5.3,top). Nevertheless, they are consistent with each other and specially to the ones that start at a particle location.

The Arc river second sequence

In contrast to the first Arc sequence above, the imaging quality of this sequence is degraded. The tracked particle is very small which made even the manual tracking difficult. The sequence is composed of 55 images and the average displacement between two images is approximately 9 pixels. In Figure (5.5), the mean SGSD vector field is superimposed on one of the images. One could notice non-zero vectors in the non-fluid area in the bottom left part of the image. These are due to the wind moving the grass and not due to estimation errors. In Figure (5.6,top) we can see that SGSD was able to recover a better vector field in terms of magnitude and direction than PIV. With (HS), both vector fields seem to have the same magnitude, SGSD however was able to recover much better direction. In the same figure, we show the normalized distance error of different methods (middle) and for SGSD using different values for Sc t (bottom). A larger diffusion coefficient gives better results. However, after 40 images it starts to diverge completely to a poor result with approximately 0.13 normalized error. The default SGSD with Sc t = 1 is more stable throughout the whole sequence and most of the time better than the other two cases of no diffusion (HS) or with a small diffusion coefficient. Figure (5.7,left) shows the derived divergence field and cross-section trajectories are shown in the same figure on the right side. Contrary to the previous sequence, there is no clear visible tracers but the trajectories are still consistent.

The Gave de Pau river sequence

The Gave de Pau river, located in the french Pyrenees, is more challenging than the previous cases. The sequence has many uniform regions with weak gradient signal. Furthermore, the tracked particle is too big in size which increases the uncertainty of the reference trajectory reconstruction. This is because many pixels in the vicinity of the original selected pixel resemble each other. In addition, the river area where the particle resides (and a great deal of other parts) suffered from severe and persistent 3D motion which causes the particle to move up and down for extended period of time. This clearly violates the 2D incompressibility condition used to derive the optical flow model. The sequence is composed of 49 images with an average displacement in the streamwise direction of approximately 6 pixels. In Figure (5.8), we get to have a general idea of the sequence motion patterns by plotting the mean SGSD vector field. The 3D deformations mentioned earlier are readily noticeable when visualizing the dense vector field. SGSD was able to recover a better trajectory compared to other methods, Figure (5.9,top). Testing different Sc t values other than the default value 1 showed similar performance for all variations until around the 20th image where the default SGSD showed better performance while SGSD with a larger diffusion coefficient gives a trajectory that diverges drastically to a bad local minimum with approximately 0.5 normalized error, Figure (5.9, bottom).

Results on ortho-rectified images

In river velocimetry, LSPIV is the benchmark image-based method to estimate free surface velocity. It applies PIV to a grid of points defined in ortho-rectified images. Image ortho-rectification is a way to transform the image to remove the perspective effect so that a direct relationship between image space and the real world could be established, Figure (5.11). The resulting transformed sequence is used to compute river velocity using image-based techniques like PIV. Unfortunately, the philosophy of LSPIV is different from optical flow which makes the comparison difficult. LSPIV relies on estimations obtained on many images and use statistics (the mean mostly) to average the final flow field to reduce the effects of outliers. This makes trajectory reconstruction difficult because it uses successive individual estimations and only one outlier is enough to cause the whole trajectory to diverge. Optical flow is different in this regard as it integrates outlier rejection in the process (i.e. even between only two images) via median filtering and/or regularization. In addition, the interpolation step needed for trajectory reconstruction will be more accurate on optical flow since the interpolation is applied between image pixels themselves while in LSPIV one needs to interpolate between grid points that are many pixels apart, Figure(5.11,bottom). We used FUDAA-LSPIV software [START_REF] Le Coz | Imagebased Velocity and Discharge Measurements in Field and Laboratory River Engineering Studies Using the Free FUDAA-LSPIV Software[END_REF] to ortho-rectify the sequence and compute real world velocity on the defined grid. The sequence is extracted with 5 images per second frame-rate. We run SGSD on the ortho-rectified sequence, Figure (5.12) shows SGSD trajectory reconstruction in comparison to the manual reference trajectory. To compare SGSD to LSPIV, we superimpose the same LSPIV grid on SGSD dense field and compare velocity values only for grid points, Figure (5.13). Both methods give similar results. A noticeable difference is observed on the top and bottom sides of the grid. However, this is not surprising since estimations on image boundaries are known to be prone to errors. In addition, lower and upper part of images have no much tracers to rely on and they contain a fixed shadow pattern.

Conclusion

We presented a novel motion estimation technique for the application of imagebased river velocimetry, based on optical flow scheme. We showed that optical flow in its own right is optimal and physically compatible with rivers free surface motions. In the ideal settings as in the first Arc sequence, section (5.1.1), original optical flow provided very good results without any additional terms. Our enhanced model is based on the decomposition of the scalar transport equation which is equivalent to the original optical flow in addition to a new term. This term, considered to be related to turbulent viscosity, models the small scale contributions that are dropped out in the image acquisition phase. We suggested a turbulent diffusion coefficient based on Prandtl Mixing Length model and a Sc t value of 1. This choice of Sc t gave more stable results. Other values (bigger or smaller than 1) either gave less accurate results or good results at the beginning of the sequence but then drastically diverge with time. We presented results on different possible cases of rivers with or without tracers, good or poor imaging conditions, high or low surface velocity and high or low surface perturbations. Our method outperformed all other approaches in all image sequences in the recovered motion magnitude and/or direction. Both assessed visually or statistically via trajectory reconstruction of particles of interest. We draw the attention here to the sensibility of trajectory reconstruction to outliers, any wrong value at any instance of the estimated vector fields sequence will ruin the part of the trajectory that follows it. It is also sensitive to systematic errors that would accumulate in time to reconstruct erroneous trajectories. Nevertheless, we have been able to reconstruct trajectories that are similar to their reference trajectories. It is also shown how the dense estimation of optical flow facilitates the computation of other relevant quantities like divergence and vorticity fields. Since the application treated here mainly focuses on river velocimetry, our model is designed to penalize the vorticity of the vector field. However, the variational formulation is flexible and could open the door to more specialized optical flow models for complex river flows. In this case, new motion patterns will emerge that are different from the translatory pattern assumed in this paper. For example, a river flow around an object will introduce large-scale rotational motion pattern. The optical flow model could be easily changed to avoid the penalization of this quantity during the estimation process by using the regularization suggested by [START_REF] Chen | Optical Flow for Incompressible Turbulence Motion Estimation[END_REF]. Using trajectory reconstruction evaluation method feedback, new data terms/regularizers and models for the diffusion coefficients could also be tested and assessed.

Figure 1

 1 Figure 1.1: 3D motion projection to images.

Figure

Figure 4

 4 Figure 4.1: A flowchart describing SGSD estimation cycle.

Figure 5 . 1 :

 51 Figure 5.1: Manual tracking of particles: The position of the particle in the first image (top). Brightest position is chosen on pixel level after zooming in two different images (bottom).

Figure 5 . 2 :

 52 Figure 5.2: Mean dense SGSD velocity field of Arc river first sequence

Figure 5 . 5 :

 55 Figure 5.5: Mean dense SGSD velocity field of Arc river second sequence

  Figure (5.10, left) shows the divergence field and the cross-section trajectories are plotted on the right in the same figure. One could see the effect of the aforementioned 3D motion on the highlighted rectangular area plotted on the trajectories.

Figure 5 . 9 :

 59 Figure 5.9: The trajectories of different methods superimposed on the first image of Gave de Pau river sequence (top), normalized distance error of different methods (middle) and the normalized distance error of SGSD with different turbulent Schmidt number (bottom)

Figure 5 .

 5 Figure 5.11: Original image (top), ortho-rectified image with grid points in red (bottom)

Acknowledgements

This work was supported by EDF-DTG and SCHAPI. Authors also would like to thank Jérôme Le Coz and Alexandre Hauet for their insightful comments on FUDAA-LSPIV software along with the LSPIV project used in section 5.2. Adrian, R., Jan. 1991. Particle-Imaging Techniques for Experimental Fluid-Mechanics. Annual Review of Fluid Mechanics 23 (1), 261-304. Baker, S., Scharstein, D., Lewis, J. P., Roth, S., Black, M. J., Szeliski, R., Mar. 2011. A Database and Evaluation Methodology for Optical Flow. International Journal of Computer Vision 92 (1), 1-31.


Francisco Chinesta
email: francisco.chinesta@ec-nantes.fr

Adrien Leygue

Felipe Bordeu

Amine Ammar

et al.. From PGD Elías Cueto

Antonio Huerta

Roland Keunings

From PGD to iPGD

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

ICTEAM, Université Catholique de Louvain Bat. Euler, Av. Georges Lemaitre 4, B-1348 Louvain-la-Neuve, Belgium Today many problems in science and engineering remain intractable, in spite of the impressive progresses attained in mechanical modelling, numerical analysis, discretization techniques and computer science during the last decade, because their numerical complexity is simply unimaginable. We can distinguish different challenging scenarios for efficient numerical simulations:

• The first one concerns models that are defined in high dimensional spaces, usually encountered in quantum chemistry and kinetic theory descriptions of complex fluids. • The second problem category involves time-dependent problems whose spectrum of characteristic times is so large that standard incremental techniques fail. • Real time simulators are needed in many applications. Control, malfunctioning identification and system reconfiguration also need to run in real time. • In problems defined in degenerate geometrical domains (beams, plates, shells …) standard grid-based discretization methods fail because the extremely fine meshes. • Many problems in process control, parametric modeling, inverse identification, and process or shape optimization require the solution of many direct problems. • A new paradigm in the field of Applied Sciences and Engineering has emerged in the last decade. Dynamic Data-Driven Application Systems (DDDAS) allow the linkage of simulation tools with measurement devices for real-time control of simulations. • Augmented reality requires efficient (fast and accurate) simulations. The idea is supplying in real time appropriate information to the reality perceived by the user. • Light computing platforms (e.g. smartphones) are appealing alternatives to heavy computing platforms.

Our main objective is to address the modeling and simulation of "real" models encountered in science and engineering with all their complexity from the geometrical and constitutive points of view, some of them never until now solved because their computational complexity. These models should be solved very fast, in some cases in real time, by using light computing platforms. Classical simulation techniques fail to fulfill the above requirements. An appealing alternative consist of considering off-line solutions of parametric models, in which all the sources of variability -loads, boundary conditions, material parameters, geometrical parameters, etc. -will be considered as extra-coordinates. Thus, by solving only once the resulting multidimensional model, we have access to the solution of the model for any value of the parameters considered as extra-coordinates. Now, from this general solution computed only once and off-line we could perform on-line real time post-processing, optimization, inverse analysis, analysis of sensibilities, stochastic analysis … by using very light computing platforms as for example smartphones. We could also adapt the model on-line while its simulation is running within the framework of dynamic data-driven application systems -DDDAS. The price to be paid is the solution of parametric models defined in high dimensional spaces that could involve hundreds of coordinates. The use of the Proper Generalized Decomposition that we recently proposed and we are intensively developing, allows such solution, because thanks to the separated representation of the unknown fields the computational complexity scales linearly with the dimensionality, instead of growing exponentially which is characteristic of mesh-based discretization techniques [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids. Part II: Transient simulation using space-time separated representations[END_REF][START_REF] Ammar | Non-incremental strategies based on separated representations: Applications in computational rheology[END_REF][START_REF] Chinesta | Recent advances in the use of the Proper Generalized Decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | An Overview of the Proper Generalized Decomposition with Applications in Computational Rheology[END_REF][START_REF] Ladeveze | Nonlinear computational structural mechanics[END_REF][START_REF] Ladeveze | The LATIN multiscale computational method and the proper generalized decomposition[END_REF]. This off-line-on-line Proper Generalized Decomposition Based Dynamic Data-Driven Application Systems could constitute a new paradigm in computational sciences.