Gabrielle Thiébaut
email: gabrielle.thiebaut@univ-rennes1.fr

Hector Rodriguez -Perez

Reciprocal interactions between the native Mentha aquatica and the invasive Ludwigia hexapetala in an outdoor experiment

Keywords: plant performance, functional traits, plant ratio, positive interactions, competition

A high ratio of M. aquatica stimulated L. hexapetala's flowering -The total biomass and branching of L. hexapetala

Introduction

Non-native species must overcome numerous geographical and biological barriers to become invasive in their introduced range (Richardson et al., 2000;[START_REF] Blackburn | A proposed unified framework for biological invasions[END_REF]. The establishment success of invasive species relies on the assumption of niche differences to native species, involving either the exploitation of unused resources (empty niche) (Elton, 1958) or competition mechanisms to access a shared resource (niche replacement). Species sharing the same resources commonly coexist in nature. However, competition can also result in the extinction of native species (competitive exclusion).

There is a limit to niche overlap or similarity in resource use between native and nonnative species (i.e., the theory of limiting similarity) [START_REF] Macarthur | The limiting similarity, convergence and divergence of coexisting species[END_REF].

When niches overlap, the species with superior competitive ability prevails (MacDougall et al., 2009). Differences among species in how they use resources determine the outcome of this interaction.

Numerous experiments in aquatic ecosystems have repeatedly shown competitive advantages of invasive aquatic plants over native ones (Spencer and Rejmanek 2010; [START_REF] Martin | Competition between two aquatic macrophytes, Lagarosiphon major (Ridley) Moss (Hydrocharitaceae) and Myriophyllum spicatum Linnaeus (Haloragaceae) as influenced by substrate sediment and nutrients[END_REF], and even over other non-native plants ( [START_REF] Barrat-Segretain | Experiments on growth interactions between two invasive macrophytes species[END_REF][START_REF] Mony | Competition between two invasive Hydrocharitaceae (Hydrilla verticillata (L.f.) (Royle) and Egeria densa (Planch)) as influenced by sediment fertility and season[END_REF]Gerard et al., 2014). The competitive advantage of invading species is often attributed to traits related to physiology and morphology, such as leaf-area allocation, shoot allocation, growth rate, size, and fitness (van Kleunen et al., 2010;Cuda et al., 2015). However, invasive species are not always superior to native species (Daehler, 2003). The biotic resistance hypothesis (Elton, 1958) partly explains the failures of many invasion attempts. According to this hypothesis, native species can repel invasive plants through competitive exclusion (MacDougall et al., 2009). Therefore, the outcome of competition among non-native and native species depends both on environmental conditions (Daehler, 2003), such as sediment fertility and nutrient availability [START_REF] Mony | Competition between two invasive Hydrocharitaceae (Hydrilla verticillata (L.f.) (Royle) and Egeria densa (Planch)) as influenced by sediment fertility and season[END_REF]Stiers et al., 2011;Gérard et al., 2014;[START_REF] Martin | Competition between two aquatic macrophytes, Lagarosiphon major (Ridley) Moss (Hydrocharitaceae) and Myriophyllum spicatum Linnaeus (Haloragaceae) as influenced by substrate sediment and nutrients[END_REF] and on the functional group identity of the native species (Petruzzella et al., 2018).

These complex interactions may support native species' persistence and coexistence with invaders (Gurevitch and Padilla, 2004;Sax and Gaines, 2008), yet studies focusing on the reciprocal impact of native species on invasive species are rare [START_REF] Leger | Coevolution between native and invasive plant competitors: implications for invasive species management[END_REF].

This study examined a non-native versus a native species, Ludwigia hexapetala (Hook. and Arn.) Zardini, H. Y. Gu and P. H. Raven (syn. L. grandiflora subsp. hexapetala), and Mentha aquatica L. respectively, as biological models. Water primrose L. hexapetala is an invasive species in Europe, native to South America (Thouvenot et al., 2013a). Water primrose was introduced voluntarily into south-eastern France in approximately 1820 (Thouvenot et al., 2013a). It is a perennial aquatic plant forming dense mats. It grows horizontally on water (or mud) and can emerge above the water surface. Early growth consists of rosette-like clusters of rounded leaves on the water's surface. It is mainly aquatic but is also able to colonise terrestrial habitats, such as riverbanks or wet meadows (Thouvenot et al., 2013a). This plant is able to outcompete both submerged and emergent native species (Dandelot et al., 2005;Stiers et al., 2011).

Ludwigia hexapetala produces allelochemicals (Dandelot et al., 2008;Santonja et al., 2018), which could be implicated in the outcome of the interactions between water primrose and native species. In France, L. hexapetala has not been imported since 2007 because its sale and introduction in natural areas has been forbidden by law. It was included in the first European Union list of 37 invasive species, adopted by the European Commission in July 2016. According to EU criteria, species included in the list can cause damage on a scale that justifies dedicated measures applicable across the European Union.

A C C E P T E D

M A N U S C R I P T
These measures effectively prevent, minimise, or mitigate the adverse impact of these species in a cost-effective manner.

Mentha aquatica (Lamiaceae) is a native macrophyte from the northern temperate regions of Europe. It has a creeping rhizome with submerged leaves, and the erect stems possess aerial leaves. Mentha aquatica is most common along wet or waterlogged riverbanks with the shoots exposed to the air, but it can also grow fully submerged during wet seasons. According to the biotic resistance hypothesis, M. aquatica could repel L. hexapetala through competitive exclusion (MacDougall et al., 2009). The invasive L. hexapetala and the native M. aquatica may co-occur in the wild in European aquatic ecosystems.

The aim of this study was to investigate, experimentally, in a common garden, the effect of neighbour (M. aquatica versus L. hexapetala) on the performance (biomass, production of flowers and life-cycle completion) of two species at a constant density. It was hypothesised that the performance of L. hexapetala and of M. aquatica would depend on the ratio and/or on the identity of the neighbouring species.

Materials and methods

Two macrophyte species were selected for the study: the non-native Ludwigia hexapetala and the native Mentha aquatica.

Experimental design

In mid-March 2017, 150 shoots each of L. hexapetala (Lh) and M. aquatica (Ma)

were collected from a channelized stream (Marais de Mainguy, 47°07'21.7"N 2°00'14.9"W) and from a pond (Pornic, 47°07'22.5"N 2°05'22.0"W), respectively, in western France. At these two sampling sites, the two species did not co-occur. Shoots (hereafter called individuals) derived from a single stand either of L. hexapetala (Lh) or

A C C E P T E D M A N U S C R I P T temperature.
For each individual, an apical shoot, without buds or lateral stems, was cut to 7 cm in length. The mean initial fresh weights of L. hexapetala and of M. aquatica were, respectively, 0.6 ± 0.24 g FW and 1.59 ± 0.74 g FW. In the laboratory, the apical shoots (150 per species) were gently rinsed with tap water to remove invertebrates, algae, and debris in preparation for planting. Each shoot was planted in a pot (22 cm in diameter and 30 cm in height), containing a fertile agricultural soil (NPK= 14:10:18 kg/m 3 , pH=6). The effects of the experimental conditions were tested by combining different ratios of M. aquatica and L. hexapetala individuals for a total density of six individuals per pot, following a replacement series design [START_REF] Cousens | Aspects of the design and interpretation of competition (interference) experiments[END_REF]. In the standard replacement design, the same density is used for both species in their respective monocultures; mixtures are formed by substituting an equal number of plants of one species for that of the other, so that the total number of plants per unit area is constant. The bottom of the pots were kept in tap water (ca. 10 cm depth). The water chemistry was basic with moderate nutrient concentrations (mean annual value according to French Water Agency data: conductivity = 400 µS cm -1 ; pH= 8.10; [NO3 -N] = 1.44 mg L -1 ; [NH4 + N] = 0.03 mg L -1 ; [PO4 3-P] = 0.05 mg L -1 ).

Measurement of morphological traits

The number of flowers and the length of the main shoot (of all individuals) were measured in all of the pots each week for 11 weeks. The Relative Growth Rate (RGR; d -1 ) was also calculated on a weekly basis, as suggested by [START_REF] Hunt | Basic growth analysis[END_REF]:

RGR stem = (ln L2 -ln L1)/(T2 -T1)
where L1 and L2 represent total length at time T1 (beginning of the experiment) and T2 (end of the experiment).

At the end of the experiment, the main shoot length was measured, the number of lateral branches was counted, and the roots and shoots were harvested. The above-ground and below-ground parts of the plants were dried separately at 65°C for 72 hours and weighed. The below-ground/above-ground ratio was calculated. RGR and stem length are indicators of apical growth, whereas the number of lateral branches is an indicator of lateral growth and the plant's ability to regenerate [START_REF] Barrat-Segretain | Comparative abilities of vegetative regeneration among aquatic plants growing in disturbed habitats[END_REF]. The number of roots is an indicator of the plant's ability to colonise [START_REF] Barrat-Segretain | Comparative abilities of vegetative regeneration among aquatic plants growing in disturbed habitats[END_REF].

Indices of plant competition

Indices were used to quantify and express several attributes of plant competition, including competition intensity and importance, competitive effects and responses, and

Data analyses

RGRs were analysed on a repeated-measures basis by means of non-parametric testing [START_REF] Naguchi | nparLD: An R software for the nonparametric analysis of longitudinal data in factorial experiments[END_REF], since data did not meet criteria for parametric analyses (i.e., homoscedasticity and normality of residuals). Whenever the interaction between treatment (planting ratio) and time was significant, pairwise comparison among treatments was performed within a given sampling time using the Mann-Whitney-Wilcox test. Data pairwise comparisons were only run for the periods when growth rates peaked and a subsequent Benjamini-Yekutieli correction for dependent multiple tests was applied to the multiple test series [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF].

Morphological traits matrices were analysed independently for each species by means of Metric Multidimensional Scaling (MDS) based on Euclidean distance with standardised values. Since the graphical results suggested treatment effects on plant morphology, univariate analyses were used to explore treatment effects on the below/above-ground biomass ratio, total biomass, main stem length, number of flowers (total number of flowers produced throughout the whole experiment), and the number of

A C C E P T E D M A N U S C R I P T
branches. Whenever possible, a linear model approach was used for univariate analysis.

A type II sum of squares was used in each case for testing the significance of the planting ratio treatment effect, owing to the unbalanced design [START_REF] Langsrud | ANOVA for unbalanced data: use Type II instead of Type III sums of squares[END_REF]. A linear model approach was used to analyse treatment effect on the harvested biomass. In the case of the M. aquatica data set, biomass was log-transformed prior to the analysis in order to meet normality and homoscedasticity requirements. Residuals were explored graphically to assess departure from the assumptions of the analysis [START_REF] Bolker | Generalized linear models : a practical guide for ecology and evolution[END_REF]. Number of branches and number of flowers were analysed with a general linear model (GLM) using a Poisson distribution. Patterns in residuals were explored graphically to assess departure from the assumptions of the analysis [START_REF] Bolker | Generalized linear models : a practical guide for ecology and evolution[END_REF]. As the number of flowers for M. aquatica did not meet normality and homoscedasticity requirements, a Kruskal-Wallis test was run. Whenever the treatment effect of the planting ratio was significant, HSD Tukey pairwise comparisons were run to determine significant differences among the levels of the experimental factor. Pairwise comparisons based on a non-parametric Wilcox test and subsequent Bonferroni corrections were used for the different ratios when a significant p-value was obtained with the Kruskal-Wallis univariate test. All analyses were performed with R software (R Core Team, 2016). The nparLD package [START_REF] Naguchi | nparLD: An R software for the nonparametric analysis of longitudinal data in factorial experiments[END_REF] was used for nonparametric analysis of longitudinal data; the multcomp package for HSD Tukey pairwise comparisons [START_REF] Hothorn | Simultaneous inference in general parametric models[END_REF]; and the vegan package (Oksanen et al., 2017) for MDS analysis.

Results

The calculation of the RYTab index indicated that the two species did not share any common limiting resources, suggesting there was no competition between M. aquatica and L. hexapetala (Fig 1A). The highest biomass total yield (RYT) was observed in pots at the planting ratio 1Ma/5Lh and, to a lesser extent, in pots at the planting ratio

A C C E P T E D M A N U S C R I P T 5Ma/1Lh (Fig 1A).
We observed an increase in biomass production RY of L. hexapetala RY of M. aquatica (over-yielding) for both species in mixtures at the planting ratio 1Ma/5Lh, and an increase in biomass production (RY) of L. hexapetala only at the planting ratio 5Ma/1Lh (Fig. 1A).

The RCI standardises competition intensity, allowing the comparison of species with different RGR in the absence of interaction. Negative RCI values were found for both species at the ratio of 1Ma/5Lh, whereas positive values were calculated at the ratio 3Ma/3Lh (Fig. 1B). The RCI was close to 0 for M. aquatica and negative for L.

hexapetala at the ratio of 5Ma/1Lh.

The highest RGRs for both L. hexapetala and M. aquatica were documented during the third growth period (week 3). For L. hexapetala, it was significantly lower at the ratio of 5Ma/1Lh than when the individual grew in a monoculture or at the 1Ma/5Lh ratio (Table 1; Fig. 2A). On the other hand, M. aquatica exhibited highest growth rates in the 1Ma/5Lh and 3Ma/3Lh ratio treatments for the same periods (Fig. 2B).

The overlap of the treatments was greater for M. aquatica than for L. hexapetala (Fig. 3). Taking into account centroids in MDS analysis (Fig. 3A), the traits of L.

hexapetala shoots issued from 5Ma/1Lh and 1Ma/5Lh treatments differed from those coming from the monoculture treatment (0Ma/6Lh). The shoots of M. aquatica coming from the 1Ma/5Lh treatment were visually distinct from those coming from other treatments (Fig. 3B). It is worth noting that the greatest variability was consistently observed in the 5Ma/1Lh treatment for L. hexapetala and the 1Ma/5Lh treatment for M.

aquatica.

The total biomass of Ludwigia hexapetala individuals was significantly smaller in monoculture than in mixture (1Ma/5Lh and 5Ma/1Lh treatments). The individuals of L.

hexapetala had significantly fewer branches in monoculture than in the 1Ma/5Lh

A C C E P T E D M A N U S C R I P T
treatment (Table 2; Fig. 4). On the other hand, L. hexapetala individuals produced significantly more flowers in 5Ma/1Lh than in 1Ma/5Lh and in monoculture (Table 2; Fig. 4). At the ratio of 1Ma/5Lh, M. aquatica individuals had the lowest below/aboveground biomass ratio, by a significant margin (Table 2; Fig. 4).

Discussion

Impacts of Mentha aquatica on the performance of Ludwigia hexapetala

Although most studies on the interactions among invasive and native aquatic plants have focused on competition (Spencer and Rejmanek, 2010; [START_REF] Martin | Competition between two aquatic macrophytes, Lagarosiphon major (Ridley) Moss (Hydrocharitaceae) and Myriophyllum spicatum Linnaeus (Haloragaceae) as influenced by substrate sediment and nutrients[END_REF]Gérard et al., 2014), there is increasing recognition that facilitation can be important in such interactions (Rodriguez, 2006). 'Facilitation' is defined as an interaction by which one species positively impacts the performance of another. The RCI indices suggested a facilitation effect for L. hexapetala in the presence of a low or high ratio of M. aquatica.

The number of lateral branches and the total biomass of L. hexapetala were stimulated by the presence of one individual of M. aquatica. Moreover, the study findings also showed a positive effect of the presence of M. aquatica on the total biomass and on the production of flowers by L. hexapetala when water primrose was out numbered by M. aquatica. It is plausible that these positive effects were mediated by some change in the soil environment such as inhibition of denitrification to enhance release of nutrients that also benefit the water primrose (indirect facilitation). The allocation of energy to lateral growth, biomass, and flowering suggests an efficient strategy for avoiding competition with M. aquatica.

Surprisingly, no negative impact of M. aquatica on the invasive L. hexapetala was found.

A possible explanation is the low level of competition for resources in the mesocosms, as indicated by the value of the RYTab. During the experiment, L. hexapetala adopted a creeping mode of growth, whereas M. aquatica produced erect stems. These two types of The number of meristems that can potentially produce reproductive organs can be inhibited by competition for light [START_REF] Bonser | Allometry and development in herbaceous plants: functional responses of meristem allocation to light and nutrient availability[END_REF], negatively influencing the number of flowers. The stem lengths of L. hexapetala individuals in monocultures were generally small. Smaller individuals with low reproductive allocation may simply have been growing and developing more slowly, and might have reached the same biomass if they had had more time to develop. Plants were also measured after three weeks, when a lower relative growth rate of M. aquatica in monoculture and at the high planting ratio (5Ma/1Lh) suggested a competition for resources within the species. This hypothesis is also supported by competition indices (RCI) results.

Impacts of Ludwigia hexapetala on the performance of Mentha aquatica

The results did not show a negative impact of L.hexapetala on the growth of M. aquatica, whereas in the field the presence of L. hexapetala affects the abundance of native species such as Alisma plantago-aquatica, Ceratophyllum demersum, and Lycopus europaeus (Stiers et al., 2011). However, this study established that the below/aboveground biomass ratio of M. aquatica was reduced in the presence of a high planting ratio of L. hexapetala (1Ma/5Lh). In a previous study, the authors showed that M. aquatica root length decreased at a higher L. hexapetala ratio, and no effect was found of the

A C C E P T E D M A N U S C R I P T
presence of L. hexapetala on the relative growth rate of M. aquatica (Thouvenot et al., 2013b). These results are congruent with Vila and Weiner's ( 2004) review of pairwise experiments on non-native and native plant species, which supports the general perception that the effect of invasive species on native species is usually stronger, and not vice versa. The lower root biomass of M. aquatica observed in the presence of five individuals of L. hexapetala could be due to the synthesis of allelochemicals by the water primrose roots (pedunculagin and an ellagic acid) (Marcellin-Gros, 2015) and a subsequent inhibition of M. aquatica's root growth. Indeed, ellagic acid is a rooting inhibitor (Qin et al., 2006). However, in our study, we did not specifically test this hypothesis. Surprisingly, the presence of five individuals of L. hexapetala could benefit the growth rate of M. aquatica. However, this result was only measured once, during the most active growing period (after three weeks). At this phenological stage, there was no competition for light between water primrose and M. aquatica.

The absence of a significant effect of different ratios on the traits of the native species related to growth (stem length, number of lateral branches, total biomass) could be explained by (a) similar competitive abilities of the two species; (b) high variability of competitive abilities within species; or (c) because the intensity of competition was not sufficiently high to cause significant differences in species performance within the time scale of the study. This study has shown that a native plant species can coexist with an invasive plant species at a certain ratio, suggesting the possibility of adaptive evolutionary responses in this plant to the invasive plant. Indeed, it has been hypothesised that strong competition from invasive plants eliminates native plant genotypes that cannot resist or tolerate strong competition, resulting in an accumulation of native plant genotypes that can resist or tolerate such competition within the native populations [START_REF] Leger | Coevolution between native and invasive plant competitors: implications for invasive species management[END_REF]. The greater fitness of these native populations could be a result of their evolved

A C C E P T E D M A N U S C R I P T
tolerance or competitive ability against non-native species [START_REF] Leger | Coevolution between native and invasive plant competitors: implications for invasive species management[END_REF][START_REF] Oduor | Evolutionary responses of native plant species to invasive plants: a review[END_REF]. In the field, the two species can co-exist, they can develop a common interaction history, allowing the native species to partially resist the invasive L.

hexapetala.

The results of this study suggest that the outcome of the interactions (Competition, Facilitation) between L. hexapetala and M. aquatica was partially determined by the planting ratio of the species. The study has shown strong intraspecific competition within L. hexapetala individuals and facilitation interactions between L. hexapetala and M.

aquatica, indicating that species coexistence is favoured over competitive exclusion. The findings can be extrapolated to the broader context of how native and non-native species can exist together (or not), although this study investigated only two species, one example of a native and non-native species interacting.

Thus, the strength of competitive interactions between invasive and native plants that might result in evolutionary responses from the native species may be contextdependent, varying with both habitat productivity and the inherent traits of the invasive and native plants under consideration (Funk and Vitousek, 2007;[START_REF] Leger | Coevolution between native and invasive plant competitors: implications for invasive species management[END_REF].

  competition (Weigelt and Joliffe, 2003). Two indices were selected: the Relative Yield Total index (RYTab) (de Witt andVan den Bergh, 1965) reflects how plant performance is influenced by population density or neighbour size; the Relative Competition Intensity (RCI) index(Grace, 1995) quantifies the intensity of the competition. The RYTab (de Witt andVan den Bergh, 1965) was calculated for each ratio treatment, using the following equation:RYTab = RYa + RYbWhere RYa and RYb corresponded to the individual relative yields of each species for the given magnitude.The relative yield of a species is:RY = Ymix/Ycontr in which Ymix and Ycontr are the biomasses in mixture and in monoculture. Individual relative yield divides the individual biomass for each plant of the given species in the given treatment by the average yield of the biomass of the same species cultivated in monoculture conditions. According to Snaydon (1991), an RYTab value of 2 would indicate that the species did not share any common limiting resource, i.e. it did not compete and showed full resource complementarity, while RYTab values greater than 1 indicate that the species did not fully share common limiting resources, i.e. it competed partially and showed partial resource complementarity, and RYTab values of 1 indicate that the species fully shared the same limiting resources, i.e. competed fully and showed no resource complementarity. the RCI (Grace 1995), with the following formulation: RCI = (Pmono -Pmix)/Pmono Where Pmono represents the performance of a plant (hereafter total biomass) in a monoculture and Pmix represents the performance (i.e. total biomass) of a plant in a mixture. RCI compares the performance of plants growing in a mixture with control plants growing alone.The RCI score was calculated for each species at each planting ratio in mixtures (i.e. 1Ma/5Lh, 3Ma/3Lh, 5 Ma/1Lh). According toGoldberg et al. (1999), a negative score indicates a better performance of the plants growing in mixture than individuals growing in monoculture, suggesting facilitation interactions between the two species. A positive value indicates a weaker performance of the plants in mixture in comparison with plants coming from a monoculture, suggesting competitive interactions between the two species. There is no minimum RCI value for facilitation, but there is a maximum value of 1 indicating maximal competition. An RCI of 0 indicates null interaction effects.

  versus erect, meant competition for space and light was avoided, and suggests niche partitioning. Competition is expected to occur mostly between species of similar growth form and occupying similar niches. The individuals of L. hexapetala had less biomass in monoculture than at 1Ma/5Lh and 5Ma/1Lh treatments, and produced fewer lateral branches in monoculture than at 1Ma/5Lh treatment and less flowers in monoculture than at 5Ma/1Lh treatment, suggesting strong intraspecific competition for resources between the individuals of L. hexapetala. Competition for light is likely to have led the plants in monoculture to reduce their investment in branching (lateral growth).

  . C., 2003. Performance comparisons of co-occurring native and alien invasive plants: Implications for conservation and restoration. Annu. Rev. Ecol. Evol. Syst. 34, 183-211. Dandelot S., Matheron, R., Le Petit, J., Verlaque, W., Cazaubon, A., 2005. Temporal variations of physicochemical and microbiological parameters in three freshwater ecosystems (south-eastern France) invaded by Ludwigia spp. C.R. Biol. 328, 991-999. Dandelot, S., Robles, C., Pech, N., Cazaubon, A., Verlaque, R., 2008. Allelopathic potential of two invasive alien Ludwigia spp. Aquat. Bot. 88, 311-316 de Wit, C.T., Van den Bergh, J.P., 1965. Competition between herbage plants. NJAS 13 (2), 212-221. Elton, C.S., 1958. The ecology of invasion by animals and plants, Chapman and Hall, London Funk, J.L., Vitousek, P.M., 2007. Resource-use efficiency and plant invasion in lowresource systems. Nature 446, 1079-1081 Gérard J, Brion, N., Triest, L., 2014. Effect of water column phosphorus reduction on competitive outcome and traits of Ludwigia grandiflora and L. peploides, invasive species in Europe. Aquat Invasions 9 (2), 157-166 Goldberg, D.E., Rajaniemi, T., Gurevitch, J., Stewart-Oaten, A., 1999. Empirical approaches to quantifying interaction intensity: competition and facilitation along productivity gradients. Ecology 80, 1118-1131. Grace, J.B., 1995. On the measurement of plant competition intensity. Ecology 76

Figure caption Figure 1 .

 caption1 Figure caption

Figure 3 .

 3 Figure 3. MDS plot based on Euclidean distance and morphological traits after

Figure 4 .

 4 Figure 4. Mean values + standard errors of morphological traits (Total biomass,

  Fig 1

  Fig 3

  Fig 4 Total biomass (g DW)

  Köppen climate classification, the climate type is categorised as a Cfb temperate oceanic climate zone, in which the mean annual temperature is 13°C and the mean annual precipitation is 694 mm.

	U S C R I P T
	N
	per pot (six individuals per pot, i.e. six individuals/0.038m²) corresponds to the density M The overall density A
	typically observed in the field (from 100 to 200 individuals/m²) (G. Thiebaut, unpubl. data). Experimental conditions consisted of monospecific pots (6Ma/0Lh or 0Ma/6Lh; containing six individuals of one plant species, Ma or Lh, respectively), mid-ratio pots (3Ma/3Lh; including three individuals of each species), high ratio (5Ma/1Lh; including five individuals of the focal species Ma, and one single individual of the alternative species Lh), and a low ratio treatment (1Ma/5Lh; including one single individual of the C C E P T E D
	focal species Ma, and five individuals of the alternative species Lh). Each treatment had A
	ten replicates. The experiment was conducted outdoors in the experimental garden
	ECOBIO Research Facility at the University of Rennes 1 (48° 7'4.50"N; 1°38'22.96"W)
	from April (14 th week of the year) to July (26 th week of the year). According to the

  Oksanen, J.,. Blanchet, F.G., Kindt, R., Legendre, P., Minchin, P.R., O'Hara, R.B., Simpson, G.L., Solymos, P., Stevens, M.H.H., Wagner H., 2017. vegan: Community Spencer, D.F., Rejmanek, M., 2010. Competition between two submersed aquatic macrophytes, Potamogeton pectinatus and Potamogeton gramineus, across a light gradient. Aquat. Bot 92, 239-244.

	Ecology Package. R package version 2.4-3.
	Petruzzella A., Manschot J., van Leeuwen C.H.A., Grutters B.M.C., Bakker E.S., 2018. Stiers, I., Njambuya, J., Triest, L., 2011. Competitive abilities of invasive Lagarosiphon
	Mechanisms of Invasion Resistance of Aquatic Plant Communities. Front. Plant Sci. 9, major and native Ceratophyllum demersum in monocultures and mixed cultures in
	1-11. Qin, B., Perry, L.G., Broeckling, C.D., Jiang, D.J., Stermitz, F.R., Paschke, M.W.,Vivanco, J.M., 2006. Phytotoxic allelochemicals from roots and root exudates of leafy spurge (Euphorbia esula). Plant Signal Behav. 1(6), 323-327. R Core Team, 2016. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. Rodriguez, L.F., 2006. Can invasive species facilitate native species? Evidence of how, M Vila, M., Weiner, J., 2004. Are invasive plant species better competitors than native plant species? Evidence from pair-wise experiments. Oikos 105, 229-238. M 6, 93-107. A A 2000. Naturalization and invasion of alien plants:concepts and definitions. Divers. Distrib N N Richardson, D.M., Pysek, P., Rejmanek, M., Barbour, M.G., Panetta, F.D., West, C.J., U S C R I P T relation to experimental sediment dredging. Aquat. Bot. 95, 161-166 Thouvenot, L., Haury, J., Thiébaut, G., 2013a. A success story: water primroses, aquatic plant pests. Aquat Conserv. 23, 790-803. U S C R I P T
	when, and why these impacts occur. Biol Invasions 8, 927-939 Santonja, M., Le Rouzic, B., Thiébaut, G., 2018. Seasonal dependence and functional implications of macrophyte-phytoplankton allelopathic interactions. Freshw Biol. 63 (9), 1161-1172 Sax, D.F., Gaines, S.D., 2008. Species invasions and extinction: the future of native biodiversity on islands. PNAS 105, 11490-11497 C C E P T E D Weigelt, A., Jolliffe, P., 2003. Indices in plant competition. J. Ecol. 91, 707-720 C C E P T E D
	Snaydon, R.W. 1991. Replacement or additive designs for competition studies?. J Appl A A
	Ecol 28(3), 930-946.

Thouvenot, L.,

Puech, C., Martinez, L., Haury, J., Thiébaut,, G. 2013b

. Strategies of the invasive macrophyte Ludwigia grandiflora in its introduced range: competition, facilitation or coexistence with native and exotic species? Aquat. Bot 107, 8-13 van Kleunen, M., Weber, E., Fischer, M., 2010. A meta-analysis of trait differences between invasive and non-invasive plant species. Ecol.

Lett. 13,[235][236][237][238][239][240][241][242][243][244][245]

Table 1 .

 1 Longitudinal data analysis results based on non-parametric repeated measures. RGR is Relative Growth Rate. ATS denotes Anova-type Statistic. Significant p-values at 5% significance level are in bold.

				L. hexapetala			M. aquatica
			ATS	df	p	ATS	df	p
	RGR	Density	2.51	2.54	0.07	1.62	1.64	0.09
		Time	31.28	4.80	<0.001	68.31	6.33	<0.001
		Density x Time	2.49	8.41	0.01	3.05	10.93	<0.001

Table 2 .

 2 Univariate analyses of morphological traits results. The first three variables were analyzed with a linear model and type II sum of squares. Lateral Branches and Flowers were analyzed by means of glm following a Poisson error distribution, except for Number of flowers variable in the case of M. aquatica analyzed with a Kruskal-Wallis test. Significant p-values at 5% significance level are in bold.

		L. hexapetala	M. aquatica
		SS df F	p	SS df F	p
	Total Biomass	28.88 3 4.55 0.004 6.92 3 0.56 0.6
	Stem Length	4.98 3 2.18 0.09 271.4 3 2.27 0.08
	Below/ above ground biomass 0.27 3 2.10 0.1	4.71 3 5.27 0.002
		𝜘 2	df	p	𝜘 2	df	p
	Lateral Branches	13.45	3	0.004	2.77	3	0.4
	Number Flowers	13.06	3	0.005	1.96	3	0.6

1

Acknowledgements

The authors would like to express their gratitude to Quentin Gaudet, Milin Morgane, and Leo Barthelemie for their assistance in setting up the experiment and measuring the traits.

We thank also the staff of the EcolEx platform. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. A C C E P T E D M A N U S C R I P T


domain: sdu sdu.stu info

text:
Introduction

As oil reserves are declining, a strong effort is required in petroleum industry to propose efficient extraction techniques able to increase field production. In this context, developing efficient simulation-based optimization methods for well placement and monitoring is a key challenge. Thanks to intelligent fields or Smart fields (van den [START_REF] Van Den Berg | Business value from intelligent fields[END_REF]; [START_REF] Abdulkarim | Overview of saudi aramco's intelligent field program[END_REF]), a large amount of information can be used to improve our knowledge of the reservoir, in particular to improve our numerical models for generating oil production forecasts. These models, calibrated with history matching (Schulze-Riegert and Ghedan (2007); [START_REF] Oliver | Inverse theory for petroleum reservoir characterization and historymatching[END_REF]; Le Ravalec et al. ( 2012)), have better forecasting capabilities and can be used to make decisions related to the development of oil fields and thus have a more efficient field exploitation.

One of these decisions is the location and geometry (i.e. trajectory) of wells. Optimizing the number of wells, their trajectory, and their type (injector or producer) is a difficult problem and can be very costly in terms of fluid flow reservoir simulations. Optimizing well location has been observed to yield non-smooth cost functions [START_REF] Bouzarkouna | Using evolution strategy with meta-models for well placement optimization[END_REF]; [START_REF] Isebor | Generalized field development optimization using derivative-free procedures[END_REF]). Hence, in optimal well placement one should perform some kind of global search that prevents the optimization from being trapped in solutions that are not acceptable in practice (e.g. "low" Net Present Value (NPV) or total oil production).

In the first section we propose a two-step methodology to solve this challenging problem. The first step consists in solving a Black-Box Mixed Integer Non-Linear Problem (MINLP) with a Pattern Search method [START_REF] Audet | Pattern search algorithms for mixed variable programming[END_REF]; Le Digabel (2011)). This method belongs to the class of heuristic optimization methods that were adapted for Black-Box MINLP problems, but may require a large number of evaluations. The variables of this first optimization problem are the number of wells, their locations and types. We optimize the configuration of vertical wells, this simplification allowing to reduce the number of variables. In order to limit the number of evaluations, the Pattern Search method is coupled with metamodels based upon a Gaussian process approach [START_REF] Gramacy | The mesh adaptive direct search algorithm with treed gaussian process surrogates[END_REF]). Simulation outputs of the obtained solution are used in the second step of the methodology.

In this second step, starting from the optimized vertical well configuration, we solve a simplified problem to optimize trajectories of vertical wells by adding branches. Only two simulations are required in this phase, one to collect reference data from the initial well configuration (optimized result of the first step), and the second one to validate the solution obtained in this second step. The simplifications we propose lead to a convex MINLP problem that can be solved by an exact method as Branch-and-Bound [START_REF] Dakin | A tree search algorithm for mixed-integer programming problems[END_REF]; [START_REF] Gupta | Branch and bound experiments in convex nonlinear programming[END_REF]).

The two phases of the methodology are described in the first section. The second section focuses on the mathematical model built for the simplified MINLP problem of well branching. Finally, we present preliminary computational results of the second step of the methodology obtained for a two-dimensional case.

Methodology

In this section, we describe the proposed methodology: first, the Black-Box MINLP optimization problem for vertical well placement, then, the analysis of the obtained solution and the definition of the input of the second optimization problem for well branching, and lastly, the objectives of this second optimization problem.

Optimization of vertical wells: type and location

In the first optimization, we limit the number of optimized variables by considering a limited number of wells and only vertical trajectories. The optimized objective function is the Net Present Value (NPV) ECMOR XIV -14 th European Conference on the Mathematics of Oil Recovery Catania, Sicily, Italy, 8-11 September 2014 defined as follows

f NPV = N I +N P ∑ n=1 (C F l n log(l n ) +C D ) (1) + T ∑ t=1 1 (1 + R) t N P ∑ n=1 (C PO q PO tn -C PW q PW tn ) - N P +N I ∑ n=N P +1 (C IW q IW tn ) ,
where

• N I and N P are the number of injectors and producers,

• C D is the drilling cost and C F l n log(l n ) accounts for the extra cost related to the length of well n (namely l n ),

• q PO tn , q PW tn and q IW tn are, respectively, produced oil, produced water (for producer well n) and injected water (for injector well n) cumulated quantities during production period t, computed by a fluid flow reservoir simulator,

• R is a discount factor (it implies that early oil production is assigned a stronger weight in the NPV function than the late production) and T is the number of production periods,

• C PO , C PW and C IW are positive weights associated with each term.

Maximizing this objective function with respect to the number of active wells (producers and injectors), their location and their length is a Black-Box MINLP optimization problem coupled with a computationally expensive simulator, namely the fluid flow reservoir simulator. Since the NPV function is generally non-smooth and non-convex and gradient information is rarely available, its minimization requires derivative-free global optimization methods [START_REF] Bouzarkouna | Using evolution strategy with meta-models for well placement optimization[END_REF][START_REF] Echeverría Ciaurri | Application of derivative-free methodologies to generally constrained oil production optimisation problems[END_REF]; [START_REF] Isebor | Generalized field development optimization using derivative-free procedures[END_REF]; [START_REF] Onwunalu | Application of a particle swarm optimization algorithm for determining optimum well location and type[END_REF]). As these methods require a large number of simulations, we propose to couple the optimization method with meta-modeling techniques to replace costly simulations by an approximated response surface. [START_REF] Gramacy | The mesh adaptive direct search algorithm with treed gaussian process surrogates[END_REF] proposed a direct search method coupled with Gaussian process surrogates that permits the introduction of dedicated strategies in the exploration phase, based on surrogate predictions. Gaussian processes provide not only a surrogate model of the objective function to be minimized/maximized, but also an estimate of the predictive variance. This can be useful for selecting the most promising candidates to be simulated in the poll step of the mesh adaptive direct search.

This step provides a first optimized configuration with production and injection vertical wells. The outputs of the fluid flow reservoir simulation obtained with this configuration are analyzed to determine the areas that were not produced.

Analysis of optimized vertical well configuration

In practice, reservoir engineers make the decision to add wells or branches from the analysis of the outputs of the fluid flow simulator. The outputs that are commonly considered are physical attributes, such as oil saturation or oil thickness. We combine this pragmatic approach with optimization to better determine the trajectory of wells with branches. We analyze the outputs of the simulator and attribute a score to distinguish regions of the field: a high score means that an area has good drilling potential. For simplicity, we focused on a single attribute called H u φ S o which is the product of the utility thickness (H u ), the porosity (φ ), and the oil saturation (S o ). H u φ S o quantifies the volume of oil still not produced.

An example of a two-dimensional map, obtained with H u φ S o is displayed in Figure 1. This map is obtained from the outputs of a fluid flow reservoir simulation after producing the field for 4 years with one vertical well at the center of the domain. In this map, we observe that there is still oil not produced.

ECMOR XIV -14 th European Conference on the Mathematics of Oil Recovery Catania, Sicily, Italy, 8-11 September 2014

We also see that some cells have a high H u φ S o but are spatially close to cells with low H u φ S o values.

Clearly, such cells are not interesting targets compared to zones corresponding to a group of cells with more homogeneous and relatively high H u φ S o values. In order to reduce the number of variables, and to avoid isolated valued, scores are computed for groups of cells. Such groups of cells are called areas in the remaining of this paper.

Figure 1 H u φ S o map associated with a single vertical well configuration considered later in the optimization case

Another property that could be taken into account in the definition of high potential areas is the connectivity between cells. Based on connectivity consideration, we can introduce here of forbidden zones, i.e, parts of the reservoir we cannot physically drill or that we do not want to drill because, in general, they are not well connected.

From this analysis, we determine the input data to be used within the following step of the proposed methodology, which is the branching optimization problem.

Optimization of well trajectory: branching problem

Starting from a vertical well configuration, given as the solution of the first optimization problem described above, we want to improve oil production by adding branches to the existing vertical wells. A branch sub-problem is introduced: given a well configuration without branches, we aim at finding the potential branches that can be drilled from each well. We focus on branches from production wells only.

The vertical injection wells obtained in the first optimization solution are fixed and kept unchanged in this optimization.

In order to obtain a very efficient methodology, in this optimization we do not use reservoir fluid flow simulation. For this reason we consider a simplified objective function compared to the NPV function of the first optimization problem. This simplified objective function is an estimation of the oil production based on the H u φ S o values of the cells belonging to all areas connected to the production wells by a branch. An additional constraint is introduced to account for the branch construction cost, which in turn depends on the branch length, to avoid well configuration with an unacceptable drilling cost.

As described in the next section, for modeling purposes we introduce the concept of cluster. A cluster is a set of oil areas that are "close enough" to be reached by a single branch (close enough has to be understood in terms of connectivity). The mathematical model below can be used to aggregate areas into clusters and to determine how to reach a cluster with a branch from a well, with an acceptable construction cost. We define the end of a branch as a cluster reference point, i.e. a point that is close enough to each area of the cluster (the branch thus obtained can be expected to produce a significant amount of the oil in the cluster).

ECMOR XIV -14 th European Conference on the Mathematics of Oil Recovery Catania, Sicily, Italy, 8-11 September 2014

Mathematical model of well branching problem

This section describes the mathematical model designed to solve the branching problem. It is organized as follows: we first present the notations, namely the sets (represented by Greek letters), the input data or parameters (upper capitals), the variables (lower capitals), and the simple bounds on variables. We conclude the section by discussing the objective function and the set of constraints, that are properly reformulated to obtain a convex MINLP problem that can be solved by exact optimization methods as discussed in last section. For the sake of simplicity, we restrict the model to a two-dimensional case (however this formulation can be readily generalized to three-dimensional problems).

Input Data

Sets:

• Φ: set of wells in the initial solution.

• Π: set of areas of the reservoir representing the cells with a high score at the end of the simulation with the initial well configuration.

• K: maximum number of clusters.

• Γ = {1,... ,K}: set of clusters.

Parameters:

• X id : coordinates of the center of oil area i (i ∈ Π, d ∈ {1, 2}).

• O i : volume of oil in area i (i ∈ Π), computed from H u φ S o values in the corresponding area.

• [X d , X d ]: oil reservoir domain for dimension d (d ∈ {1, 2}).

• L, L and H: minimum and maximum lengths of a branch, and maximum length of all branches built. The third parameter is necessary to limit the construction cost based on the target budget.

• R: maximum distance between oil area and branch end, i.e., maximum distance of a point in Π with respect to the cluster reference point.

• V nd : coordinates of the top of well n for dimension d (d ∈ {1, 2}, n ∈ Φ).

• W n : height of well n (n ∈ Φ), we assume that there are vertical wells only.

• B: maximum number of branches per well.

Variables

Continuous:

• z jd : d-th coordinate of the reference point for cluster j ( j ∈ Γ, d ∈ {1, 2}).

• w jn : height of the junction point on well n connecting cluster j through a branch ( j ∈ Γ, n ∈ Φ).

• s jn : square of the distance between the reference point of cluster j and the junction point of well n ( j ∈ Γ, n ∈ Φ).

ECMOR XIV -14 th European Conference on the Mathematics of Oil Recovery Catania, Sicily, Italy, 8-11 September 2014

Binary:

• q jn : activation of connection between cluster j and well n ( j ∈ Γ, n ∈ Φ): [1 if active, 0 otherwise].

• x i j : inclusion of oil area i in cluster j (i ∈ Π, j ∈ Γ): [1 if included, 0 otherwise].

• y j : activation of cluster j ( j ∈ Γ), an inactive cluster is not linked to any well: [1 if active, 0 otherwise].

Some additional parameters and variables will be introduced later to formulate the optimization as a convex MINLP problem.

Simple Bounds

Simple bound on continuous variables are defined to ensure that the well configurations obtained are specified within the physical domain considered.

0 ≤ w jn ≤ W n ∀ j ∈ Γ, n ∈ Φ, (2)
Z d ≤ z jd ≤ Z d ∀ j ∈ Γ, d ∈ {1, 2}, (3) 0 ≤ s jn ≤ ||X d -X d || 2 2 ∀ j ∈ Γ, n ∈ Φ. (4) Constraints ∑ j∈Γ,n∈Φ q jn s jn ≤ H 2 , ( 5
)
q jn L 2 ≤ q jn s jn ∀ j ∈ Γ, ∀n ∈ Φ, (6)
q jn s jn ≤ q jn L 2 ∀ j ∈ Γ, ∀n ∈ Φ, (7)
∑

j∈Γ x i j ≤ 1 ∀i ∈ Π, (8)
x i ≤ y j ∀i ∈ Π, j ∈ Γ, (9)

y j ≤ ∑ i∈Π x i j j ∈ Γ, (10)
q jn ≤ y j ∀n ∈ Φ, j ∈ Γ, (11)
∑ j∈Γ q jn ≤ B ∀n ∈ Φ, (12)
∑ n∈Φ q jn ≤ y j ∀ j ∈ Γ, (13)
∑ d∈{1,2} (X id -z jd ) 2 ≤ R 2 + M(1 -x i j ) ∀ j ∈ Γ, i ∈ Π, (14)
z j2 ≤ V n2 + w jn + (1 -q jn )Z 2 ∀ j ∈ Γ, n ∈ Φ, ( 15
)
s jn ≥ (V n1 -z j1 ) 2 + ((V n2 + w jn ) -z j2 ) 2 ∀ j ∈ Γ, n ∈ Φ, (16)
z j1 ≤ z k1 ∀ j < k ∈ Γ. (17)
Equation ( 5) represents a budget constraint. As anticipated at the beginning of this section, this constraint is introduced to limit the cost associated with the branches. The bilinear terms q jn s jn will be reformulated at the end of the section to obtain a linear constraint plus auxiliary variables and linear constraints.

Equations ( 6) and ( 7) represent bound inequalities for the branch length.

ECMOR XIV -14 th European Conference on the Mathematics of Oil Recovery Catania, Sicily, Italy, 8-11 September 2014

By means of equations ( 8) each area i is assigned at most to one cluster. Equations ( 9) ensure that area i is assigned to cluster j only if this cluster is active. Equations (10) imply that cluster j is not active if no area is assigned to it.

Equations ( 11) allow cluster j to be connected to well n only if this cluster is active, while equations ( 12) limit the number of branches for each well to B. Additionally, constraint (13) indicates that cluster j can be connected to at most one well, only if the cluster is active.

As mentioned earlier, the oil-bearing areas have to be close enough to the reference point of a cluster to be assigned to it. In such a case, we can expect that a significant amount of the oil contained in these areas is produced through the branch corresponding to the cluster. Equations ( 14) ensure that, if area i is assigned to cluster j, the distance between the cluster reference point and the area cannot be larger than R.

Equations ( 15) represent an angle constraint between the well mainbore and a branch (if cluster j is connected to well n). Note that since all wells are vertical, forcing the height of the well junction point to be greater than the height of the cluster reference point is the same as setting a 90 o upper bound on the branch angle. Equation ( 16) is related to the square distance between cluster reference point j and its junction point on well n. Equations ( 17) establish symmetry breaking constraints, i.e., they impose an order for the lateral coordinates of the clusters1 .

Objective function

The model that is proposed here aims at maximizing the total amount of oil produced. Formally, this quantity can be estimated through the linear objective function: max ∑ i∈Π, j∈Γ,n∈Φ

x i j q jn O i , which accounts for the oil associated with areas belonging to a connected cluster.

Reformulations

We can simplify a number of constraints with the help of auxiliary variables: r jn and u jni replace the bilinear terms q jn s jn ( j ∈ Γ, n ∈ Φ) and x i j q jn , ( j ∈ Γ, n ∈ Φ, i ∈ Π) respectively. The variable r jn is the square of the distance between a connected cluster j and well n (if cluster j is not connected, r jn = 0). The binary variable u jni indicates that area i belongs to cluster j which is connected to well n (if oil area i is not produced, u jni = 0).

The relation between the original and the new variables is given by the following linear constraints:

r jn ≤ s jn + M(1 -q jn ) ∀ j ∈ Γ, n ∈ Φ, (18)
r jn ≥ s jn -M(1 -q jn ) ∀ j ∈ Γ, n ∈ Φ, (19)
r jn ≥ q jn L 2 ∀ j ∈ Γ, n ∈ Φ, (20)
r jn ≤ q jn L 2 ∀ j ∈ Γ, n ∈ Φ, (21)
u jni ≤ q jn ∀n ∈ Φ, j ∈ Γ, i ∈ Π, (22)
u jni ≤ x i j ∀n ∈ Φ, j ∈ Γ, i ∈ Π, (23)
u jni ≥ q jn + x i j -1 ∀n ∈ Φ, j ∈ Γ, i ∈ Π. ( 24
)
where M is defined as the square of the maximum distance between two points in the reservoir.

According to equations ( 18), ( 19), ( 20) and ( 21), if cluster j is connected to well n (i.e. q jn = 1), then the corresponding branch length cannot be larger than L and smaller than L. If cluster j is not connected to well n (i.e. q jn = 0), equations ( 20) and ( 21) imply r jn = 0, and ( 18) and ( 19) are satisfied.

In terms of the new variables, the budget constraint (5) can be rewritten as the linear constraint:

∑ j∈Γ,n∈Φ r jn ≤ H 2 .
In addition, the objective function is now linear:

max ∑ i∈Π, j∈Γ,n∈Φ u jni O i .
Constraints ( 6) and ( 7) are also rewritten as linear constraints, which are constraints ( 20) and ( 21).

As a consequence, the optimization problem formulated in terms of the new variables is a convex MINLP.

It is important to note that in the formulation above we impose no restriction on the crossing or proximity of branches. These constraints are very important in practice but would complicate the problem considerably, and transform it into a non convex MINLP, which is much harder to solve. These constraints could be dealt by means of validation and post-processing stages: identifying violation of crossing constraints and recovering feasibility is easy in two-dimensional case.

Computational results

This section shows preliminary results that illustrate the second phase of our methodology. We investigate the potential of the proposed method using a relatively simple reservoir model based on a twodimensional discretization. The reservoir corresponds to an anticline trap with an aquifer at the bottom to ensure pressure support. The reservoir model has of 830 × 1 × 50 grid blocks, and the dimensions of each grid block are DX = 2.5 m, DY = 10 m, and DZ = 1 m. The reservoir presents two facies with volume fractions 40% and 60%. The mean and variance of the corresponding porosities are 0.4 and 0.005 for the first facies and 0.2 and 0.001 for the second one. A facies realization is first generated using a Fast Fourier Transform Moving Average (FFTMA) algorithm [START_REF] Le Ravalec | The fft moving average (fft-ma) generator: An efficient numerical method for generating and conditioning gaussian simulations[END_REF]) and the Truncated Gaussian method [START_REF] Chiles | Simulation of a nickel deposit: problems encountered and practical solutions[END_REF]; [START_REF] Matheron | Conditional simulation of the geometry of fluviodeltaic reservoirs[END_REF]) with a Gaussian variogram characterized by a horizontal range of 1250 m and a vertical range of 5 m. Thereafter, the two resulting facies are populated with porosity values applying again the FFTMA algorithm. The porosity realizations for each facies are characterized by an exponential variogram. The porosity model is converted into a permeability model using deterministic relationships between the logarithm of permeability K and porosity φ . For the first facies, the relation is log(K) = 5φ + 0.4, and for the second facies, the relation is log(K) = 7φ + 0.3. The resulting porosity and permeability models are shown in Figure 2 and Figure 3.

The fluid flow simulator used is PumaFlow T M (PumaFlow ( 2012)). The initial configuration considered here is not the solution of a previous optimization. It consists of a single vertical producer located at the top of the anticline. Figure 1 represents H u φ S o after a production time-frame of 4 years. The cumulated oil produced is 45.39 thousand cubic meters (285 Mbbl).

In order to limit the number of variables in the model, we selected the high production potential areas, by using an approximation of H u φ S o based on spline interpolation, for areas of 2 × 20, cells and we kept areas with an H u φ S o score above 2.75. These areas are represented in Figure 4.

The mathematical model presented in the previous section is written in AMPL, A Mathematical Programming Language, [START_REF] Fourer | AMPL : a modeling language for mathematical programming[END_REF]), an algebraic modeling language widely used in the MINLP ECMOR XIV -14 th European Conference on the Mathematics of Oil Recovery Catania, Sicily, Italy, 8-11 September 2014 In Figure 5 we show the optimized well trajectory where 5 branches were added. In Figure 6 cubic meters (324 Mbbl). This means an improvement of 13.7% with respect to the initial configuration with one vertical well.

Despite the addition of branches, we can observe on Figure 8, that there is still oil not produced in regions too far to be reached by a branch. Thus, we should revise the model and add the possibility to add vertical production wells when an area is too far from existing wells.

On this very simple case, we analyze the fluid flow simulator outputs corresponding to the end of the production time, but it may be interesting to compute scores at different periods when the oil production starts to decline.

Conclusions

In this paper we have proposed a two-step optimization methodology, to better determine both locations and trajectories of wells (branches) to maximize oil production. Promising results were obtained for a relatively simple example where the number and trajectory of the branches of a production well were optimized.

In future work, we will test the methodology on more realistic cases, the Brugge Benchmark [START_REF] Peters | Results of the brugge benchmark study for flooding optimization and history matching[END_REF]), based on a three-dimensional model with 10 injectors and 20 producers. For this more complex latter problem, heuristics methods, such as Variable Neighborhood Search (VNS, [START_REF] Liberti | A recipe for finding good solutions to minlps[END_REF]) may perform better than Branch-and-Bound in terms of computational speed.

From our first results, the branching problem model should be review to have the possibility to add vertical production wells if remaining oil cannot be reached by a branch associated with an existing well in the neighborhood. This would increase the size of the branching problem, thus we have to better manage the number and relative location of branches for each well. The model could also be improved by accounting for oil along a branch for the computation of produced oil associated with branches.

Another enhancement of the approach presented in this paper could be related to the connectivity measure considered (a distance criterion in this work). Alternative measures could be determined using the permeability field distribution, or a simplified/fast-to-evaluate fluid flow simulator.

Figure 2 Figure 4

 24 Figure 2 Porosity map Figure 3 Permeability map (log)

Figure 5

 5 Figure 5 Optimized well trajectory with high production potential areas.

Figure 6

 6 Figure 6 Optimized well trajectory with initial H u φ S o map.

Figure 7

 7 Figure 7 Oil saturation for initial vertical well configuration after 4 years of production.

Figure 8

 8 Figure 8 Oil saturation for optimized well trajectory with 5 branches after 4 years of production.

Note that these symmetry breaking constraints do not remove any feasible solution (they keep only one configuration when several are equivalent). This limits the feasible space of the solution to be explored by the optimization algorithm considered.

Acknowledgments

We would like to acknowledge David Echeverría Ciaurri for his precious advices. We also thank Eric Flauraud, Guillaume Enchéry, and Steven Thomas for their help in reservoir simulation.

ECMOR XIV -14 th European Conference on the Mathematics of Oil Recovery Catania, Sicily, Italy, 8-11 September 2014