Hervé Serieys
email: serieys@supagro.inra.fr

Hamidreza Nooryazdan

François Kaan

Roberto Bacilieri

Jacques David

Marie-France Ostrowski

Marie-Hélène Muller

André Bervillé
email: berville@supagro.inra.fr

Helianthus annuus natural populations to increase the whole genetic diversity of domesticated sunflower: the concept of neodomestication

Keywords:

To broaden genetic basis of cultivated crop, seventy-seven Helianthus annuus natural ecotypes from North America were grown in Mauguio near Montpellier and crossed with three domesticated sunflower lines. The crossed F 1 progenies were conducted in maternal lineage to maintain 10 plants per lineage in isolation for pollination and the intermixing was continued till G4 generation under low to high selection pressure for some domestication traits. The G0 generation enabled to phenotype thirteen morphologic descriptors. Some traits: phenology, plant height, branching, oil content, seed size and yield were evaluated though the G1, G2, G3 and G4 cycles. The studies of the first generations enabled us to release some main trends which are promising for the future of sunflower breeding. Indeed, we showed that the neo-domesticated materials gather adaptative traits from the natural populations. The first immediate effect of the intraspecific crossing between wild H. annuus and cultivated lines (G0 -> G1) was to increase 1000 seeds weight (1000SW) by 3.6 (9.7 g to 35.2 g), oil content by 20% (26.8 % to 32.3 %), and to reduce apical, intermediate and basal branching, respectively, by 53%, 40% and 13%. In the last G3 and G4 cycles, strong selection pressures were applied on the combined (1000SW*oil content) trait. Oil content increased by 3% in G3 and continued by 11% in G4. In the same time 1000SW, after small decrease, exhibited strong seed size increase (13%) in G4 cycle. The other visible effect was a clear evolution towards domestication traits, i.e., decrease in the different types of branching, reaching 76% non branchy plants in G4. Qualitatively, the different types of branching declined in the same way in G4 cycle: 14% plants were apical branched, 11% intermediate and 17% basal branched.

INTRODUCTION

Wild sunflowers spread naturally in North America where they thrive as populations covering most of the environments from sea to mountain. [START_REF] Heiser | Hybridization between the sunflower species Helianthus annuus and H. petiolaris[END_REF] reported that wild sunflowers were found throughout most of North America on disturbed, mesic, heavy soils that are wet in the spring but dry out by midsummer. The H. annuus natural populations reveal a particular interest for adaptation to biotic and abiotic environments (drought, cold tolerance). Rust, Alternaria and downy mildew resistances are frequently observed in that wild compartment. In screening downy mildew tests performed at INRA we observe that most of the wild accessions (61/72) carried resistance factors to downy mildew (race 710).

Broadly, wild sunflowers did not exhibit meiotic abnormalities in crosses with cultivars like other annual Helianthus species that carry genomic rearrangements in comparison to sunflower. Consequently, one can expect efficient intermingling of the cultivated and wild genomes along recurrent crosses. Moreover, the wild sunflower is strictly self-incompatible, whereas crop lines are self-compatible, which also favours intermingling of genomes.

Sunflower is an oil crop and breeding has dealt mainly with increasing oil yield. Oil yield has been increased drastically in Russia [START_REF] Vrânceanu | Floarea-soarelui hibrida[END_REF] and oil sunflower spread worldwide. In Russia, oil improvement has been accompanied by broomrape resistance since broomrape was prevalent in this country [START_REF] Vrânceanu | Floarea-soarelui hibrida[END_REF]. To become an industrial crop sunflower diversity has been reduced for crop adaptation to hybrid seed production through a CMS [START_REF] Leclercq | La stérilité mâle cytoplasmique du tournesol. I. Premières études sur la restauration de la fertilité[END_REF]) that infers splitting lines in A (male sterile), B (maintainer) and R (restorer) classes, to diseases, mainly downy mildew [START_REF] Leclercq | La stérilité mâle cytoplasmique du tournesol. I. Premières études sur la restauration de la fertilité[END_REF] and after Phomopsis spread (Škorić, 1985). Recently, the introduction of Pervenets mutation in the classic sunflower to obtain seed oil with high oleic acid content has revealed the narrowness of genetic diversity in the crop [START_REF] Tang | Microstellites uncover extraordinary diversity in Native American landraces and wild populations of cultivated sunflower[END_REF].

Until now, the method used to enlarge the genetic diversity for agronomic traits in sunflower crop was to cross sunflower with a wild sunflower carrying the useful trait and to evaluate progenies for the trait under breeding. Similarly, other Helianthus annual species have been used in this way to introgress useful traits in sunflower such as Rf genes: (Jan and Seiler, 2007); Phomopsis [START_REF] Škorić | Sunflower breeding for resistance to Diaporthe/ Phomopsis helianthi[END_REF][START_REF] Griveau | Resistance evaluation of interspecific and cultivated progenies of sunflower infected by Diaporthe helianthi[END_REF][START_REF] Serieys | Etude de la variabilité d'Helianthus annuus sauvages[END_REF], and fertility (Quillet et al., 1992). This method has, as a main advantage, the limitation of the material to be examined in the field, but it also has as a main defect that only one main trait can be considered in each cross. Moreover, introgression lines keep only a few of the foreign genomes, which are eliminated to restore male fertility required in agreement with cultivation. We therefore conceived a method to enlarge genetic diversity of sunflower crop by accumulating, in a series of crosses between wild sunflower and crop lines, neutral diversity from seventy-seven wild sunflower accessions and the main already domesticated traits from three sunflower elite lines. During four successive cycles, we applied selection pressures both to improve the agronomical value of the population and maintain at its highest level the genetic diversity coming from the wild H. annuus compartment. After four generations of mixing in isolated plots, but always keeping separated the female progenies, we evaluated the diversity. We present here the experimental design, the characterization of the main agromorphological traits at each cycle and the analysis of response of important phenotypical and phenological traits in the HAS population as a consequence of the breeding constraints exerted on that genetic material.

MATERIALS AND METHODS

Plant Material

Seventy-seven H. annuus populations (HAS) were chosen among the set of 350 maintained at INRA Breeding station (Mauguio). They were screened for their spreading in the USA by covering most media under different climates and by the diversity of the morphological traits (Table 1). Between 1 to 8 populations were sampled in each of the following 15 US states. 378;388;970;1042;1047;1055;1136 TEXAS 209;509;511;647;648;650;652;734 UTAH 658;774;775;822;826;1149 WYOMING 358;361;955;963;966;974;975

Construction of an intraspecific H. annuus genepool with wide genetic background

The wild and cultivated H. annuus were grown in the field at the INRA breeding station of Montpellier. Sowing dates were on 1 April 1st 1996, in greenhouse with transfer in the field April 15 for G0, April 1 st and 11 April 1997 for G1, April 1998 G2, May 17, 2006 G3, April 17 2007 G4. Five plants per population were grown in our nursery (G0) and each of the plants received separately pollen from the 3 testers, (branching allowed different crosses on the same plant). Testers were oil sunflower lines 89HR2, 90R19 and RT1B11. The two first are top branched (CMS-PET1 restorers) and the third is single headed (CMS-PET1maintainer) We grew 10 G1 plants per progeny (77 accessions x 3 testers x 10 individuals) the next year in a field isolated from gene flow from other sunflowers and we left the plants intercrossing. Each set of 10 plants was harvested separately and the seeds were bulked. At the next generation (G2) we grew up to 52 plants per progeny to enable a smooth elimination of too many branched individuals before they contributed to pollination. At the end of the G2 cycle, an average of 8 plants with reduced or no branching was harvested in each of the cultivated 228 progenies (3 combinations were missing). In order to induce fast increase for the domesticated traits and the agronomical value of the population, the following G3 and G4 generations were obtained after choosing the best individual in each progeny for maximum (OIL * 1000SW) combination value. Like in previous cycles, 5 plants per progeny were grown.

The measured variables

Each G0 plant was phenotyped for 13 morphological and phenological traits reported in Table 2. 1 H1=89HR2, H2=90R19, H3= RT1B11

In the G1, G2, G3 and G4 generations the measured variables mainly concerned flowering period, branching, plant height, seed size, yield and oil content. Branching was defined according to two components (I) types of branching architecture including, basal, medium, apical, and (II) branching intensity scale related to the size of the branches. We used the grid shown in Fig. 1 to note the plants. Analyses were performed according to qualitative and quantitative notation grid presented in Table 3.

G0 G1 G2 G3-BL G4-BL SOWFLO Sow-flowering duration days X X X x X HEIGHT Plant Height cm X X X x X BRAPIC Basal branching 0: absent, 1: present X X X x x BRMED Medial branching 0: absent, 1: present X X X x X BRBAS Apical branching 0: absent, 1: present X X X x x NR Unbranched plants 0: branched, 1: unbranched X X X X X INT_BR Branching Intensity 0 (null) to 4 (high) X X X X X SW Seed weight g X x X x X TSW 1000 seeds weight g x x X x x OIL Oil content % x x x x x

RESULTS

Construction and phenotypic variability in the HAS population

Our principle was to use enough HAS ecotypes to handle a genetic diversity representative of the HAS compartment in North America, to catch from the wild sunflower most of its favourable genes to enhance stress resistance and adaptation in sunflower and to drive the population to decrease the unfavourable traits such as branching, low oil content and seed size. As indicated in Table 4, an important variability is observed for the sow-flowering period among the plant population (more than 2 months), and plant height (from 71 to 308 cm).The branching was also characterised and all the wild Helianthus observed were extremely branched at all levels (apical, intermediate, basal), associated with high branching intensity, Average seed weight (TSW) was estimated at 9.9 g with a range of between 5.7 and 29 g. (close to cultivated types). On the other hand, average oil content of the 77 HAS was low (26.6%) and varied between 20.4 and 32.0%. Seed yield / head in open pollination varied between 0.1 and 9.5 g. The G0 plants were evaluated for 15 traits in the Mauguio environment although they originated from different locations and climates in the USA. A PCA analysis (data not shown) gives the relations between measured phenotypic traits. Plotting of the geographical origin (State) on the graph provides identification of large geographical areas (northern, eastern, south -southwest US), grouping wild H. annuus with similar characteristics. Correlation analysis between main strengths of the climate NR characteristics and morphology of the accessions showed that temperature was highly correlated (r>50%) with height, leaf and petiole lengths. Cycle length and plant height were significantly linked to precipitation (r=0.4) and temperature levels (r= 0.28 and 0.57 respectively). Additionally, plant height is connected to important rainfalls in short periods (coastal Texas climate). The apical and intermediate branching also seems to be affected by rainfall and to a lower number of rain days. Increased temperature range (maxi-mini) enhanced architectures with intermediate branching. The native accessions from geographical areas with a low rainfall (New Mexico, Arizona) were more branched. On the contrary, the different types of branching did not clearly respond to temperature.

Changes observed from G0 to G4 for main traits (Fig 2):

Changes in sowing-flowering period and plant height, plant branching, in the average oil content for seed size, seed yield were plotted. Although the experiments were not done the same year, we observed some general trends. DISCUSSION Systematic series of crosses between HAS and elite lines of domesticated sunflower have never been performed. We did not have any problems in the intercrosses, except some male sterility that segregated in some HAS. Here we reported several main features that tended to show that the material could be used for breeding sunflower and we tried to quantify our conclusions as much as possible. Some traits not retained in our notations were strongly selected such as seed dormancy, which was rapidly eliminated. Although HAS all displayed strong seed dormancy, we did not observe any problems in our progenies. This means that the trait is probably monofactorial and the locus efficiently eliminated. Branching behaved as a complex trait: we distinguished branching at 3 levels: basal, medium and top (Fig. 2). Basal and medium levels are unfavourable in the crop whereas top branching is used in most restorer lines. We surmised that the length of the branches was not inherited but due to the environment; the later the sowing, the shorter the branches for some lines and ecotypes. Non branching plants appeared in the G1 and due to selection pressure for the traits, it increased continuously. We used the grid (Fig. 1) to note the plants. We observed all the combinations of the elementary traits, which suggested that they are controlled by different genetic factors mostly independent. Seed weight was much lower in HAS (4-15g for 100 seeds) than in the crop (50-70g). However, when we chose the HAS we retained some with higher seed weight up to 30g. Two explanations are possible: 1) this is the range of variation for seed weight in HAS or, alternatively, 2) they have already been introgressed by domesticated sunflower. Due to the origin of the HAS with a large seed size in areas where the sunflower is cropped, we favoured the second hypothesis. However, in our experiment we used a combination of seed weight and oil content to avoid screening progenies with large seeds and poor oil content. Seed oil content did not increase rapidly, but we did not exert any strong selection pressure on the trait. Plant height decreased rapidly without any selection pressure. The variability is wide and height reduction will probably be rapidly feasible.

Fig. 1 .

 1 Fig. 1. Notation grid mixing 'Type' and 'Intensity' for branching of HAS and sunflower.

Fig. 2 .Fig

 2 Fig. 2.Changes in phonological and morphological parameters along five generations in the neodomesticated progenies. SOBLO, sowing -blossoming stages in days. NR: non branching. Two G3 progenies were obtained in two different fields: One with ten plants per family with all 228 families (G3BL+), which has also led to G4BL+; and one with only 100 families (G3-Syspro)

Table 1 .

 1 Distribution, by states, of the 77 collected accessions in North America

	US State	H. annuus accessions with INRA CODE
	ARIZONA	383; 386; 519; 654
	CALIFORNIA	410; 421; 435; 437; 446; 458; 468; 833
	COLORADO	363; 660; 980; 989; 996; 1147; 1148
	ILLINOIS	211
	IOWA	665; 829
	KANSAS	733; 997; 998; 999; 1000
	MISSOURI	351
	MONTANA	943; 945; 948; 954; 1150
	NORTH DAKOTA	928; 929; 931; 933; 939
	NEW MEXICO	461; 463; 649; 661
	OKLAHOMA	646; 651; 662
	SOUTH DAKOTA

Table 3 .

 3 Quantitative and qualitative characters measured on the HAS population through G0 to G4 generations

	HAS

Table 2 .

 2 History of selection pressures applied on the wild H. annuus Pool (HAS), built on 1996

	Cycle /year	Operation at Melgueil	Size of populations	Traits examined	Selection pressure
	Cycle 0	F1 construction in field	77 HAS	Description origin,	None
	(G0)	HAS= female parent 5	pollinated by 3	agromorphology +	Harvest; only F1 hybrid kernels,
	1996	plants from each	sunflower lines	phenology + kernel size
		accession
	Cycle 1	Intermixing in isolation	228 F1 hybrid	Branching type intensity,	Harvest: Whole plants 1
	(G1)	(=77 ecotypes x 3	combinations	leaf shape, pollen fertility,	* H1 (=10 plants separately)
	1997	testers x 10 plants)	with 10 plants	head diameter, earliness;	* H2, H3 (= 10 plants in bulk)
			per progeny	kernel yield, oil content,	G1 to G2: weak selection (24 %
			2280 F1 plants	TSW, seed yield.	plant strongly branched were
					discarded for H1). None for H2
					and H3
	Cycle 2	Intermixing in isolation	228 maternal	Height + blooming date	Before blossom Elimination of
	(G2)	228 maternal progenies	progenies	(start, medium, end)	most branched plants (bush form).
	1998	(776 X 3 testers	11800 plants	On harvested plants:	At harvest: choice of 8 plants less
				Branching, high, TSW, oil	branched. On the average, 84 %
				content, Total W	plants were eliminated in G2
	Cycle 3	Intermixing in isolation	228 maternal	Notations: blooming date,	For each G2 progeny: the best
	(G3-	(77 x 3 testers)	progenies with	height, branching, TSW,	OIL content * TSW plants were
	BL)		5 plants	oil content, seed yield	selected for G3. 80 % plants
	2006.		1140 plants		eliminated
	Cycle 3	Intermixing in isolation	228 maternal	Notations: blossoming	For each G3 progeny: the best
	(G4-	(77 x 3 testers)	progenies with	date, height, branching,	OIL content * TSW plants were
	BL)		5 plants	TSW, oil content, seed	selected for G4. 80 % eliminated
	2007		1140 plants	yield	plants.

Table 4 .

 4 Agro-morphological characteristics of the selected 77 accessions, in the G0 HAS population

	Trait

Proc. 17 th International Sunflower Conference, Córdoba,Spain (2008)


Magdalena Kristiawan

Tomáš Jirout

Alioune El

Vaclav Sobolik
email: vsobolik@univ-lr.fr.

Tomas Jirout

Vaclav Sobolík

Wall shear rate components for wavy Taylor-Couette flow

Introduction

Several techniques can be used for the measurements of shear rates at the solid-liquid interface. The measurement of the current of an electrochemical reaction, i.e. the electrodiffusion (polarographic) method [START_REF] Cognet | Utilisation de la polarographie pour l'étude de l'écoulement de Couette[END_REF][START_REF] Lebouché | Relation entre les fluctuations pariétales du transfert massique et du gradient de vitesse dans le cas d'un nombre de Schmidt grand[END_REF], is the only method able to follow the space distribution and time evolution of wall shear rates. This method is based on the measurement of the electric current limited by the transfert of actives species to the working electrode. Using a redox reaction, repetitive measurements can be carried out with a system which does not require addition of active compounds as in the case of the above mentioned method. The measurements are confined to the diffusion boundary layer which is very thin (several micrometers) in the electrolyte solutions (high Schmidt number). This fact distinguishes the electrodiffusion method from the other methods like laser Doppler anemometry, hot film anemometry and particle image velocimetry. The three-segment probes (see Fig. 1) are composed of three insulated circular sectors that act as working electrodes [START_REF] Sobolik | Method of manufacture of segmented probe with circular cross-section[END_REF]. For the evaluation of flow direction, we make use of the non-linear dependence of the limiting current density on the distance from the forward electrode edge, i(x) ~ x -1/3 [START_REF] Lévêque | Les lois de la transmission de la chaleur par convection[END_REF]). The forward sector has higher current than the sector that lies in his shade. Hence, shear rate and its components can be calculated from the limiting diffusion currents passing through the three segments [START_REF] Wein | Theory of direction sensitive probes for electrodiffusion measurement of wall velocity gradients[END_REF]. In this work, the electrodiffusion method with a three-segment and simple microelectrodes was used for the study of the axial ( z ) and azimuthal ( ) components of shear rate at the outer wall of Taylor-Couette flow.

Experimental

Experiments were carried out in the flow between concentric cylinders with the inner cylinder rotating and outer cylinder fixed. The inner cylinder was driven by a stepper motor with electronically controlled rotation rate. The diameters of the outer and inner cylinders were equal to 62 and 49.6 mm, respectively. This corresponds to the radius ratio, R1/R2 = 0.8. The aspect ratio was about 44. For simultaneous observation of flow pattern and measurement of shear rate and vortex characteristics a Kalliroscope solution was added to the electrolyte. The flow regime with wavy Taylor vortices, i.e. vortices with steady azimuthal waves was studied. The vortices were swept past the electrodes by a slow axial flow. Using three-segment probe, see Fig. 1, complete maps of axial and azimuthal components of shear rate and their fluctuations at the wall of the outer cylinder were obtained as a function of axial co-ordinate and time. The vortex size (height) and the number and celerity of azimuthal waves were measured by an array of simple electrodes mounted in the wall of the outer cylinder. The arrangement of the electrodiffusion probes along the perimeter of the outer cylinder is shown in Fig. 2. The azimuthal ( ) and axial ( z ) components of the shear rate ( ) were measured by the three-segment probe 3s with a diameter of 0.5 mm. Simple probes e4, e5 and e6 were formed by circular electrodes with the same diameter (0.5 mm). The celerity of azimuthal waves was calculated from the time lag (correlation) of currents delivered by e4 and e5. The wave period was obtained from autocorrelation of the current through the electrode e4. The wavelength was calculated from wave period and celerity. The vortices were driven by a slow axial flow (z-direction) along the probes. The drifting vortex velocity was calculated from the time necessary for a vortex to cover the distance between e5 and e6. The height of vortices was obtained from their period and drifting velocity. The correlations and fast Fourier transformation of measured currents were calculated using Scilab software

L 3 L 1 L 2 R 2 3s e4

Results and discussion

The limiting diffusion current histories measured by the three-segment probe (I1, I2, I3) and three simple probes (I4, I5, I6) at Ta=323 are shown in Fig. 3. The measurements were taken over about one and half vortex pairs with a sampling frequency of 40 1/s. The vortex drifting velocity was 0.14 mm/s. Knowing the drifting velocity, the time axis can be transformed to the axial coordinate z. The time td55 corresponds to the passage of a vortex pair along the probe e5 in axial direction, whereas the time td56 indicates the passage of vortex crest between the probes e5 and e6. The time td56 was used for calculation of vortex drifting velocity and td55 for calculation of vortex height. The azimuthal waves manifest themselves in the signal oscillations. These waves were studied by means of the current correlations R4i, calculated with respect to I4. The time lag tw45 between the autocorrelation of R44 and correlation R45 (see Fig. 4) corresponds to the passage of a selected wave point between e4 and e5, and tw44 to the passage of one wave past e4. The time tw45 was used for the calculation of azimuthal wave celerity and the wave period tw44 for calculation of wavelength. The number of azimuthal waves was obtained by dividing the perimeter of the outer cylinder by the wavelength. The shear rate components calculated from currents I1, I2, I3 are shown in Fig. 5a. The z-coordinate concerns the height of vortices but it is irrelevant for the wavelength of azimuthal waves. The period of azimuthal waves shown in Fig. 5a is a function of number of waves, angular velocity of waves and drifting velocity of vortices. The amplitude of the shear rate components at different position z is given by the envelope of the oscillations which is not shown. The meaning of different symbols is also shown in this figure. The notion of filtered shear rate is introduced. This means the shear rate or its components from which the azimuthal waves were filtred out by a Butterworth filter, see f and zf in Fig. 5a. The course of filtered wall shear rates and their components is similar to the wall shear rates in Taylor-Couette non wavy flow (Sobolik et al.). From the point of view of the axial velocity, the outflow creates a forward critical point on the outer cylinder, see Fig. 5a. The azimuthal component of shear rate exhibits a maximum at this point and a minimum at the rear critical point (inflow). The axial component of shear rate is zero at these points. In fact, there is a small finite value of axial shear rate due to the drifting velocity of vortices. The amplitude of azimuthal waves, given by the envelope of oscillations, has azimuthal a , and axial component a z . The azimuthal component is zero in forward and rear critical points. The axial component has a maximum in forward critical point and it is very small at the extremes of zf . Like in the classical forward critical point, the amplitude of axial oscillations, a z , exhibits maxima and the variation of zf is steep and almost linear. It is much steeper than in the inflow. The outflow has the form of a narrow jet with a high radial and azimuthal velocity. Even if the amplitude of the axial component of shear rate has its maximum at the forward critical point, the oscillations of absolute value shear rate are very small, see Fig. 5b. It is due to the high value of azimuthal shear rate which is steady at this point. We observed the transition from Taylor vortex flow into wavy vortex flow at Ta= 80.5. This value is much higher than the values calculated by Jones 1985. From his Fig. 3, we deducted values of 55.0, 57.6 and 64.5 for set-up of 1, 2 and 3 waves, respectively. We found always three waves at the transition which is in agreement with the simulations carried out by Jones. The components of wall shear rate were measured at different rotational rate (Ta=97,202 and 483). At the lowest rotational rate the regime of three azimuthal waves was stable. The characteristic values of the vortices are shown in Table 1. The axial component of the wall shear rate can be compared to the wall shear rate in classical impinging jet even if the azimuthal component in the forward stagnation point exhibits a maximum. According to Prandtl's boundary theory [START_REF] Schlichting | Boundary Layer Theory[END_REF]), it holds for the stagnation region of an impinging jet that the wall shear rate is proportional to the distance from the stagnation point. The maximum of the wall shear rate is at the transition from the stagnation zone to the wall jet region. The spatial distribution of the wall shear rate fluctuations depends on the jet Reynolds number and the distance of the orifice from the plate (Alekseenko and Markowich 1994). It can have a maximum or minimum in the stagnation point or a constant value in the jet region. At Ta=97, the amplitude of z had maximum at the forward stagnation point and with increasing Ta the amplitude became almost constant in the outflow zone. [START_REF] Jones | The transition to wavy Taylor vortices[END_REF] concluded that the strong azimuthal jets at the outflow region destabilize the flow making the vortices wavy. [START_REF] Coughlin | Modulated waves in Taylor-Couette flow Part 2. Numerical simulation[END_REF] supposed that the both radial and azimuthal jets and axial gradient of the azimuthal velocity are responsible for the waviness. The strong azimuthal jet manifests itself by the maxima of f (see Fig 5a). It increases with the Taylor number. The strength of radial jet can be deduced from the axial gradient of zf . Its maximum value is also in the forward critical point.

0 1 2 -1 0 1 R 4i R 4 R 42 R 41 R 45 R 46 R 44 t w45 t w44 t, s
The shear rates and their components were calculated from the measured limiting diffusion currents using the quasi-steady approach, i.e. the Lévêque solution. Taking into account the inertia of the diffusion boundary layer, the amplitudes of oscillations will be even higher.

Conclusions

The three-segment ellectrodiffusion probes are convenient for the mapping of the shear rate components on the wall of the wavy Taylor-Couette flow.

The form of the axial distribution of filtered wall shear rate components is similar to the distribution in steady Taylor vortices. However, the maximal values are more important. Jet like outflow characterized by forward critical point was identified. The oscillations of axial component of shear rate exhibit maximum oscillations whereas the oscillations of the azimuthal component are negligible in the forward critical point.

The effect of the inertia of the boundary layer on the oscillation amplitudes should be studied in the future

  Fig. 1. Three-segment probe Fig. 2. Arrangement of the three-segment probe 3s and with a diameter of 0.5 mm. simple electrodes e4, e5 and e6 along the perimeter.

Fig. 3 .

 3 Fig. 3. Current histories: I1, I2, I3 -three-segment probe; I4, I5, I6 -simple electrodes; Ta =323

Fig. 4 .

 4 Fig. 4. Correlations with respect to the current I4

Acknowledgements

This work was supported by the "Agence National de la Recherche", France, by the grant number ANR-08-BLAN-0184-01.


Rodolphe Turpault
email: rodolphe.turpault@univ-nantes.fr

Mathematical Modelling and Numerical Analysis Will be set by the publisher Modélisation Mathématique et Analyse Numérique AN HYBRID MODEL FOR RADIATIVE TRANSFER

Keywords: 1991 Mathematics Subject Classification. 35Q99, 80A99, 85A25. April 24, 2008 Radiative transfer, kinetic models, moments models, entropic closure, asymptotic preservation, mathematical modelling

come

An hybrid model for radiative transfer

Introduction

The radiative transfer equation (RTE) describes the evolution of the radiative intensity I ν (Ω) = I(t, x, Ω, ν) where Ω is the photons' direction of propagation and ν is the frequency. This radiative intensity I ν (Ω) is linked to the photons' distribution function. Assuming local thermal equilibrium and neglecting scattering, the RTE writes:

1 c ∂ t I ν (Ω) + Ω.∇I ν (Ω) = σ ν B ν (T ) -I ν (Ω) , (0.1)
where c is the speed of the light, σ ν is the emission opacity, T is the material temperature and B ν (T ) is Planck's blackbody function given by:

B ν (Ω) = 2hν 3 c 2 exp hν kT -1 -1 . (0.2)
The constants h and k are respectively Planck and Boltzmann constants. It is to note that the energy of B (ie its integral over all directions and frequencies) is proportional to T 4 . The constant of proportionality is denoted by a( 7.56.10 -16 SI). Due to emission and absorption process, the radiative energy is not conserved. To preserve the conservation of the total (radiative+material) energy, the following simplified temperature evolution equation will be used in this paper:

ρC v ∂ t T = ∞ 0 S 2 σ ν I ν (Ω) -B ν (T ) dΩdν. (0.
3)

The direction of propagation Ω and the frequency ν play a very different role in the RTE. Indeed, the opacity σ ν is a critical parameter. When σ ν = 0 the RTE is nothing but a convection equation and the associated regime is called the transport (or free-streaming) regime. On the other hand, when σ ν → ∞, the RTE degenerates into the diffusion equation (see [START_REF] Larsen | Asymptotic analysis of radiative transfer problems[END_REF], [START_REF] Pomraning | The Equations of Radiation Hydrodynamics[END_REF] and [START_REF] Mihalas | Foundations of Radiation Hydrodynamics[END_REF]):

∂ t (ρC v T + 4aT 4 ) -∇ x 4acT 3 3σ ∇ x T = 0. (0.4)
A lot of simulations involve free-streaming, diffusion as well as in-between regimes altogether. It is therefore crucial to insure that the models we use correctly handle all regimes. Moreover, σ ν has huge frequencial variations due to phenomena such as bound-bound transitions which produce emission lines that are very narrow. Thus grey models (ie integrated over all frequencies) are often too coarse and introduce significant errors due to frequency averaging. The models we use hence have several frequency degrees of freedom. The frequencies are split into groups and bands respectively indexed by q and k:

[0, +∞[ ν = Q q=1 Q q K k=1 Q k .
As a notation and to avoid confusion, every frequency group in this paper will be indexed by q and every frequency band will be indexed by k. Inside the bands Q k , a quadrature formula is used to approximate integrals over directions and frequencies. The following notations will be used throughout this paper:

< • > q = 1 c Qq S 2 • dΩdν, < • > = K k=1 L l=1 • ω l ξ k 1 c ∪Q S 2 • dΩdν, < • > = Q q=1 < • > q +< • >.
Throughout this paper, frequency and direction indexes will be omitted inside integrals and quadrature formulas except when necessary.

• Bands refer to frequency intervals that are small enough to make several assumptions. For instance, Planck's function may be supposed to be constant inside a band. Among others, these assumptions allow to use several useful techniques to compute mean values of σ ν (eg CK methods). A typical simulation may take into account thousands of frequency bands. • Groups refer to large frequency intervals where the previous assumptions cannot be made. A particular attention must then be paid to the way to compute mean values of σ ν . A typical multigroup simulation may consider a few frequency groups (up to several dozens).

The computation of the opacities' mean values inside groups may be problematic for modelling issues as we will see later. Among the classical choices to do so are Planck and Rosseland mean values, respectively denoted by σ P and σ R and given by:

σ P = σ P (T ) = < σB(T ) > < B(T ) > , (0.5) σ R = σ R (T ) = < ∂ T B(T ) > < 1 σ ∂ T B(T ) > . (0.6)
A lot of physical applications involve radiation in a wide range of regimes. It is possible to distinguish between two categories considering the influence of radiation on the flow:

• The radiation impact on the flow is neglectable. This is the case for instance of climatology applications.

• The radiation effects have an impact on the flow. This happens in superorbital atmospheric reentry, astrophysics (radiative shocks,. . . ) and ICF (fusion obtained by impacting laser beams on a deuterium/tritium target). These two categories imply very different levels of modelling. For the latter, a full coupling between hydrodynamics and radiation may be mandatory and requires the use of a robust and cheap radiative model. The objective of this paper is to introduce an hybrid model which mixes a kinetic and a moments model. This hybrid model shall be supple enough to naturally adapt to any simulation, even though the main target applications requires a full coupling between radiation and material and a specific geometry (eg convergence of lasers for ICF). This article will be organized as follows. In the next section, we will briefly recall a kinetic model suitable for a narrow-band level of frequencial resolution. We will mainly focus on its construction and its main properties. Then a section will be dedicated to the M 1 model. The multigroup version will be recalled, we will also see that it has several important properties that make it robust and a very interesting candidate for multigroup simulations. We will focus on one of its most essential characteristic which is the possibility to consistantly handle the opacities' mean values. Finally, we will introduce the way of coupling these two models inside an hybrid model. We will then show that this model still preserves two fundamental properties: the asymptotic behaviour in the diffusion regime and the local decreasing of the total entropy. A simple numerical example will illustrate its interest.

A kinetic model (Q = 0)

This model was first introduced in [START_REF] Charrier | Discrete-velocity models for numerical simulations in transitional regime for rarefied flows and radiative transfer[END_REF]. It is designed for using frequency bands. One of the consequences of this choice is that the mean values of the opacity σ k do not depend on variations of B or I. Therefore, σ k has the same value that Planck (0.5) or Rosseland (0.6) mean values. Aside from these frequency bands, we also consider a finite number of directions (Ω l ) l=1...L and look for an approximation of I ν k (Ω l ) = I(t, x, Ω l , ν k ) denoted J k,l which is solution of the following discrete system:

1 c ∂ t J k,l + Ω l ∇J k,l = σ k B k (T ) -J k,l ], (k, l) ∈ [1 . . . K] × [1 . . . L] (1.1) ρC v ∂ t T = < cσ(J -B(T )) >. (1.2)
The main difficulty is to carefully define the discrete equilibrium function B k (T ). The "naive" choice B k (T ) = B ν k (T ) is not sufficient enough to insure the conservation of several important physical properties. For example, it is easy to see that in general 4 . Therefore, some energy is artificially created (or lost) only due to this form of the discrete equilibrium function. A solution is to use Levermore's technique of minimum entropy principle [START_REF] Levermore | Moment Closure Hierarchies for Kinetic Theories[END_REF] in a fashion related to Mieussens' choice for Boltzmann-BGK [START_REF] Mieussens | Discrete velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries[END_REF]: the discrete equilibrium function is chosen to be the minimum of the discrete radiative entropy:

K k=1 L l=1 B ν k (T ) ω l ξ k = aT
B = (B k,l (T )) k,l = argmin{< h d (J) >, < J > = aT 4 }, (1.3)
h d (J) = 2kν 2 k c 3 n k,l ln(n k,l ) -(n k,l + 1) ln(n k,l + 1) , n k,l = n(J) k,l = c 2 2hν 3 k J k,l .
This choice is consistent in the sense that Planck's function is the minimum of the radiative entropy among all the functions whose energy is aT 4 . Therefore (B k,l ) k,l can be considered as an approximation of B ν .

Theorem 1.1. If T > 0 then there exists a unique B such that:

B k,l = B k = 2hν 3 k c 2 exp hν k α k -1 -1 , (1.4)
where α is the unique solution of < B > = aT 4 .

Proof. This results comes from a straightforward computation, which is briefly (and formally) recalled here.

Let us introduce the Lagrangian of the problem (1.3):

L(J, λ) = < h d (J) > -λ(aT 4 -< J >). (1.5)
According to Lagrange's theorem, if B is a solution of (1.3) then there exists α ∈ IR such that (B, α) is a saddlepoint of L. Moreover, thanks to the convexity of the discrete entropy, we can say that if (B, α) is a saddle-point of L, then it is a solution of (1.3). At the saddle-points of L, if L is derivable, we have ∂ J L(B, α) = 0 hence:

∀ϕ ∈ D(S 2 × IR + ), < < h d (J) >, ϕ > D ,D = -< < α >, ϕ > D ,D .
In particular, the last equality stands for all ϕ ∈ D(S 2 l × Q k ), where S 2 l is the part of the unit sphere around the direction Ω l . Thus we can say that:

∀ k, l, h d (B k,l ) = -α, k hν k ln(1 + 1 n(B) k,l ) = -α, hence B k,l = 2hν 3 k c 2 exp hν k α k -1 -1
.

There are two points that have to be emphasized here. First, the discrete equilibrium function -as Planck's function-is isotropic: B k,l = B k . Among others, a consequence of this is that only the sum over k is required to compute α. Lagrange's multiplier α = α(T ) plays an important role. It is a scalar regardless of the number of frequency bands considered and its determination requires to solve the nonlinear scalar equation < B(α) > = aT 4 . If the quadrature formula was exact, then we would have α = 1/T . Using 1/T as a starting value for any suitable efficient numerical method allows to converge quickly, hence the determination of α does not involve expensive computations.

Proposition 1.2. The system (1.1)-(1.2) is hyperbolic symmetrizable. Moreover, if the quadrature formula < • > is wisely chosen then it has the following properties:

• the total energy is conserved:

∂ t (< J > + ρC v T ) + c ∇ x < Ω l J > = 0, • if J(t = 0) ∈ L ∞ then ∀t > 0, J(t) ∈ L ∞ ,
• the total entropy is locally decreasing:

∂ t < h(J) > c + ρC v U + ∇ x < Ω l h d (J) > ≤ 0.
where U = -T 0 α(T )dT is a convex function of T . • it is asymptotic preserving in the diffusive limit.

The first two results above were proved in previous papers (see [START_REF] Charrier | Discrete-velocity models for numerical simulations in transitional regime for rarefied flows and radiative transfer[END_REF] and [START_REF] Turpault | Modélisation, approximation numérique et applications du transfert radiatif en déséquilibre spectral couplé avec l'hydrodynamique[END_REF]). The last two will be proved for the more general hybrid model -which includes this kinetic model-later in this article.

Unfortunately, the number of bands required for any kinetic model implies that these models are way to expensive in a fully coupled context involving radiation, hydrodynamics and other processes such as chemical reactions or turbulence.

Multigroup-M 1 model (K = 0)

Since our main objective is to perform fully coupled simulations, we have to consider cheaper models ie models that are somehow integrated over directions and frequencies. There are two categories of such models: flux-limited diffusion and moments models (see [START_REF] Su | An analytical benchmark for non-equilibrium radiative transfer in an isotropically scattering medium[END_REF] and references therein for a list of the classic choices). Our choice is to select the M 1 model, which belongs to the category of moments models. It was introduced in [START_REF] Dubroca | Entropic Moment Closure Hierarchy for the Radiative Transfer Equation[END_REF] and has several variations among which [START_REF] Turpault | A consistant multigroup model for radiative transfer and its underlying mean opacities[END_REF], [START_REF] Dubroca | Half Moment closure for radiative transfer equations[END_REF], [START_REF] Turpault | Multigroup Half Space Moment Approximations to the Radiative Heat Transfer Equations[END_REF] and [START_REF] Frank | Partial moment entropy approximation to radiative heat transfer[END_REF]. To built it, the first step consists in obtaining the moments equations from (0.1). To do so, an integration over directions and inside each frequency group Q q is performed for (0.1) and Ω × (0.1). It leads to the following system: q = 1 . . . Q, ∂ t E q + ∇F q = σ e q aθ 4 q -σ a q E q , (2.1)

∂ t F q + c 2 ∇P q = -σ f q F q . (2.
2)

The radiative energy E q , the radiative flux vector F q and the radiative pressure tensor P q are the first three moments of the radiative intensity inside the q th frequency group. Moreover, σ e , σ a and σ f are mean values of the opacity and their choice will be discussed later. These variables are defined as:

E q =< I > q , F q =< cΩI > q , P q =< Ω ⊗ ΩI > q , aθ 4 q =< B(T ) > q , σ e q = < σB(T ) > q aθ 4 q , (2.3)
σ a q = < σI > q E q , (2.4)
σ f q =     < σΩ x I > q F x q 0 0 < σΩ y I > q F y q     .
(2.5)

The system (2.1)-(2.2) is the generic multigroup moments system. It is not closed. To close it, one has to make an hypothesis on the radiative pressure P q so that it can be expressed as a function of E q and F q . The choice of the closure determines the model. There are several ways to choose the closure, including the classic P 1 model where P q := I d E q /3. This shape is based on the fact that at the radiative equilibrium (ie whenever I ν (Ω) = B ν (T )), the radiative pressure is indeed one third of the radiative energy times the identity tensor. The main drawback of the P 1-model is that it does not predict physical solutions far from the radiative equilibrium.

Definition 2.1. A solution (E q , F q ) will be said to be physically relevant if E q > 0 and F q ≤ cE q .

In fact, this definition is simply the consequence that E q and F q are the first two moments of some (positive) radiative intensity. The second property is called the (radiative) flux limitation. It is very important yet hard to preserve in models. For instance, the P 1 model violates this property and thus sometimes predicts unphysical solutions. There are several other options to close the system (2.1)-(2.2) (see [START_REF] Su | An analytical benchmark for non-equilibrium radiative transfer in an isotropically scattering medium[END_REF] and references in [START_REF] Turpault | A consistant multigroup model for radiative transfer and its underlying mean opacities[END_REF]). One of which is to use a minimum entropy principle in a fashion similar to what Levermore did for fluid mechanics [START_REF] Levermore | Relating Eddington Factors to Flux Limiters[END_REF]. It is to note that other models based on a minimum entropy principle exist to close this system (see for example [START_REF] Minerbo | Maximum entropy Eddington Factors[END_REF]), however the M 1 model is based on the (physical) radiative entropy:

h(I) = 2kν 2 c 3 n I ln(n I ) -(n I + 1) ln(n I + 1) , (2.6
)

n I = c 2 2hν 3 I ν (Ω).
(2.7)

The first step is to build an underlying radiative intensity I thanks to the minimum entropy principle:

I = argmin{< h(I) > / ∀q = 1 . . . Q < I > q = E q , < cΩI > q = F q }. (2.8)
The solution of this minimization problem is:

I = q 1 q 2hν 3 c 2 exp( hν kT α q .(1, Ω) ) -1 -1 , (2.9)
where 1 q is the characteristic function of the q th group and α q is the Lagrange multiplier of the minimization problem (fully determined by the constraints).

The M 1 form of the radiative pressure is simply the corresponding moment of the underlying radiative intensity:

P q :=< Ω ⊗ ΩI > q .
P q may be expressed in Eddington form as follows:

P q = D q E q , (2.10)
D q = 1 -χ q 2 I d + 3χ q -1 2
F q ⊗ F q F q , (2.11)
where the Eddington factor χ q is the eigenvalue of D q associated with F q .

Proposition 2.2. The system (2.1)-(2.2) with the M 1 choice of the closure is hyperbolic symmetrizable. Moreover, it has the following properties:

• the total energy is conserved,

• the total entropy is locally decreasing,

• it is asymptotic preserving in the diffusive limit,

• it has a natural flux limitation inside each group (ie

F q ≤ cE q ).
Once again, these properties were proved in previous articles ( [START_REF] Turpault | A consistant multigroup model for radiative transfer and its underlying mean opacities[END_REF], [START_REF] Turpault | Modélisation, approximation numérique et applications du transfert radiatif en déséquilibre spectral couplé avec l'hydrodynamique[END_REF]) or will be proved for the more general hybrid model in the next section of this article. One of the edges the M 1 model has compared to other related models is the possibility to chose consistent mean opacities. As a matter of fact, the opacities' mean values required for moments or flux-limited diffusion models are problematic since their expression should be given by (2.3), (2.4) and (2.5). On the one hand, (2.3) is nothing but Planck's mean value of the opacity which is relatively easy to compute. But on the other hand, (2.4) and (2.5) are mean values relatively to functions of I, which is not reachable (the unknowns of our model are its moments). As a consequence, one has to make an additional assumption to compute these mean opacities. Classic choices include taking Planck, Rosseland or Chandrasekhar mean values but this additional assumption has a cost in terms of precision of the model. Hopefully for the M 1 model, even if the real radiative intensity cannot be used, we have built an underlying one I which may be considered to approximate (2.4) and (2.5). For instance:

σ a q < σI > q E q .
(2.12)

Since I is fully determined by the knowledge of E q and F q , it is possible to make all the computations. This choice is consistent in the sense that no additional assumption is made. It allows to have a good precision even with a few number of groups: see the test-cases in [START_REF] Turpault | A consistant multigroup model for radiative transfer and its underlying mean opacities[END_REF].

Using the consistent mean values of the opacity, the multigroup-M 1 model is equivalent to the following formulation:

< 1 c ∂ t I + Ω∇I > q =< σ B(T ) -I > q , q = 1 . . . Q (2.13) < 1 c ∂ t (ΩI) + ∇(Ω ⊗ ΩI) > q = -< σΩI > q ,
(2.14)

I q = 1 q I q (E q , F q ), (2.15)
where I q is the underlying radiative intensity used as a closure function for the M 1 model (2.9).

The multigroup-M 1 is usually a good choice for full coupling purposes. It is relatively cheap, very robust, conserves the fundamental physical properties and even has a natural way to approximate mean opacities. However, it has a directional drawback since it has to cope with only one direction of anisotropy per group (colinear to F q ): it cannot handle the case when at a given point two or more comparable energies converge from opposite directions (see [START_REF] Dubroca | Half Moment closure for radiative transfer equations[END_REF] for an example). Depending on the geometry of the problem, this drawback may be problematic or not. For instance, for computing the fluxes on the front shield of a probe or to simulate an explosion, this directional effect has nearly no impact. The M 1 model has been declined into half and partial space M 1 models [START_REF] Dubroca | Half Moment closure for radiative transfer equations[END_REF], [START_REF] Turpault | Multigroup Half Space Moment Approximations to the Radiative Heat Transfer Equations[END_REF], [START_REF] Frank | Partial moment entropy approximation to radiative heat transfer[END_REF]. For 1D simulations, the half-M 1 model (ie integrated over half of the unit sphere for directions) nearly annihilate the directional effect, but for 2D and 3D simulations, even partial M 1 might not be sufficient. An example of tricky simulation would be a ICF-type problem where dozens of lasers converge on 1 mm-radius zone.

Hybrid Model

As we have seen in previous sections, both the kinetic and the multigroup M 1 models may have some drawbacks. The former is often too expensive and the latter sometimes has directional problems. Starting from the observation that these two models are complementary, we derive an hybrid model which mixes the multigroup-M 1 model inside frequency groups and the kinetic model inside frequency bands:

∀q = 1 . . . Q, ∂ t E q + ∇.F q = c σ e q aθ 4 q -σ a q E q , (3.1)
∂ t F q + c 2 ∇P q = -cσ f q F q , (3.2)
∀k = 1 . . . K, ∀l = 1 . . . L, 1 c ∂ t J k,l + Ω l ∇J k,l = σ k B k (T ) -J k,l ], (3.3)
ρC v ∂ t (T -T ) = Q q=1
(σ a q E q -σ e q aθ 4 q ) (3.4)

ρC v ∂ t T = < cσ(J -B(T )) > (3.5)
I q = argmin < h(I) > q /∀q = 1 . . . Q, < I > q = E q and < cΩI > q = F q (3.6)

B = argmin < h d (I) >/< I > = aT 4 - Q q=1 aθ 4 q }. (3.7)
Where T = T (T ) is such that:

T = < B(T ) > a 1/4 . (3.8)
This model may be seen as an improvement of either the M 1 or the kinetic model. Typically, groups may be chosen in regions were either very little energy exists or there is no directional problem whereas bands are chosen where a large amount of energy lies inside a few narrow bands (eg lasers) or may converge at a given point.

Of course, it is crucial to insure that this hybrid model behaves well. The conservation of the energy is straightforward and the next two theorems respectively show that it is still asymptotic preserving in the diffusive limit and locally decreases the total entropy. Theorem 3.1. If (Q) k are narrow bands, σ ν is constant inside each frequency group Q q and if the quadrature formula < • > is at least 2nd-order, then the model (3.1) -(3.7) is asymptotic preserving in the diffusive limit ie when σ → ∞ it degenerates into the equation:

∂ t (ρC v T + 4aT 4 ) -∇ 4acT 3 3σ R ∇ x T = 0, (3.9)
where σ R is Rosseland's mean value of the opacity given by (0.6).

Proof. Adding (3.4) and (3.5) together, and using the Eddington form of P q (2.10), the system (3.1) -(3.7) may be written as:

∀q = 1 . . . Q, ∂ t E q + ∇.F q = c σ q aθ 4 q -σ q E q , ∂ t F q + c 2 ∇(D q E q ) = -cσ q F q , < 1 c ∂ t (ΩI) + ∇(Ω ⊗ ΩI) > q = -< σΩI > q , 1 c ∂ t J k,l + Ω l ∇J k,l = σ k B k (T ) -J k,l ], k = 1 . . . K, l = 1 . . . L, ρC v ∂ t T = Q q=1
(σ q E q -σ q aθ 4 q ) + < cσ(J -B(T )) >.

To exhibit the asymptotic limit of this model in the diffusion regime, a Chapmann-Enskog-like expansion is used. First, the equations are scaled setting: t := εt, σ := εσ.

The system thus becomes:

ε 2 c ∂ tE q + ε∇F q = σq aθ 4 q -E q , q = 1 . . . Q, (3.10)
ε 2 c ∂ tF q + ε∇(D q E q ) = -σ q F q , (3.11)
ε 2 c ∂ tJ k,l + εΩ l ∇J k,l = σk B k (T ) -J k,l ], k = 1 . . . K, l = 1 . . . L, (3.12
)

ε 2 ρC v ∂ tT = Q q=1 (σ q E q -σ q aθ 4 q ) + < cσ(J -B) >. (3.13)
Then an asymptotic expansion of E, F and J is performed:

E = E 0 + εE 1 + ε 2 E 2 + . . . , F = F 0 + εF 1 + ε 2 F 2 + . . . , J = J 0 + εJ 1 + εJ 2 + . . . .
Let us now sum over all l the second equation. Since B is isotropic and the quadrature formula < • > exactly integrates 2nd order polynomials, we get:

4π c ∂ tB k (T ) -div 4π 3σ k ∇B k (T ) = - L l=1 σk J 2 k,l ω l .
Hence,

4π c ∂ tB k (T ) -div 4π 3σ k ∂ T B k (T )∇T = - L l=1 σk J 2 k,l ω l . (3.30)
Finally, summing all the equations (3.28) and (3.30) over all groups and bands:

∂ t(aT 4 ) - c 3 div < ∂ T B(T σ > ∇T = - Q q=1 σq E 2 q -< cσJ 2 >, (3.31)
This last equality is true thanks to the narrow band hypothesis that allows us to correctly handle the σk . The last step is performed by adding (3.25) to this last equation:

∂ t (ρC v T + aT 4 ) -∇ x 4acT 3 3σ R ∇T = 0.
∂ t < ȟ > c + ρC v U + ∇ < Ω ȟ >≤ 0, (3.32)
where U = U(T ) is a convex function and ȟ is given by:

ȟ = h if ν ∈ ∪Q q h d if ν ∈ ∪Q k (3.33)
Remark 3.3. In the above theorem α is Lagrange's multiplier that defines B, see (1.4) and (3.7).

Remark 3.4. The condition on (ω k ) k is not restrictive for the applications considered in this paper. Indeed, if the quadrature formula was exact, then α = 1/T . Only the quadrature errors may prevent α to be positive, decreasing and convex. Whenever only a few narrow bands are considered, it is easy to find a quadrature formula that allows to conserve these properties.

Proof. First, let us recall that the form of I given in (2.9) depends on a parameter per group denoted by α q which is a vector of IR 1+d (d is the space dimension). In this proof, we will write α q = (α 0 q , α x q , α y q ) , assuming that d = 2 (extensions to other dimensions are straightforward). We will also use the following notations:

A q = α x q 0 0 α y q , β q = α q .(1, Ω) .
It is to note that by construction of I, β q is nothing but:

β q = ∂ I h(I q ).
(3.34)

Let us now multiply (??) by α 0 q and (??) by A q :

< α 0 q ∂ t I > q + < cα 0 q ∇(ΩI) > q =< cσα 0 q (B(T ) -I) > q , < A q Ω∂ t I > q + < cA q ∇(Ω ⊗ ΩI) > q =< cσA q Ω(B(T ) -I) > q , Adding these two equations we get:

< β q ∂ t I > q + < cβ q Ω∇I > q =< cσβ q (B(T ) -I) > q .
Now, replacing β by its value according to (3.34), this last equation becomes:

∂ t < h(I) > q +∇ < cΩh(I) > q =< cσ∂ I h(I) B(T ) -I > q , (3.35)
which can be written as:

∂ t < h(I) > q +∇ < cΩh(I) > q =< cσ ∂ I h(B) -∂ I h(I) B(T ) -I > q (3.36) + < cσ∂ I h(B) B(T ) -I > q , (3.37)
and since h is a convex function and ∂ I h(B) = 1/T , we have:

∂ t < h(I) > q +∇ < cΩh(I) > q ≤< c σ T B(T ) -I > q . (3.38)
Similarly, multiplying (3.3) by ∂ I h d (J k,l ) leads to:

1 c ∂ I h d (J k,l )∂ t J k,l + ∂ I h d (J k,l )Ω l ∇ l J k,l = σ k ∂ I h d (J k,l ) B k -J k,l , (3.39)
then:

1 c ∂ t h d (J k,l ) + Ω l ∇ h d (J k,l ) = σ k ∂ I h d (J k,l ) B k -J k,l , (3.40)
and using the quadrature formula:

∂ t < h d (J k,l ) > + ∇< cΩ l h d (J k,l ) > = < cσ k ∂ I h d (J k,l ) B k -J k,l >. (3.41)
Proceeding as above and since ∂ I h d (B) = α we have: If α has the required behavior, all four terms inside the integral are positive, decreasing and convex functions of T , then U is a (negative, decreasing) convex function of T .

∂ t < h d (J k,l ) > + ∇< cΩ l h d (J k,l ) > ≤ < cσ k α B k -J k,l >. ( 3

Numerical Example

This numerical example intends to emphasize the behavior of the hybrid model. A 1 m 2 wide square domain is considered. Radiative equilibrium at T 0 = 1000 K is initially assumed. For t > 0, four "beams" are enforced at the NW, NE and SE corners as well as at the middle of the bottom of the domain (see figure 4.1 below). Each "beam" has the following profile: where T b = 6000 K, ν b = 1.10 14 Hz and ∆ν = 0.10 14 Hz. The simulation is made with two groups and one narrow band. Figure 4.1 shows the predicted radiative intensity integrated over all the frequencies and located at coordinates (0.5, 0.5) compared with the same computation carried out with a multigroup-M 1 model: As figure 4.1 shows, the hybrid model perfectly deals with this simple test-case: each beam is fully preserved with the correct amplitude and direction. Moreover, since there is only one narrow band, the computation is very cheap (almost immediate for a 100 × 100 mesh). On the other hand, the M 1 model, as any direction-integrated model, cannot properly deal with these four beams of the same amplitude converging. This is due to the fact that it has to cope with only one direction of anisotropy. This is a typical case where moments or flux-limited diffusion models fail to predict a realistic solution.

Even though this test is very simple, it is representative of what can be done for simulating a ICF-type problem where dozens of lasers coming from the same source (and hence having the energy concentrated inside the same frequencies) converge.

Theorem 3 . 2 .

 32 If (ω k ) k are chosen such that α(T ) is a positive, decreasing convex function of T , then the model (3.1) -(3.7) has a total entropy which is locally decreasing:

1 T

 1 .42) Now, let us respectively multiply(3.4) and (3.5) by 1/T and α (Lagrange's multiplier that defines B, see (1.4) ) to get:ρC v T ∂ t (T -T ) =< c σ T I -B >,(3.43)ρC v α∂ t T = < cσα J -B >. (3.44)Finally, summing (3.38) over q and adding (3.42), (3.43) and (3.44) leads to:∂ t < ȟ > c + ρC v U + ∇ < Ω ȟ >≤ 0, (3.45) U = -T ∂ T (T -T ) + α∂ T T . (3.46)

I

  ν (Ω) = B ν (T b ), if |ν -ν b | ≤ ∆ν, B ν (T 0 ), otherwise.

acknowledgement

  The author wish to thank Christophe Berthon for his enlightened pieces of advice.

Figure 4 . 1 .

 41 Figure 4.1. (above ) Directionnal energies given by the hybrid (l.) and M 1 (r.) models (below) illustration of the considered domain.

The use of this expansion in (3.10)- (3.13) and an identification of the coefficients lead to:

F 0 q = 0, (3.15)

(3.17)

The expression of F 1 q and J 1 k,l are respectively taken from (3.19) and (3.20):

The expressions (3.27) and (3.26) may be inserted in (3.22) and (3.24) to get:

Then E 0 q and J 0 k,l can be replaced by aθ 4 q andB k (T ) thanks to (3.14) and (3.16):


domain: sdv sde

text:
e-FLORA-sys, a website tool to evaluate the agronomical and environmental value of grasslands

e-FLORA-sys is a free website developed to provide a research and a diagnosis tool for grassland researchers and technicians. The system is based on databases describing grassland species, floristic composition of grasslands, agricultural practices, soils and climates, and vegetation associations. From this information, the system calculates numerous indexes to evaluate the agronomical and ecological value of grasslands. Users can freely record their own observations (botanical relevés, agricultural practices, ...), which are protected by a login device. E-FLORA-sys is currently used in a national program studying grassland agronomical and ecological value (CASDAR)

Diagnosis tool

The species table contains data on plants (currently 3000 species) found in most European grasslands : identification (translations and synonyms), agronomical value (potential production level, quality for cattle, sheep, goats and horses), patrimony value (rarity, inscription on red lists), reaction to abiotic factors (soil and climate) and agricultural practices (i.e. adaptation to frequent cutting or trampling) including Ellenberg indices, biological and ecological characteristics (aerial part, root and diaspore morphology, germination, reproduction and dispersion, life traits, Raunkier types).

The relevé table includes the floristic composition, and all the calculations, graphs and texts for an agronomical and an ecological diagnosis. Several methods of relevés are available. Each relevé is geolocated (Google maps®) and linked to a station which is itself linked to 3 tables describing the climate (temperature, precipitation, radiation), the soil (physical and chemical parameters) and the agricultural practices (fertilization, grazing, mowing, other practices). The associations table enables the users to link the relevés to recorded known types of grassland vegetation or European habitats. Relevés can be gathered in studies, in order to study grasslands of a particular region, to compare evolution in time, or to focus on a specific aspect (i.e. effect of nitrogen fertilization amount). Pictures of species and grasslands can be stored. The access to information is controlled by a system of user rights defined in the table of 'people and rights". A dictionary allows a full english and french translation of all technical and software command terms. A table is describing decision rules involved in agronomical and ecological interpretation, and in prediction of floristic composition (Flora-predict model) .

Botanical composition of the grassland

Forage value Ecological indexes

Biodiversity

Functional types

Introduction

  Aknowledgements : M. Pernot, B. Noirtin, M. Brunias, M. Vanschepdael, A. Chapon, F. Radet eFLORAsys database


Patrick Lambert

Guillaume Roch

Aurore Graillat

Jean Leonetti

Eric Martin

Anne Marie Ferréol

Alain Blanc

Jean Marc Audergon

Analysis of the resistance to PPV in a F1 progeny derived from the apricot cultivars 'Goldrich' and 'Moniqui'. Identification of the QTLs involved

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


Chiraz Jridi

Sarah Leon

Jean Luc Macia

Michel Yvon

Rémy Froissart

Specialization of a phytovirus to its environment: experimental evidence for a cost of adaptation

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

Specialization of a phytovirus to its environment: experimental evidence for a cost of adaptation

Chiraz JRIDI 1 , Sarah LEON 1 , Jean-Luc MACIA 1 , Michel YVON 1 and Remy FROISSART Studying patterns and genetic basis of adaptation is crucial to understand viral emergence and to anticipate suitable strategies of control. In this work, we propose to detect adaptive mutations in different environments (constant and variable) and to evaluate their phenotypic expression (measured by evaluation of fitness and virulence) in different environments.

We designed an experimental evolution of Cauliflower mosaic virus by transmitting viral populations from plant to plant in either two homogeneous environments (Arabidopsis thaliana or Nicotiana bigelovii) or a variable environment (alternation of both species). For each of these three treatments, 10 independent viral populations were evolved in parallel. After five and ten passages (approx. 150 and 300 viral generations, respectively), we sequenced complete genome length of consensus viral populations and look for mutations appearing at the same locus in more than two independent populations. These parallel mutations have a high probability to be a signature(s) of adaptive changes.

Genomic data revealed that parallel mutations appeared differentially depending on environmental contexts: (i) most of them were accumulating in constant environments; (ii) none were common to several environments. Evaluation of phenotypic effects of viral evolution in different environment gave the evidence for a cost of adaptation and suggest a positive correlation between virulence and virus accumulation. These genomic and phenotypic results will be discussed in light of evolution of specialist/generalist strategies.


Emmanuel Bertrand

Cyril Chevarin

Stéphane Portanguen

Pascal Tournayre

Saïd Abou

El Karam

Jean-Louis Berdagué

Alain Kondjoyan

Anne-Sophie Guillard

Identification of odour active compounds generated in a processed cheese-like model system

à 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.

The elaboration of food and especially their thermal treatment lead to major changes in food rheology, colour, nutritional value, microbiological stability and flavour properties. In the case of processed cheese, the first four points are relatively well known and taken into account for its production. However there is still a need to integrate the developpement of flavour in the multi-objectives optimisation strategies. This work aims at the identification of odour-active compounds which are influenced by the heating conditions. Cooking and analytical strategies enabled to get samples with very different odour active profiles. In order to understand the thermal effect on the odour active fraction, these compounds as well as other non odorant key intermediates will now be kinetically measured during various time and temperature treatments.

 Conclusion and prospects



  n.i. not identified A mass spectrum, retention index and odour identical with those of an authentic sample B mass spectrum, retention index and odour in agreement with literature data C tentative identification by mass spectra D odour and retention index in agreement with literature data, no mass spectrum found Overview of the Temperature profile inside the processed cheese sample at the end of the treatment Comparison between experimental and simulated temperature profile at the center of the cheese matrixOverview of GC-MS-Material and MethodsProcessed Cheese-like model formulation: Native Caseinate, skimmed milk powder, citric acid, emulsifying salts and natrium chloride were hydrated overnight in water at 4°C. Anhydrous milk fat was emulsified at 25°C for 10 min using a Janke and Kunkel T25 Ultra-Turrax (20500 rpm). The Cheese matrix was then stored at 4°C and cooked during its day of fabrication.Cooking process: a 10 g sample was cooked in an hermetically closed cooking cell designed to heat uniformly as quickly as possible the cheese matrix without any volatile compound leak. This cell is subjected to a succession of thermal treatments such as boiling water and/or superheated steam. Details concerning the pilot plant can be found in [1]. active compounds: An eight-way gas chromatograph-olfactometer (GC-MS-8WO) allows to detect a maximum of odorants. It consists of a chromatograph coupled to a divider that distributes the volatile component outflow evenly and synchronously to eight sniffing ports located in separated booths. The GC-MS-8WO analytical conditions are described in detail in [2]. An independent GC-O-MS device was used to phase olfactory and mass spectrometry signals in order to propose structures for odour-active compounds. GCxGC-Mstof analyses: A SPME Carboxen/PDMS, 75 μm fiber was equilibrated with the headspace of an 1g sample for 30 min at 30°C. A "LECO Pegasus IV" GC x GC-MStof instrument (LECO Corporation, St. Joseph, MI, USA) was used to identify the volatiles compounds. time-of-flight mass spectrometry (GCxGC-MStof) Euro Food Chem XV -July 5-8, 2009 -Copenhagen

1

  Institut National de la Recherche Agronomique (INRA), Unité de Recherche sur la Qualité des Produits Animaux (UR 370), F -63122 Saint Genès Champanelle, France; 2 Fromageries Bel D.R.A.G., 7 Bd de l'Industrie , F -41100 Vendôme cedex, France; Evaluation of the thermal treatment : Evaluation of odour active compounds:

  Odour active compounds detected by GC-O and identified by GCxGC-MStof Both experimental measurements and numerical simulations demonstrate a good temperature reproductibility at the end of the treatment (less than +/-5°C) Aromagrams obtained with the 8W-GCO device for the raw sample (left) and the heated one (right) *email: emmanuel.bertrand@clermont.inra.fr (1)Kondjoyan, A. and Portanguen, S. (2008). Prediction of surface and "under surface" temperatures on poultry muscles and poultry skins subjected to jets of superheated steam. Food Research International, 41(1):16-30. (2) Berdagué, J.L.; Tournayre, P.; Cambou, S. 2007.Novel multi-gas chromatographyolfactometry device and software for the identification of odour-active compounds. Journal of Chromatography A, 1146 (1), 85-92


domain: math info

text:
Pierre Biver

Denis Allard

Dimitri D'or

COMBINING DIRICHLET DISTRIBUTIONS FOR MODELING LITHO-TYPE WITH SOFT PROBABILTIES FROM SEISMIC ATTRIBUTES AND EXTERNAL INFORMATION

We consider the problem of deriving consistent 3D proportion trends from seismic attributes. From well data, i.e. where seismic attributes and litho-types are known, we derive conditional probabilities with a kernel smoothing algorithm applied for each litho-type from which we define a probability model with the Dirichlet distribution. In this framework, we show how these probabilities can be combined with external geological information. Some implementation issues are discussed and the approach is illustrated on a real field case.

INTRODUCTION

In many applications of earth modelling, it is necessary to perform a 3D lithology model because petrophysical properties of the medium are strongly dependent on litho-type. In stochastic models, lithology assignment is controlled by the probability of occurrence (or proportion) assigned as input. Due to conceptual geology and indirect measures (seismic campaign,...), a non stationary 3D proportion trend is often needed. This basic information is needed for all geostatistical algorithms addressing lithology simulation (Sequential Indicator Simulation, Sequential Categorical Simulation (D'Or et al., 2008), Truncated Gaussian Simulation, Boolean models, and even Multiple Point Statistics Simulation); it is not easy to obtain it and a rigorous methodology has to be used to build these 3D trends.

In this paper, we suggest here below a formalism to incorporate indirect measures in proportion evaluation. The outline of the paper is the following. We first present the methodology: from well data (where a couple of seismic attributes and the litho-types are known) conditional probabilities for each litho-type are estimated using a kernel smoothing algorithm. If necessary, the global marginal distribution of seismic attributes from the exhaustive grid can also be used as weighting Comment citer ce document : Allard, D. (2008). Litho-type modeling using soft-probabilities from seismic attributes and other sources of information. In: Geostats 2008. 2 volumes (p. 319-328). Presented at 8. International Geostatistics Congress, Santiago, CHL (2008-12-01 -2008-12-05). Santiago, CHL : Universidad de Chile.

P. BIVER, D. ALLARD AND D. D'OR correction. Then, these conditional probabilities are used as parameters for a soft probability model using a Dirichlet distribution. It is shown that the Equivalent Number of Independent Data (ENID) plays an important role, and we present how it can be comptued. In this framework, we show how this soft probability model can be combined with external geological information and we propose indices to assess the uncertainty related to the probability model. A real field case is then analysed using our approach, which illustrates that it allows us to obtain at the same time a reliable trend surface and a consistent soft-probabilty model with a measure of its associated uncertainty.

WEIGHTED KERNEL REGRESSION

Basic kernel regression algorithm

Using the seismic information and the lithotypes measured at the wells, each well data can be reported on the cross-plot of seismic variables and labeled with the litho-type. Now, we would like to obtain the probability distribution of the litho-types for every possible pair of seismic values on the cross-plot area. We thus have to perform some kind of regression on the litho-types probabilities. To achieve this, we propose to use the kernel regression approach [START_REF] Nadaraya | Nonparametric Estimates of Probability densities and Regression Curves[END_REF][START_REF] Watson | Smooth regression analysis[END_REF].

Let us consider the case of bivariate data (e.g. two seismic attributes), where s 1 , .

Y (s) = f (s) + ε(s) (1)
where f (s) is the deterministic trend we wish to estimate, and ε(s) is a random residual. The kernel methods propose a linear estimator for f (s):

f (s) = n ∑ i=1 w i Y (s i ) (2)
for which the weights w i are functions of the distances between s and the data locations s i , and where the sum of the weights is equal to one. The weights are written:

w i = K h (s -s i ) ∑ n i=1 K h (s -s i ) (3)
where K h (u) is a compact writing for In our situation, the random variables are the indicator functions Y k (s) for each one of the K litho-type, i.e. Y k (s i ) = 1 if litho-type k is observed at location s i and 0 otherwise. Eq. ( 3) automatically ensures that all the estimated probabilities belong to [0; 1] and sum up to one at each location, yielding valid probability distributions of litho-types at location s. As primary result, the kernel regression procedure yields a probability map per litho facies over the whole cross-plot area.

K h (u) = 1 h x h y K u x h x ; u y h y (4) in which K(•,

Weighting procedure

The kernel regression procedure may be easily modified to become weighted. The regression weights w i are then modified from Eq. ( 3) to incorporate some weights p(s i ):

w i = K h (s -s i )p(s i ) ∑ n i=1 K h (s -s i )p(s i ) (5)
Two kinds of weights p(s i ) may be used: (1) declustering weights, and (2) bias correction weights.

Declustering weights

Declustering weights are obtained from a declustering of the data performed on the reservoir area and allow to take into account a possible redundancy of information brought by data that are clustered on this area.

Bias correction weights

As shown above, kernel regression is performed using well data only. In some situations, those data may not be perfectly representative of the true bivariate distribution of seismic variables. This distribution may be estimated more accurately using the whole seismic data avalaible on the reservoir grid. The bias is then defined as the difference between the probability density functions estimated from the whole seismic data set on one hand, denoted by fIR (s), and from the well data on the other hand, denoted by fIW (s). The kernel smoothing approach can be used to estimate both densities with:

f (s) = 1 n n ∑ i=1 K h (s -s i ) (6)
where K h (ss i ) is a kernel as defined in Eq. ( 4). For each data set, the data are used all together to estimated the density, whatever their litho-type.

Bias correction weights may then be computed as:

p(s i ) = 1 + ∆ I (s i ) n (7) with ∆ I (s i ) = a fIR (s i ) -fIW (s i )
where a is the area of a grid cell on the cross-plot. Doing so,

∆ I (s i ) ∈ [-1; 1]. Eq. (7) ensures that: i) p(s i ) → 0 when ∆ I (s i ) → -1; ii) p(s i ) → 1 when ∆ I (s i ) → 1; iii) p(s i ) = 1/n if ∆ I (I i ) = 0
(with n, the number of pairs of seismic variables measured at the wells); and iv) ∑ n p(s i ) = 1.

Uncertainty indices

In diagnostics, a first interesting result is the map of the most likely litho-type on the cross-plot area. This map is easy to compute from the litho-types probability maps. It may be accompanied by two maps of measures of uncertainty on the estimation. A first type of uncertainty arises from the estimated probability distribution itself. It is related to the uncertainty on the choice of the most likely facies. To measure this uncertainty, we can compute a quality index derived from the Gini index, which is defined as [START_REF] Gini | Variabilita e mutabilita: Contributo allo studio delle distribuzioni e delle relazioni statistice[END_REF]:

G(s) = ∑ i = j p i (s)p j (s) = 1 -∑ i p i (s) 2 (8)
The Gini index G is 0 if one of the litho-types has a probability equal to 1 (no uncertainty), and maximum (equal to 1 -(1/K)) when all the litho-types are equilikely (maximum uncertainty). The maximum goes towards 1 when the number of litho-types increases. In order to allow for comparisons, an uncertainty index belonging to [0; 1], equal to 0 when the uncertainty is maximum and to 1 when there is no uncertainty, can be computed from the Gini index as follows:

Q(s) = 1 - G(s) G max with G max = 1 - 1 K (9)
The second uncertainty measure is related to the uncertainty on the estimation of the individual proportions. It can be quantified by the Equivalent Number of Independent Data (ENID), ñ:

ñ = 1 ∑ w 2 i ( 10
)
with the weights w i obtained using Eq. ( 5) in which all the data points are considered, wathever their litho-type. The ENID map allows to identify areas that less informed by the data than others.

DIRICHLET DISTRIBUTIONS

Description

Using a Bayesian formalism, [START_REF] Haas | Uncertainties in facies proportion estimation i. theoretical framework: the dirichlet distribution[END_REF] and [START_REF] Biver | Uncertainties in facies proportion estimation ii. application to geostatistical simulation of facies and assessment of volumetric uncertainties[END_REF] have shown that the uncertainty on the vector of litho-types probabilities may be modeled using a Dirichlet distribution. The Dirichlet distribution with K parameters α k is defined as

f n (z) = f n (z 1 , . . . , z K-1 ) = Γ(θ ) Γ(α 1 ) × . . . × Γ(α K ) z α 1 -1 1 × . . . × z α K -1 K with K ∑ k=1 z k = 1 and θ = K ∑ i=1 α i = n + K (11)
where z k is the proportion (probability) of litho-type k, Γ(q) is the gamma function for the positive variable q. It is easy to see that when K = 2, the Dirichlet distribution is in fact a Beta distribution.

Comment citer ce document : Allard, D. ( 2008). Litho-type modeling using soft-probabilities from seismic attributes and other sources of information. In: Geostats 2008. 2 volumes (p. 319-328). Presented at 8. International Geostatistics Congress, Santiago, CHL (2008-12-01 -2008-12-05). Santiago, CHL : Universidad de Chile.

MODELING LITHO-TYPES WITH SOFT PROBABILITIES"

The α k parameters can be computed as

α k = ñz k + 1 (12)
where ñ is the ENID, zk is the estimated proportion. If necessary, Eq. ( 12) can be made non stationary by indexing its components with location s i .

Using this formalism, vectors of proportions may be randomly generated by drawing in the Dirichlet distribution with parameters α k . Note that the ENID plays a crucial role as it is the parameter that controls the uncertainty on the distribution of proportions. A large ENID induces low uncertainty and vice-versa.

Combination

In some situations, additional information may be available. It can be either another independent data set or an expert knowledge. In both cases, it is interesting to incorporate it in the analysis. It is easy to do it in the Dirichlet framework.

Additional data set

Let z k1 and z k2 respectivly denote the proportions of the two data sets, and ñ1 and ñ2 , the associated ENIDs. The two Dirichlet distributions may be combined by computing a combined ENID and a combined vector of proportions with ñc = ñ1 + ñ2 and

z kc = ñ1 z k1 + ñ2 z k2 ñ1 + ñ2 (13)
Those combined parameters are then plugged into Eq. ( 12) to obtain the parameters of the combined distribution.

Expert knowledge

The combined ENID and vector of proportions are here obtained using ñkc = ñ(1 + r) and

z kc = z k + r.y k 1 + r ( 14
)
where r is the relative confidence in the expert judgement (varying between 0 and +∞) and y k is the proportion for litho-type k provided by the expert. Again, those two parameters are plugged into Eq. ( 12) to obtain the parameters of the combined distribution.

CASE STUDY

Data set description

To test the methodology described in the previous paragraphs, a turbiditic reservoir has been chosen. attributes seem to be informative for lithology (Figure 1b); grid sampling at well data location can be performed and the crossplot PR residual versus IP residual at wells can be labeled by litho-types (Figure 1a); a first look at this crossplot confirms the prediction potential of seismic attributes for litho-type modelling.

Crossplot smoothing results

To provide a 3D proportion trend based on seismic attributes available exhaustively on the reservoir grid, it is necessary to associate probabilities of litho-types to each value of seismic attributes and therefore to perform the smoothing of the crossplot.

Comment citer ce document : Allard, D. ( 2008). Litho-type modeling using soft-probabilities from seismic attributes and other sources of information. In: Geostats 2008. 2 volumes (p. 319-328). Presented at 8. International Geostatistics Congress, Santiago, CHL (2008-12-01 -2008-12-05). Santiago, CHL : Universidad de Chile.

MODELING LITHO-TYPES WITH SOFT PROBABILITIES"

Several parameters have to be chosen for this operation: kernel type, data weights (same weights for all data in our case as there are neither bias nor clusters), kernel size to control smoothing. Impact of kernel size is illustrated on Figure 2 for massive sand probability of occurrence: a 100*100 kernel size leads to a regression surface containing irrelevant variations (interpreted as non significant noise), while a 200*200 kernel size reflects the trend surface, without erasing significant variations. The latter kernel size is thus chosen. Once this crossplot smoothing is validated, it is possible to assign the information to reservoir grid since now seismic attributes are define over the entire domain. This operation is illustrated on Figure 3 (first part) for the probability of occurrence of massive sands and for the most likely litho-type; the main channel shape can be distinguished on these maps.

Some additional diagnostics can be performed, see Figure 3 (second part):

-from the computation of the ENID, a region where proportions are known with the best precision can be identified; it is in an area of the crossplot dominated by massive clay lithology; -from the computation of the Gini index, a region where maximum mixing is occurring can be identified; the most likely facies is not relevant in this region;

Combination with geological information

The seismic campaign has been previously interpreted qualitatively by geologists; they have picked polygons corresponding to major geological events and they have computed vertical proportion curves (from well data and conceptual geology) for each interpreted polygon. This information could be used to correct 3D proportion

VIII International Geostatistics Congress

Comment citer ce document : Allard, D. ( 2008). Litho-type modeling using soft-probabilities from seismic attributes and other sources of information. In: Geostats 2008. 2 volumes (p. 319-328). Presented at 8. International Geostatistics Congress, Santiago, CHL (2008-12-01 -2008-12-05). Santiago, CHL : Universidad de Chile.

P. BIVER, D. ALLARD AND D. D'OR trend established from seismic attributes; in order to achieve this goal, an equivalent number of independent data has been assigned to this geological trend (based on well data number and confidence in the geological interpretation).

If an equivalent number of independent data is assigned on both seismic trend and geological trend, it is possible to calculate the resulting Dirichlet distribution mixing both informations with their relevant weights. This operation is illustrated on Figure 4 for the probability of occurrence of massive sands and for the most likely litho-type. In the region delimited by the ellipse, the ENID of the geological interpretation is larger than the ENID of the seismic data analysis. This results in an improved visibility of the channel meander and its lateral overbank on the combined proportion trend. In other regions, the seismic data analysis has the largest ENID and thus is better represented on the combined proportion trend.

CONCLUSIONS

A rigorous methodology has been established to build 3D proportion cubes including seismic and geological information with their relevant weight (associated to uncertainties). It has been tested on a real field case with satisfactory results. However some points need further investigations: i) testing the efficiency of bias correction on other field cases with less well data sampling, and ii) comparing the combination of probabilities using Dirichlet formalism with other theories for combining information (e.g. Tau model from [START_REF] Journel | Combining knowledge from diverse sources: an alternative to traditional data independence hypotheses[END_REF], or minimum variance combination).

Comment citer ce document : Allard, D. ( 2008). Litho-type modeling using soft-probabilities from seismic attributes and other sources of information. In: Geostats 2008. 2 volumes (p. 319-328). Presented at 8. International Geostatistics Congress, Santiago, CHL (2008-12-01 -2008-12-05). Santiago, CHL : Universidad de Chile.

P

  . BIVER, D. ALLARD AND D. D'OR

Figure 1 :

 1 Figure 1: Data set used to test the methodology. a) Well data: on the basis of Facies log data interpretation at 8 wells, and seismic attributes extraction at these wells location, a crossplot of residual Poisson ratio versus residual acoustic impedance of compression waves is labeled with litho-type. b) Seismic data: on the reservoir grid, two seismic attributes are available in each cell: residual acoustic impedance of compression waves and residual Poisson ratio (interpreted from the seismic campaign).

Figure 2 :

 2 Figure 2: Influence of the Kernel size on the smoothing results: the probability of occurrence for each litho-type is evaluated from experimental crossplot by kernel smoothing; with a small size of the kernel (100*100) unrealistic fluctuations are occurring on probability map; these fluctuations are filtered out with a larger kernel size (200*200).

Figure 3 :

 3 Figure 3: Smoothing results (crossplot and reservoir grid) and diagnosis (crossplot): (a) probability of occurrence of massive sands in the bivariate space of seismic attributes; (b) most probable litho-type in the same bivariate space; (c) and (d) same variables directly on the reservoir grid; (e) equivalent number of independent data in the bivariate space of seismic attributes (ellipse indicates the most representative area, see also (a)); (f) Gini index in the same bivariate space (ellipse indicates the area with an important mixing between litho-types, see also (b)).

  It is sampled with a set of wells reasonably representative (no sampling bias problem); the well data have been interpreted in five different litho-types. The seismic campaign has an intermediate quality but seismic

VIII International Geostatistics Congress

Comment citer ce document : Allard, D. (

2008

). Litho-type modeling using soft-probabilities from seismic attributes and other sources of information. In: Geostats 2008. 2 volumes (p. 319-328). Presented at 8. International Geostatistics Congress, Santiago, CHL

(2008-12-01 -2008-12-05)

. Santiago, CHL : Universidad de Chile.

GEOSTATS 2008, Santiago, Chile

ACKNOWLEDGEMENTS

The authors would like to acknowledge R. Froidevaux (FSS Consultants) for positive discussions and debates on this subject and Total SA for sponsoring this research and allowing this publication.


text:
Ali Chahine

Sylvain Dupont

Yves Brunet

Carole A Sinfort

A large-eddy simulation model for pesticide dispersal during spray application

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