E Caffau

R Faraggiana

H.-G Ludwig

P Bonifacio

M Steffen

The solar photospheric abundance of zirconium

Keywords: abundances -Stars, abundances -Hydrodynamics -Line, formation

Zirconium (Zr), together with strontium and yttrium, is an important element in the understanding of the Galactic nucleosynthesis. In fact, the triad Sr-Y-Zr constitutes the first peak of s-process elements. Despite its general relevance not many studies of the solar abundance of Zr were conducted. We derive the zirconium abundance in the solar photosphere with the same CO5BOLD hydrodynamical model of the solar atmosphere that we previously used to investigate the abundances of C-N-O. We review the zirconium lines available in the observed solar spectra and select a sample of lines to determine the zirconium abundance, considering lines of neutral and singly ionised zirconium. We apply different line profile fitting strategies for a reliable analysis of Zr lines that are blended by lines of other elements. The abundance obtained from lines of neutral zirconium is very uncertain because these lines are commonly blended and weak in the solar spectrum. However, we believe that some lines of ionised zirconium are reliable abundance indicators. Restricting the set to Zr II lines, from the CO5BOLD 3D model atmosphere we derive A(Zr)=2.62 ± 0.06, where the quoted error is the RMS line-to-line scatter.

Introduction

The photospheric solar abundances of some elements have been studied more often than others; this is mainly due to the importance of the elements to explain nucleosynthesic processes. Moreover, the difficulty of extracting suitable lines with accurately known transition probabilities in the visible range of the solar spectrum can also explain why some elements have been studied less extensively.

For the triad Sr-Y-Zr, there exist only a few detailed studies of their solar abundances (an ADS 1 search yields 4 papers for Sr, 5 for Y, and 10 for Zr) in spite of their importance in Galactic chemical evolution. Sr-Y-Zr, Ba-La, and Pb are located at three abundance peaks of the s-processes producing the enrichment of these elements in the Galaxy. Knowledge of the present Sr-Y-Zr abundances in stars of different metallicity and age is required to understand the complicated nucleosynthesis of these elements (for details see Travaglio et al. 2004).

Zirconium is present in the solar spectrum with lines of Zr I and Zr II. The dominant species is Zr II. According to our 3D model, and in agreement with our 1D reference model, about 99% of zirconium is singly ionised in the solar photosphere. Departure from local thermodynam- 1 The Astrophysics Data System, ads.harvard.edu ical equilibrium (LTE) probably affects Zr I lines, producing a too low Zr abundance under the assumption of LTE. In fact NLTE effects, even on weak unsaturated Zr I lines, have been found by Brown et al. (1983) for G and K giants. However, both Biemont et al. (1981) and Bogdanovich et al. (1996) analysed in LTE Zr I lines (34 and 21 lines, respectively) and Zr II lines (24 and 15 lines, respectively) in the solar photosphere, and found an excellent agreement between the abundances derived from both ionisation stages. Very recently Velichko et al. (2010) performed a NLTE analysis of zirconium in the case of the Sun and late type stars. According to their computations the NLTE abundance is larger than the LTE one, by up to 0.03 dex for Zr II and by 0.29 dex for Zr I.

Owing to the small number of Zr I and Zr II lines in metal poor stars (Gratton & Sneden 1994), it is important to analyse both of them in the solar photosphere to derive the solar abundance from both ionisation stages and to assess the agreement between the derived abundances. In spectra of cool stars it is easier to observe Zr I lines. For example, Goswami & Aoki (2010) realised that none of the Zr II lines were usable in their analysis of the cool Pop. II CH star HD 209621, and the Zr abundance is derived from the only Zr I line in their spectrum at 613.457 nm. A similar situation had been encountered by Vanture & Wallerstein (2002) in their study of the Zr/Ti abundance ratio in cool S stars, where only Zr I lines in the red part of the spectrum can be used.

Lines and atomic data

According to Malcheva et al. (2006), zirconium has five stable isotopes. Four of them ( 90 Zr, 91 Zr, 92 Zr, and 94 Zr) are produced by the s-process; the fifth, 96 Zr, is produced in the r-process. Zirconium is a very refractory element, it is difficult to be vaporised by conventional thermal means and, consequently, has been relatively little investigated in the laboratory.

In Table 1 and Table 2 the Zr I and Zr II lines, with the log gf values available in the literature are collected. Biemont et al. (1981) selected lines that lie at λ < 800 nm. They derived log gf using the technique developed by Hannaford & Lowe (1981), suitable for highly refractory elements like Zr, to determine lifetimes of 34 levels of Zr I and 20 levels of Zr II. From these measurements, coupled with the measurements of branching ratios, they derived the log gf values of 38 Zr I and 31 Zr II lines. Bogdanovich et al. (1996) computed and used log gf values for 21 Zr I lines. They give new log gf values for 15 Zr II lines among those selected by Biemont et al. (1981). Their log gf are given in Table 2. Sikström et al. (1999) derived the abundance of zirconium in HgMn star χ Lupi, finding a disagreement when using Zr II or Zr III lines. They measured the f-values for several Zr II lines in the UV at λ < 300 nm. Vanture & Wallerstein (2002) extended the search for Zr I lines in the near-IR, at wavelengths longer than 800 nm, a region important in the study of cool stars, because this is the region were cool stars emit most of their flux. The complete sample of these near-infrared lines is also given in Table 1. Because absorption bands of molecules are weak in the warmer S stars, these near-IR zirconium lines are better abundance indicators than the Zr II lines in the blue part of the spectrum. These lines are, however, not necessarily good for the Sun.

We have also checked the 13 Zr I lines identified by Swensson et al. (1970) in the Delbouille et al. (1981) atlas, but 12 of them appear severely blended. Only the line at 784.9 nm is retained in Table 1, as it was by Biemont et al. (1981) and Vanture & Wallerstein (2002).

In the present analysis we adopt the log gf values of Biemont et al. (1981) for all Zr I lines.

According to the NIST database, Zr II lines lie only in the near UV-blue (241.941-535.035nm); only two weak lines (at 667.801 and 678.715 nm) are present in NIST with λ >540 nm.

The line identifications by Moore et al. (1966) on the Utrecht solar atlas have been used by Biemont et al. (1981) to select 31 Zr II lines, seven of which were later discarded because the derived abundance exceeded the mean by more than 3σ, so they are likely blended with other unknown species. Gratton & Sneden (1994) used the same line list as Biemont et al. (1981). Ljung et al. (2006) derived new oscillator strengths for 263 Zr II lines and studied 7 lines, the lines that they judged to be the best and unperturbed in the photospheric spectrum, to derive the solar Zr abundance based on both 1D and 3D model atmospheres.

We extracted from the 243 lines by Malcheva et al. (2006), those in common with Biemont et al. (1981); for these lines the best values of log gf , according to Malcheva et al. (2006), are the same as in Biemont et al. (1981). We compared the new log gf values determined by Ljung et al. (2006) with those by Biemont et al. (1981). There is a good agreement for most of the lines, with a few exceptions (see Table 2). We report all these values in Table 2.

In the present analysis we used the log gf values of Ljung et al. (2006) for all the Zr II lines.

Solar Zr abundance in the literature

Several analyses of the solar Zr abundance, made before 1980, used [START_REF] Corliss | Experimental transition probabilities for spectral lines of seventy elements; derived from the NBS Tables of spectral-line intensities[END_REF] log gf [START_REF] Aller | Advances in Astronomy and Astrophysics[END_REF][START_REF] Wallerstein | A preliminary analysis of the abundances of the rare Earths in the Sun[END_REF][START_REF] Grevesse | Origin and distribution of the elements[END_REF] or older ones (Goldberg et al. 1960) or have been made to derive a temperature correction for the [START_REF] Corliss | Experimental transition probabilities for spectral lines of seventy elements; derived from the NBS Tables of spectral-line intensities[END_REF] data (Allen 1976). The use of these transition probabilities, which are known to contain errors depending on excitation and temperature, coupled with the use of old solar spectral atlases, explains discrepancies of these Zr abundance determinations based either on Zr I and/or Zr II.

We concentrate on the more recent determinations of the photospheric solar Zr abundance, that we summarise below.

1. Biemont et al. (1981) analysed the Delbouille et al.

( 1973) atlas. The equivalent widths (EWs) were independently measured by E. Bièmont and N. Grevesse, and the results examined and discussed to produce the final list. They used the Holweger-Müller solar model (Holweger 1967;Holweger & Müller 1974) with a microturbulence, ξ, of 0.8 km s -1 . They note that their analysis is independent of line-broadening parameters since most of the lines are faint or very faint. The log gf are from their own measurements. They find that the abundance from Zr I lines is independent of ξ and that the abundance from Zr II lines is insensitive to the model. These authors rejected 4 Zr I lines and 7 Zr II lines from their selected sample in their final solar abundance analysis because they imply a Zr abundance which is more than 3σ higher than the average. These lines appear after the horizontal line in Tables 1 and 2. They conclude that Zr II lines are better abundance indicators. They obtain a zirconium abundance of A(Zr)=2.57 ± 0.07 from Zr I (34 lines) and A(Zr)=2.56 ± 0.07 from Zr II (24 lines). 2. Gratton & Sneden (1994) studied the Zr behaviour in metal-poor stars. They used as model atmosphere a grid from [START_REF] Bell | [END_REF] for giants, and a similar grid for dwarfs provided by Bell. They also analysed the solar spectrum, to have a reference A(Zr)⊙. For this purpose they used the same 24 lines as Biemont et al. (1981), the Holweger-Müller model with microturbulence of 1.5 km s -1 , and the solar abundances from [START_REF] Anders | [END_REF]. They remark that in their metal poor stars Zr abundances from Zr I are, on average, lower than those from Zr II lines, a result similar to that found by Brown et al. (1983) for Pop. I stars. For the analysis of the metal poor stars they used 7 Zr II lines and 3 of them (407.1, 414.9, and 416.1 nm) are among those discarded by Biemont et al. (1981), suggesting that these are no longer significantly blended in metal poor stars. They obtain for the solar Zr abundance 2.59±0.04 from Zr I (34 lines) σ=0.22 dex, 2.53±0.03 from Zr II (24 lines) σ=0.14 dex, respectively. 3. Bogdanovich et al. (1996) computed theoretical log gf , and used the Holweger-Müller model with ξ=0.8, and a model computed with a code by Kipper et al. (1981) which is a modified version of the ATLAS 5 code (Ku-rucz 1970). They used the damping constants modified by [START_REF] Galdikas | [END_REF], and the EW measurements of Biemont et al. (1981). The results are: A(Zr)=2.60±0.07 from Zr I (21 lines) and A(Zr)=2.61±0.11 from Zr II (15 lines). They adopted A(Zr)= 2.60 ± 0.06. 4. Ljung et al. (2006) (347.9, 350.5, 355.1, 405.0, 420.8, 425.8, 444.2, 449.6, and 511.2 nm).

Analysing the Kurucz (2005a) solar spectrum with a MAFAGS model atmosphere (5780 K/4.44/0.0) and a microturbulence of 0.9km s -1 they derive a LTE abundance of 2.33 and 2.61 from Zr I and Zr II lines, respectively. They constructed a Zr model atom and derived the NLTE abundance with different assumptions about the rates of collisions with H-atoms. In their scenario, the NLTE abundance is always higher than the LTE abundance varying in a range from 0.21 to 0.32 dex for Zr I and from 0.01 to 0.08 dex for Zr II, depending on the adopted cross sections for collision with H-atoms. Their best estimate of the NLTE abundance (correction) from Zr I and Zr II lines is 2.62 (0.29) and 2.64 (0.03) dex, respectively. Their adopted Zr abundance is 2.63 ± 0.07.

The abundance of Zr in the solar system is derived from the analysis of meteorites. Zirconium is not a volatile element, and its abundance in meteorites should agree with the one in the solar photosphere. To compare the abundances derived from the meteorites, on the cosmochemical abundance scale relative to 10 6 Si atoms, with the ones derived form the solar photosphere, on the astronomical scale of 10 12 H-atoms, a coupling factor must be derived. [START_REF] Lodders | [END_REF] and [START_REF] Lodders | Abundances of the elements in the solar system[END_REF] to link the two scales introduce an average factor by looking at a sample of refractory elements well-determined in the photosphere and at the same elements in meteorites. In this way the factor is not affected by the individual uncertainty and, if the abundance of one element derived from the photosphere changes, the conversion factor is not much affected.

The meteoritic Zr abundance is A(Zr)=2.57±0.04, [START_REF] Lodders | Abundances of the elements in the solar system[END_REF]. Other Zr meteoritic values that can be found in the literature are: 2.61 ± 0.03 according to [START_REF] Anders | [END_REF], 2.61 ± 0.02 according to Grevesse & Sauval (1998), 2.60±0.02 according to [START_REF] Lodders | [END_REF]. The value of 2.53 ± 0.04 given in Asplund et al. (2009) is based on the same data from [START_REF] Lodders | Abundances of the elements in the solar system[END_REF], but they give a lower value for Zr on the astronomical scale because they use a lower coupling factor for the cosmochemical and astronomical scales. In the same way the value of 2.57 ± 0.02 in [START_REF] Asplund | Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis[END_REF] and Grevesse et al. (2007) is based on data from [START_REF] Lodders | [END_REF], but their value is lower because again a different scaling factor (which is 0.03 log units smaller than the one recommended in Lodders ( 2003)) was used.

Model atmospheres

We base our abundance analysis on the same model atmospheres that we used in the previous analysis of solar abundance determinations, summarised in Caffau et al. (2010). We rely on a time-dependent, 3D hydrodynamical model atmosphere, computed with the CO 5 BOLD code; see Freytag, [START_REF] Freytag | [END_REF] and [START_REF] Freytag | CO5BOLD User Manual[END_REF] for details. This model has a box size of 5.6 × 5.6 × 2.27 Mm 3 , resolved by 140 × 140 × 150 grid points. Its range in Rosseland optical depth is -6.7 < log τ Ross < 5.5. The complete time series is formed by 90 snapshots covering 1.2 h of solar time. For the spectral synthesis computations we selected 19 representative snapshots out of the complete series of 90 snapshots.

To compute 3D-corrections for the zirconium abundance, we make use of two 1D reference atmospheres: the 1D model obtained by horizontally averaging each 3D snapshot over surfaces of equal (Rosseland) optical depth (henceforth 3D model), and the hydrostatic 1D mixing-length model computed with the LHD code, that employs the same micro-physics and radiative transfer scheme as the CO 5 BOLD code, (henceforth 1D LHD model). We define 3D-corrections as in Caffau et al. (2010), ∆ gran = A(Zr) 3D -A(Zr) 3D , to isolate the effects of horizontal fluctuations (granulation effects), and ∆ LHD = A(Zr) 3D -A(Zr) LHD , to measure the total 3D effect. For details see Caffau et al. (2010) and Caffau & Ludwig (2007).

We also considered the semi-empirical Holweger-Müller solar model (Holweger 1967;Holweger & Müller 1974, hereafter HM) for comparison.

Observed spectra

We analysed the same four solar spectral atlases, publicly available (two disc-centre and two disc-integrated atlases), that we used in all our previous works on solar abundances, such as in Caffau et al. (2010). These high resolution, high signal-to-noise ratio (S/N) spectra, are: 1. the disc-integrated spectrum of Kurucz (2005a), based on fifty solar FTS scans, observed by J. Brault and L. Testerman at Kitt Peak between 1981 and 1984;

2. the FTS atlas of Neckel & Labs (1984), observed at Kitt Peak in the 1980ies, providing both disc-centre and discintegrated data; 3. the disc-centre atlas of Delbouille et al. (1973), observed from the Jungfraujoch.

Abundance determination methods

In principle, the best way to derive the abundance of an element from a spectral line is by line profile fitting. Once the atomic data are known, this technique allows to take into account the strength of the line and its shape at the same time.

The limiting factor is that the observed solar spectra are of such a good quality that the synthetic line profile is often not realistic enough to reproduce the shape of the line in the observed spectrum. Synthetic line profiles derived from 1D model atmospheres do not take into account the line asymmetry induced by convection, they are symmetric. In many cases NLTE effects modify the line profile (see Asplund et al. 2004). In case the line is not clean, but some contaminating lines are present in the range, the line profile fitting should take into account these blending components. Good atomic data for these blending lines are then necessary, and, as for the line of interest, granulation and NLTE effects can be limiting factors.

To fit a line blended with other components, one has to optimise the agreement between synthetic and observed spectra also for the blending lines, allowing, in the line profile fitting process, the abundance of the blending lines to change as well. Even though the 3D spectrum synthesis is still computationally demanding, we have attempted to derive A(Zr) by properly taking into account also the blending lines whenever possible. In the fitting procedure, the Zr abundance and the abundances of the blending lines are adjusted independently until the best match of the observed spectrum is achieved (in the following referred to as method 'A').

We also tested a simplified fitting approach (in the following method 'B'): we synthesised one typical 3D profile of each blending components, in addition to a grid with different abundances for the Zr line of interest. Then, in the fitting procedure, a profile is interpolated in the grid for Zr, while the profiles (more precisely the relative line depressions) of the blending lines are simply scaled, and then the blending lines and the Zr line are added linearly to obtain the composite line profile. In the fitting procedure, the Zr abundance and the scaling factors (and shifts) of the blending lines are adjusted until the best agreement with the observed spectrum is found. An advantage of this fitting procedure is that it works even if the atomic data of the blending components are poorly known. Note, however, that this procedure is theoretically justified only if all the lines are on the linear part of the curve of growth. If the lines are partly saturated, the line strengths of all components are underestimated, and so is the derived Zr abundance. Nevertheless, method 'B' was also applied to partly saturated lines in order to understand its limitations.

Methods 'A' and 'B' provide the Zr abundance without the need to measure an equivalent width. However, we can formally derive the equivalent width of the Zr component from the synthetic line profile computed for Zr only (ignoring all blends) with the Zr abundance that provides the best fit of the observed spectrum. For method 'B', this is simply the equivalent width of the Zr line that gives the best agreement between added-synthetic and observed profile.

A third alternative (method 'C') is the measurement of the equivalent width of the line of interest, which is then translated into an abundance via the 3D curve-of-growth. For this purpose, we employ the IRAF task splot2 for fitting Gaussian or Voigt profiles to the observed line. Even in the case that the line of interest is contaminated by unidentified components, which cannot be taken into account with methods 'A' and 'B', it may still be possible to use splot with the deblending option to determine the equivalent widths of all components. This procedure may also be useful if the number of blending lines is too large to be treated conveniently by the methods 'A' or 'B' based on synthetic spectra. However, as method 'B', this approach is only justified if all involved lines are weak.

Zirconium abundance from Zr I lines

We looked at the Zr I lines in the sample of Biemont et al. (1981). All these lines are rather weak and blended, or in a crowded spectral region. After inspection, we analysed a subsample of 11 lines, discarding the lines we think to be too heavily blended and/or in a too crowded range. We found our EW measurements to be affected by rather large uncertainties, but for the majority of the 11 lines we agree with the results of Biemont et al. (1981). We think that the cleanest lines are the ones at 480.9, 481.50, 614.0, and 644.5 nm, but also these lines are too poor to give precise information on the solar abundance. We found it difficult to use the flux spectra to measure the EWs of so weak and blended lines, and decided to use only the disc-centre spectra, which are sharper. From the four clean lines mentioned above, we obtain, A(Zr)= 2.65 ± 0.07. If we remove also the 644.5 pm line we obtain A(Zr)= 2.62 ± 0.04. This result is in perfect agreement with the abundance obtained from Zr II lines (see below), but we consider this agreement fortuitous. From the 1D LHD model we find a slightly higher value of A(Zr)=2.67±0.07, and from the 3D model A(Zr)=2.74 ± 0.08.

Investigation of Zr II lines

We picked up all the Zr II lines of Biemont et al. (1981) and the line at 402.4 pm analysed in Ljung et al. (2006). Our Fig. 1 Temporal and horizontal average of the temperature profile of the 3D model (solid line) and the temperature profile of 1D LHD model, shown as a function of (Rosseland) optical depth. In addition, the equivalent width contribution functions (lower, roughly Gaussian-shaped curves) for the 405.0 nm Zr II line are plotted on the same optical depth scale for disc-centre (solid line) and disc-integrated spectra (dashed line).

starting sample was formed by 32 Zr II lines. After inspection of the observed data, and comparison with 1D synthetic spectra, we discarded 11 lines that we think are too heavily blended. After analysis, we discarded four other lines that also Biemont et al. (1981) rejected because they provided a too high abundance. Our final sample is formed by the seven lines of Ljung et al. (2006), that we discuss one by one in the next section, and 10 lines that are used for the abundance analysis in Biemont et al. (1981).

These Zr II lines have lower level energies below 1.7 eV. If we look at the equivalent width contribution functions for disc-centre of the 3D model, we can see that they are formed in a range in τ Ross between -4 and 0, with the maximum in the range -1.0 < τ Ross < -0.5 (see Fig. 1). Only two lines (350.5 and 355.1 nm) are formed over a larger range, with a maximum around τ Ross ≈ -2.0, because these two lines have smaller lower level energies, and are significantly stronger than the others.

The abundance results we obtained from the Zr II lines are summarised in Table 3. We averaged the EW measured in the two disc-centre and the two disc-integrated spectra, respectively. The last column of the table indicates if the result refers to disc-centre or disc-integrated spectra. From the complete sample of 17 lines we obtain A(Zr)=2.612±0.093, with 2.638 ± 0.078 and 2.586 ± 0.100 from disc-centre and disc-integrated spectra, respectively. The line at 347.9 nm gives a very low abundance, while the line at 355.1 nm gives a rather high one. Both these two lines lie in a very crowded region of the spectrum and it is very difficult to place the continuum. When we remove these two lines from the sample, we have 15 Zr II lines and we obtain A(Zr)=2.616±0.060, with 2.643 ± 0.050 and 2.589 ± 0.057 from disc-centre and disc-integrated spectra, respectively.

We note that the 3D abundance from flux is systematically lower than the one from the intensity spectra.

When we look only at the seven lines from Ljung et al. (2006), we find that the line-to-line scatter is reduced with respect to the sample of 15 lines. The 3D abundance we find is 2.653 ± 0.022 and 2.623 ± 0.031 for disc-centre and disc-integrated data, respectively. The abundance averaged over both disc-centre and disc-integrated spectra is 2.638 ± 0.031. We find a very low line-to-line scatter, but again the disc-centre spectra give a systematically higher abundance than the disc-integrated spectra, by an amount comparable to the line-to-line scatter of ≈ 0.03 dex.

We found similar results for Fe, Hf, and Th, where the abundances derived from disc-integrated spectra are lower by 0.02, 0.05, and 0.02 dex, respectively, than the disccentre abundances (see Caffau et al. 2008 andCaffau et al. 2010 for details). A similar behaviour is seen when using 1D model atmospheres. We have presently no obvious explanation for this result.

The zirconium abundance derived from 1D models depends on the choice of the microturbulence parameter ξ. The microturbulence we find for the sample of 15 Zr II lines, by requiring that the relation between EW and A(Zr) obtained with the 1D model has the same slope as that of the 3D model (method 3a in Steffen et al. (2009)), is consistent with the results obtained by Steffen et al. (2009) from a sample of Fe II lines. We adopt ξ = 1.0 km s -1 for the analysis of disc-integrated spectra, and obtain A(Zr) 3D = 2.597 ± 0.053 and A(Zr) LHD = 2.575 ± 0.057. For disccentre data, we choose ξ = 0.7 km s -1 and the result is A(Zr) 3D = 2.652 ± 0.044 and A(Zr) LHD = 2.635 ± 0.045 (for ξ = 1.0 km s -1 we get A(Zr) 3D = 2.602 ± 0.067, A(Zr) LHD = 2.587 ± 0.072).

3D corrections are small in absolute value and depend on the value of the microturbulence. The average 3D correction with respect to the reference 1D LHD model is positive, ∆ LHD = +0.011 ± 0.012 dex, while, with respect to the 3D model it is slightly negative, ∆ gran = -0.008 ± 0.016 dex, again taking ξ = 0.7 km s -1 for disc-centre and ξ = 1.0 km s -1 for disc-integrated spectra.

With the HM model we obtain A(Zr) HM = 2.679 ± 0.044 and A(Zr) HM = 2.619 ± 0.053 for disc-centre and integrated disc, respectively. Taking into account both lines of disc-centre and integrated disc, we obtain A(Zr) HM = 2.649 ± 0.057, and when applying the ∆ gran correction, we have A(Zr) HM = 2.641 ± 0.057.

For our final determination of the solar zirconium abundance, we adopt the complete sample of 15 Zr II lines, and our recommended value is A(Zr) 3D = 2.62±0.06, obtained with 3D model atmosphere.

Analysis of a selected subsample of Zr II lines

We particularly investigated the Zr II lines selected by Ljung et al. (2006) for the solar abundance determination. We find Table 3 The results from the 7 lines of Zr II from Ljung et al. (2006), and other 10 lines from Biemont et al. (1981). this sample interesting because it contains: two mostly clean lines (λ 405.0, 420.8 nm), a line blended on the blue side (λ 449.6 nm), two lines blended on the red side (λ 402.4, 425.8 nm), and two lines blended on both sides (λ 444.2, 511.2 nm)

The Zr II line at 402.4 nm

This weak Zr II line is on the blue wing of a strong blend of Ti I and Fe II where the Ti I line is the dominant component.

On the blue wing of this Ti-Fe blend, there is also a Ce II line. Close by, on the red side of the Ti I line, there is a strong line of iron blended with a much weaker Nd II line.

The zirconium abundance is determined with method 'B' (see Sect. 6): to take into account the effects of the lines present in the range of the Zr II line, we computed a 3D profile for each of the Ce II line at 402.44, the Ti I line at 402.45 nm, and the strong Fe I line at 402.47 nm, to be used in the fitting process. The comparison of the best fitting 3D synthetic profile and the observed Delbouille et al. (1973) disc-centre spectrum is shown in Fig. 2. The Zr abundances and EWs resulting from the best fit of the four observed spectra are listed in Table 3. The situation of the 402.4 nm line may be compared to that of the 425.8 nm Zr II line (see Sect. 9.4), and the 444.2 nm Zr II line (see Sect 9.5), for which we show below that the abundance is underestimated by about 0.03 dex and 0.01 dex, respectively, when using the approximate method 'B'.

The use of splot for the EW measurements is difficult in this case, also with the deblending option. In fact, several lines should be introduced in the range, making the profile fitting unstable.

We also fitted the line profile with a grid of 3D synthetic spectra computed with zirconium only. We limited the fitting to the blue wing and the line core, neglecting (1973) in the region of the Zr II line at 402.4 nm (black symbols), superimposed on the best fitting 3D synthetic profile obtained with method 'B' (solid green/grey).

The graph near the bottom shows the synthetic-observed residuals. The Zr abundance derived from this fit is A(Zr)=2.602.

the red wing, which is blended with the other lines (see Fig. 3). In this kind of fit, neglecting the contribution of the other lines, the continuum placement, when considered as a free parameter, is problematic. When fixing the continuum, the abundance from the fit is in close agreement with the abundance determination from method 'B', in the case of the disc-centre spectra within 0.03 dex (cf. Fig. 3). For the disc-integrated spectra we find the fit to be ambiguous. The abundance comes out larger than from method 'B' by about 0.1 dex, but the fit is convincing only when taking the atlas of Neckel & Labs (1984), although even in this case the core of the line is not well fitted. We find a different result from the two disc-integrated and the two disc-centre spectra, respectively, that we cannot easily explain. Maybe the differences are related to the fact that lines in flux spectra are more broadened than at disc-centre.

The Zr II line at 405.0 nm

This line seems to be unblended. A weaker line is on the red side at about ∆λ 0.016 nm from the centre of the Zr II line, and a stronger Fe I line is also on the red side at about ∆λ 0.035 nm. Looking at the spectra one could conclude that the line is clean and that the abundance can be deduced from line profile fitting with a synthetic grid computed with only Zr and/or by easily measuring the EW. When we fit the line profile with a grid that takes into account only zirconium, the result is very satisfactory for all observed spectra, except perhaps for the disc-integrated spectrum of Neckel & Labs (1984) (see 'NF', Fig. 4). The results are in close agreement with the abundance determinations from the EWs measured with splot. We also compared the theoretical line bisector from the 3D synthetic line profile to the observed disccentre spectrum of Delbouille et al. (1973), and we find an excellent agreement (see Fig. 5). The shift in wavelength between observed and synthetic line bisector is mainly due to the fact that the spectrum of Delbouille et al. (1973) is not corrected for gravitational redshift.

The Zr II line at 420.8 nm

This line is strong and sensitive to the adopted damping constants. The red side is more or less clean, while on the blue side there is a weak line, and close by another weak line.

The line profile fitting is done with method 'A', using Zr only, and restricting the fit to the wavelength range 420.891 . . . 420.913 nm in order to avoid the weak blends in the far blue wing of this Zr II line. The continuum level is fixed in the fitting procedure. The fitting of the disc-integrated line profile in the Kurucz (2005a) spectral atlas is shown in Fig. 6. The abundance derived from line profile fitting and from the EW measured with splot are in very good agreement (within 0.007 dex) for the observed spectra of Delbouille et al. (1973) and Kurucz (2005a). For Neckel Fig. 4 The observed solar disc-integrated (top) and disccentre (bottom) spectra in the region of the Zr II line at 405.0 pm (black symbols), superimposed on the best fitting 3D synthetic profile (solid green/grey) obtained with a grid computed with only zirconium. (1973) for the the Zr II line at 405.0 pm (dashed black), together with its line bisector (plotted on a ten-times expanded abscissa), superimposed on the 3D synthetic profile (solid black) with its line bisector. lower with respect to the abundance derived from the EW measurements by 0.04 dex and and 0.01 dex, for the discintegrated and disc-centre spectrum, respectively.

The Zr II line at 425.8 nm

This Zr II line is on the blue wing of a much stronger Fe II line that has close by, on the red side, a very strong Fe I line, blended with a weaker Fe II line (see Fig. 7). We tried different fitting procedures.

First, we fitted the Zr II line profile with a grid of 3D synthetic profiles computed with only Zr, restricting the fit-Fig. 6 The observed solar disc-integrated spectrum of Kurucz (2005a) (black symbols) in the region of the Zr II line at 420.8 pm, superimposed on the best fitting profile (solid green/grey) derived from a grid of 3D synthetic spectra computed with only zirconium. The fit corresponds to A(Zr)=2.58. ting to the blue wing of the Zr line. For the disc-centre spectrum of Delbouille et al. (1973) we obtain A(Zr)=2.65 ± 0.01. The fit looks quite satisfactory (see Fig. 7), but it is not clear how much the strong lines close by affect the Zr II line, and hence by how much A(Zr) is overestimated.

Second, we applied method 'B', taking into account the contribution of the neighbouring iron lines as a linear combination of scaled 3D synthetic profiles, and obtain A(Zr)=2.61 ± 0.01.

Finally, we also fitted the complete range with method 'A', properly taking into account the three iron lines, with A(Zr) and A(Fe) as free fitting parameters. We obtain a Zr abundance of 2.643 ± 0.001 and 2.627 ± 0.000 from disc-centre and disc-integrated spectra, respectively. Even though the iron lines are not very well reproduced (see Fig. 8), the Zr II line is fitted very well and is not affected by the minor deficiencies of the iron lines. The equivalent widths given in Table 3 are derived from the abundances obtained with this fitting method.

The Zr II line at 444.2 nm

This line is surrounded by two stronger features, an iron line on the blue side at 444.2831 nm and a blend of three lines of iron on the red side. Of these three Fe I lines on the red side, the one at 444.3194 nm is much stronger than the other two components.

Clearly, the Zr II line is affected by the two Fe I features on both sides, and cannot be treated in isolation. Since the two iron features strongly influence the shape of the Zr II line, the Zr abundance cannot be derived from this line by fitting with a grid of 3D synthetic profiles computed with zirconium only.

The most reliable Zr abundance is obtained from the full fitting, method 'A': we fitted the observed Fe-Zr-Fe blend with A(Zr) and A(Fe) as free parameters. The Zr abundance Fig. 7 The observed solar disc-centre spectrum of Delbouille et al. (1973) in the region of the Zr II line at 425.8 pm (black symbols) is superimposed to the 3D best fit (solid green/grey), obtained with a grid of 3D synthetic profiles computed with only Zr. Fig. 8 The observed Neckel & Labs (1984) spectrum in the region of the Zr II line at 425.8 pm (black symbols) is superimposed to the best fitting 3D synthetic spectrum (solid green/grey) obtained with method 'A'. Upper panel: discintegrated spectrum, lower panel: disc-centre spectrum. derived in this way is 2.673 ± 0.004 and 2.650 ± 0.005 for the disc-centre and disc-integrated spectra, respectively (cf. Fig. 9). The equivalent widths given in Table 3 are derived from the abundances of this fitting method.

For completeness, we also derived the Zr abundance with methods 'B' and 'C'.

For method 'B', the profiles of the four iron blending features are represented by scaled 3D synthetic Fe line profiles, which are then added to the profile of the Zr II line in order to obtain an approximate profile of the complete blend (ignoring saturation effects). With this approach, the best fits are obtained with Zr abundances of A(Zr)=2.67 ± 0.01 and A(Zr)=2.64 ± 0.01 for disc-centre and disc-integrated case, respectively. The fit of the observed disc-integrated Fig. 9 The observed solar disc-centre spectrum of Neckel & Labs (1984) in the region of the Zr II line at 444.2 pm (black symbols) with the 3D synthetic best-fit-profile obtained with method 'A' superimposed (solid green/grey). Upper panel: disc-integrated spectrum, lower panel: disccentre spectrum. The residual is below each plot.

spectrum is somewhat better that in the case of the disccentre observations.

In method 'C', A(Zr) is found from the equivalent width of the Zr II line, as obtained from splot with the deblending option. Representing the Fe feature on the red side of the Zr II line either by a single or by three Voigt profiles turned out to have no influence on the resulting zirconium abundance. We find A(Zr) derived in this way to be very close to the result of method 'B' in the case of the disc-centre observed spectra. In the disc-integrated case, A(Zr) is about 0.057 dex smaller than what is obtained with method 'B'.

The Zr II line at 44nm

This Zr II line (EW ≈ 3 pm) is contaminated by a much stronger Cr I line (EW ≈ 9 pm). Certainly, saturation effects cannot be ignored in this case.

Based on the full 3D line profile fitting of the complete blend using method 'A', with A(Cr) and A(Zr) as free fitting parameters, we obtain a Zr abundance of A(r)=2.660±0.003 for disc-centre (cf. Fig. 11), and A(Zr)=2.642±0.005 for integrated disc. These Zr abundances correspond to equivalent widths of 3.11 and 3.41 pm respectively (see Table 3).

The fit resulting from method 'B', scaling the Cr profile and interpolating in the grid of 3D synthetic Zr profiles to achieve the best agreement between synthetic and observed profile, is shown in Fig. 10 for the disc-centre case. It is of similar quality as the one obtained from the full fitting procedure with method 'A' (Fig. 11). Method 'B' gives A(Zr)=2.63 and 2.60 for disc-centre and integrated disc, respectively. As expected, the approximate method 'B' tends to underestimate the Zr abundance, in this example by 0.027 and 0.042 dex for disc-centre and integrated disc,

Fig. 10

The observed solar disc-centre spectrum of Neckel & Labs (1984) in the region of the Zr II line at 449.6 nm (black symbols) with the best fitting 3D synthetic profile obtained from method 'B' with A(Zr)=2.63 superimposed (solid green/grey). The graph in the lower part shows the syntheticobserved residuals.

Fig. 11

The observed solar disc-centre spectrum of Neckel & Labs (1984) in the region of the Zr II line at 449.6 nm (black symbols) with the best fitting 3D synthetic profile obtained from method 'A' superimposed (solid green/grey). The graph in the lower part shows the synthetic observed residuals.

respectively. The average equivalent width we derive with this method for the observed disc-centre and disc-integrated spectra are EW=3.01 and 3.25 pm, respectively, 3.2 and 4.7% lower compared to method 'A'.

When we apply method 'C' using splot with the deblending option, we cannot reproduce the Cr line profile with a single Voigt profile, but we can formally have a good reproduction of the profile with two Voigt profiles, but this option cannot take into account properly the wings of the Cr II line.

The Zr II line at 511.2 nm

For this line, the placement of the continuum is difficult, but the uncertainty induces an error in the EW measure-Fig. 12 The observed solar disc-centre spectrum of Delbouille et al. (1973) in the region of the Zr II line at 511.2 nm (black solid) with the identifications of the lines in the range according the line list of Kurucz. ment of less than 4% and hence a difference in the Zr abundance of at most 0.02 dex. The line shape suggests that it is blended on both wings, and in fact three weak blend lines are expected according to the Kurucz line list3 ; two K I lines at 511.2212 and 511.2246 nm, and a Sm II line at 511.2294 nm. In addition, several weak molecular lines are expected in the range (see Fig. 12).

The average EWs derived with method 'B' for disccentre and integrated disc, respectively, are 8.4 and 9.2 pm, corresponding to zirconium abundances of A(Zr)=2.69 and 2.67. When measuring the EW with splot (method 'C'), we find similar results, but the uncertainty in the measurements is larger, giving an indetermination in the abundance determination of about 0.06 dex.

Conclusions

We analysed the zirconium abundance in the solar photosphere, investigating a selected sample of lines of both Zr I and Zr II. The Zr I lines are weak and heavily blended so that only four of them are acceptable for abundance work. However, the present analysis of the zirconium abundance relies primarily on 15 lines of Zr II that we found suitable for this purpose.

We have applied three different fitting strategies to derive the abundance from Zr lines that are blended by lines of other elements. In the case that all components making up the blend are weak, the different methods give consistent abundances. If, however, stronger lines are involved (for an example see Sect. 9.6), the methods that ignore saturation effects may severely underestimate the Zr abundance by up to 0.05 dex, even though the result of the fitting may appear pleasing to the eye, and the reduced χ 2 may be close to one.

We find a good agreement between A(Zr) derived from the Zr I and the Zr II lines, but, due to the scarcity of Zr I lines, we consider this result as fortuitous. The abundance from disc-centre spectra is systematically higher than the one from disc-integrated observations. A similar result is also found for other heavy elements like Fe, Th, Hf, both with 3D and 1D model atmospheres. This may indicate a problem with the thermal structure of the models, rather than a physical abundance gradient in the solar atmosphere. Further investigations are necessary to find an explanation for this small discrepancy.

Our recommended solar zirconium abundance is based on the 3D result for 15 Zr II lines, and is A(Zr)= 2.62±0.06, where the uncertainty is the line-to-line scatter of the selected sample of Zr II lines. This value is at the upper end of the solar zirconium abundances found in the literature. Still, within the mutual error bars, this result is in good agreement with the meteoritic zirconium abundance of A(Zr)=2.57 ± 0.04 [START_REF] Lodders | Abundances of the elements in the solar system[END_REF]).

Fig. 2

 2 Fig.2The observed solar disc-centre spectrum ofDelbouille et al. (1973) in the region of the Zr II line at 402.4 nm (black symbols), superimposed on the best fitting 3D synthetic profile obtained with method 'B' (solid green/grey). The graph near the bottom shows the synthetic-observed residuals. The Zr abundance derived from this fit is A(Zr)=2.602.

Fig. 3

 3 Fig.3In the top panel, the observed solar disc-centre spectrum ofDelbouille et al. (1973) in the region of the Zr II line at 402.4 pm (black symbols), superimposed on the best fitting 3D synthetic profile (solid green/grey) obtained with a grid computed with only zirconium. In this fit A(Zr)=2.644. From the fit of theNeckel & Labs (1984) disccentre profile, we get A(Zr)=2.669. The lower panel shows a zoom in the fitting range.

Fig. 5

 5 Fig. 5The observed solar spectrum ofDelbouille et al. (1973) for the the Zr II line at 405.0 pm (dashed black), together with its line bisector (plotted on a ten-times expanded abscissa), superimposed on the 3D synthetic profile (solid black) with its line bisector.

Table 1

 1 Lines of Zr I chosen byBiemont et al. (1981) andVanture & Wallerstein (2002)

	λ	E low	EW		log gf
	nm	eV	pm	B	VW	Bog
	350.9331 0.07	0.65 -0.11		-0.21
	360.1198 0.15	1.3 -0.47
	389.1383 0.15	1.7 -0.10
	402.893	0.52	0.06 -0.72
	403.0049 0.60	0.26 -0.36		-0.59
	404.3609 0.52	0.58 -0.37
	407.2696 0.69	0.57 +0.31		+0.24
	424.1706 0.65	0.37 +0.14		+0.07
	450.7100 0.54	0.36 -0.43		-0.46
	454.2234 0.63	0.46 -0.31
	468.7805 0.73	1.00 +0.55		+0.30
	471.0077 0.69	1.05 +0.37		+0.19
	473.2323 0.63	0.25 -0.49		-0.56
	473.9454 0.65	0.55 +0.23		+0.07
	477.2310 0.62	0.53 +0.04		-0.07
	478.494	0.69	0.16 -0.49
	480.589	0.69	0.15 -0.42		-0.63
	480.9477 1.58	0.16 +0.16
	481.5056 0.65	0.20 -0.53		-0.22
	481.5637 0.60	0.30 -0.03
	482.806	0.62	0.19 -0.64
	504.655	1.53 0.050 +0.06		-0.25
	538.5128 0.52	0.18 -0.71
	612.746	0.15	0.21 -1.06 -1.06 -0.87
	613.457	0.00	0.19 -1.28 -1.28 -1.05
	614.046	0.52 0.073 -1.41 -1.41 -0.85
	614.3183 0.07	0.21 -1.10 -1.10 -0.98
	631.303	1.58	0.11 +0.27		+0.18
	644.572	1.00 0.094 -0.83 -0.83
	699.084	0.62 0.050 -1.22		-1.44
	709.776	0.69	0.21 -0.57
	710.289	0.65 0.065 -0.84		-1.06
	743.989	0.54			-1.18
	755.149	1.58			-1.36
	755.473	0.51			-2.28
	755.841	1.54			-1.47
	756.213	0.62			-2.71
	781.935	1.82 0.065 -0.38 -0.39
	782.292	1.75			-1.14
	784.938	0.69	0.10 -1.30 -1.30
	856.859	0.73			-2.80
	857.1085 1.53			-2.07
	858.421	1.86			-1.32
	858.787	1.48			-2.12
	874.958	0.60			-2.79
	350.1133 0.07	0.7 -0.93
	357.5765 0.07	3.3 -0.03
	366.3698 0.15	2.5 +0.01
	588.5629 0.07 0.050 -2.12		-1.82
	B: Biemont et al. (1981)
	VW: Vanture & Wallerstein (2002)
	Bog: Bogdanovich et al. (1996)

Table 2

 2 Lines of Zr II chosen byBiemont et al. (1981) andLjung et al. (2006)

	λ	E low	EW [pm]		log gf
	nm	eV	B	L	B	Bog	L
	343.2415	0.93	2.1		-0.75 -0.51 -0.72
	345.4572	0.93 1.00		-1.34		-1.33
	345.8940	0.96	1.6		-0.52		-0.48
	347.9017	0.53	2.8		-0.69 -1.12 -0.67
	347.9393	0.71	5.1		+0.17 +0.12 +0.18
	349.9571	0.41	2.4		-0.81 -1.08 -1.06
	350.5666	0.16	5.1		-0.36 -0.62 -0.39
	354.9508	1.24	1.6		-0.40 -0.68 -0.72
	355.1951	0.09	5.8		-0.31		-0.36
	358.8325	0.41	2.7		-1.13 -1.25 -1.13
	360.7369	1.24	1.3		-0.64 -0.48 -0.70
	367.1264	0.71	3.2		-0.60 -0.56 -0.58
	371.4777	0.53	3.0		-0.93		-0.96
	379.6496	1.01	1.5		-0.83 -1.17 -0.89
	383.6769	0.56	4.7		-0.06 -0.22 -0.12
	403.4091	0.80 0.65		-1.55		-1.51
	405.0320	0.71 2.38 2.20 -1.00 -0.60 -1.06
	408.5719	0.93 0.54		-1.61 -1.54 -1.84
	420.8980	0.71	4.3 4.26 -0.46		-0.51
	425.8041	0.56	2.6 2.34 -1.13		-1.20
	431.7321	0.71 1.20		-1.38		-1.45
	444.2992	1.49	2.0 2.04 -0.33		-0.42
	449.6962	0.71	3.6 3.15 -0.81 -0.87 -0.89
	511.2270	1.66 0.83 0.78 -0.59 -0.58 -0.85
	363.0027	0.36	3.8		-1.11		-1.11
	407.1093	1.00	2.1		-1.60		-1.66
	414.9202	0.80	7.5		-0.03		-0.04
	416.1208	0.71	5.8		-0.72		-0.59
	426.4925	1.66 1.50		-1.41		-1.63
	444.5849	1.66 0.83		-1.35
	461.3921	0.97 2.91		-1.52		-1.54
	402.4435 0.999	1.20			-1.13
	B: Biemont et al. (1981)
	Bog: Bogdanovich et al. (1996)
	L: Ljung et al. (2006)

  considered only Zr II lines, they used log gf from experimental branching ratios measured by them, and radiative lifetimes taken from different papers. They used 1D Holweger-Müller and MARCS models, and a 3D model computed with the Stein-Nordlund code(Stein & Nordlund 1998) for line profile fitting. EWs of the best 3D fitting profiles are used to estimate the abundance. From 7 Zr II lines, six of which in common withBiemont et al. (1981), they obtain A(Zr)=2.58 ± 0.02 with the 3D model, A(Zr)=2.56 ± 0.02 with MARCS model, and A(Zr)=2.63 ± 0.02 with the HM model. Their adopted value is A(Zr)=2.58 ± 0.02. 5.Velichko et al. (2010) performed the first NLTE analysis of Zr in the solar photosphere. They investigated all the Zr lines analysed byBiemont et al. (1981) andLjung et al. (2006) and selected a subsample of two Zr I lines (424.1 and 468.7 nm) and ten Zr II lines

c 2010 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

www.an-journal.orgc 2010 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

c 2010 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.an-journal.org

http://iraf.noao.edu/ www.an-journal.org c 2010 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

& Labs (1984), the abundance from line profile fitting is c 2010 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.an-journal.org

http://kurucz.harvard.edu/linelists.html www.an-journal.org c 2010 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

Acknowledgements. We acknowledge financial support from "Programme National de Physique Stellaire" (PNPS) and "Programme Nationale de Cosmologie et Galaxies" (PNCG) of CNRS/INSU, France. We thank the referee, K. Lodders, for helpful comments on the relation between cosmochemical and astronomical abundance scale.


Amélie Mathieu
email: amelie.mathieu@agroparistech.fr

Véronique Letort
email: veronique.letort@ecp.fr

Renaud Gosset

Jacques Gignoux

Paul-Henry Cournède

Simulation of morphological plasticity of Acacia tortilis in response to herbivore attacks

Understanding tree architectural plasticity is a major challenge. Functional-structural models can be of great interest to quantify the allocation patterns and their changes according to environment. In this paper, we focus on Acacia tortilis trees in savanna to understand their responses to herbivore attacks both in terms of architecture and allocation pattern. Measurements were carried out in Hwange National Park (Zimbabwe) in 2008. Architectural and biomass data were analysed to quantify the effects of herbivore attacks on plant morphogenesis, and trees were gathered in two groups depending on the intensity of herbivore attacks they suffered. Based on these results, a simple model of herbivore attacks has been implemented and parameters have been estimated in the framework of the GreenLab model and we simulate different realistic behaviours according to the intensity of herbivore attacks.

Introduction

Savanna is a very heterogenous environment where trees are subject to diverse constraints such as fire, hydric stress or herbivore attacks. Trees express a high degree of morphological plasticity due to the variability of the constraints they are subject to. Inside a monospecific population, depending on its environmental constraints, an individual may either develop an intensively branched architecture in order to protect itself against herbivore attacks or conversely invest in underground storage organs to ensure its survival against fire [START_REF] Gignoux | Alternative fire resistance strategies in savanna trees[END_REF], [START_REF] Bond | Fire as a global 'herbivore': the ecology and evolution of flammable ecosystems[END_REF]. The observed morphological diversity cannot be explained by genetic differences [START_REF] Brain | Leaf peroxidase types in Acacia karroo hayne (acacieae, leguminosae) : a rangewilde study[END_REF] as the savanna constraints evolve too quickly compared to the time characteristic for natural selection [START_REF] Abbadie | Lamto : Structure, Functioning, and Dynamics of a Savanna Ecosystem[END_REF]. In this paper, we aim at studying the different strategies trees use to survive in such environment and we focus on the Acacia tortilis. Several studies were conducted on the architectural response of acacia to external constraints [START_REF] Archibald | Growing tall vs growing wide: tree architecture and allometry of Acacia karroo in forest, savanna and arid environments[END_REF], [START_REF] Howaida | Environmental effects on morphology of Acacia tortilis group in the red sea hills, north-eastern sudan and south-eastern egypt[END_REF], but only few of them were concerned with the modification of the allocation pattern mainly because of the experimental difficulty to measure carbon allocation to organs. Quantifying the allocation remains a challenge in ecology and functional-structural models of plant growth can be powerful tools to access these hidden variables. In the South Africa savanna, the amount of wood a tree produces depends on its environment [START_REF] Dangerfield | Shoot characteristics of Acacia tortilis (forsk.) in wildlife and rangeland habitats of botswana[END_REF], [START_REF] Gourlay | Wood production in a natural stand of Acacia karroo in zimbabwe[END_REF]. Plant is the result of balanced allocation between vegetative and reproductive parts [START_REF] Hermes | The dilemma of plants: to grow or defend[END_REF] or between aerial and underground parts [START_REF] Stave | Seedling growth of Acacia tortilis and faidherbia albida in response to simulated groundwater tables[END_REF]. For example, young Acacia karroo subject to pruning will compensate by allocating less biomass to branches but keeping up biomass allocation to leaves [START_REF] Scogings | Effect of water, grass and n on responses of Acacia karroo seedlings to early wet season simulated browsing : Aboveground growth and biomass allocation[END_REF]. Trees were measured in Zimbabwe in 2008 and characterized according to the intensity of the herbivore attacks they have suffered. We first present a statistical analysis of data to estimate the effects of attacks on plant growth. we deduce the parameters of the GreenLab model, which allows us to reproduce the growth and development of plants. A specific version of the model has been developed to reproduce herbivore attacks and we present simulations of the effects of herbivore attacks on virtual plants.

Materials and Methods

Measurements

Acacia tortillis were measured in the wildlife reserve of Hwange (18 • 30'-19 • 53'S, 25 • 48'-27 • 27'E, Zimbabwe) in April 2008. The selected tree species, Acacia tortilis, is one of the dominating species of this area. Its complex architecture varies according to tree developmental stage and environment. Leaves are relatively short and pinnate. In order to estimate the impacts of herbivore attacks on plant morphogenesis, we selected seven young trees that are more sensitive to attacks. Measured trees were characterized according to the proportion of live apices on the total number of apices and distributed in two populations (see table 2.1). We assume that this variable is correlated to the intensity of the herbivore attacks that depends on the location of the tree (distance to village, protection by fences). Both architecture and biomasses were measured on the sampled trees. For the aerial part, detailed measurements of axes of order 1 and 2 were carried out: lengths and masses of internodes, branch position. For the axes of order 3, the number of internodes and the number of lateral axes were counted. The blade surfaces and dry masses were measured for some axes. Due to its depth, the root system was only partly dug up (up to 40 cm), and its total mass was deduced from the weighted mass. Besides, five branches, denoted by E1, E2, E5, E10, E11, were sampled on other trees protected from herbivore attacks to understand variability within the population L. For each blade and internode, dry mass and dimensions were recorded.

The GreenLab model

We recall in this section the main principles of the Green-Lab model. A more detailed description can be found in [START_REF] Cournède | Structural Factorization of Plants to Compute their Functional and Architectural Growth[END_REF]. This model has already been calibrated for real trees with complex architecture, c.f. [START_REF] Letort | Parametric identification of a functional-structural tree growth model and application to beech trees Fagus sylvatica[END_REF] for an example. The tree growth is represented with a discrete time step that is the time needed to build a new phytomer. A phytomer is an internode with a main leaf, spines or hooks at its basis and small axillary leaves. A growth unit (GU ) is a cohort of organs appearing together on the same axis. At each time step, the number of new organs is given by mathematical rules [START_REF] Cournède | Structural Factorization of Plants to Compute their Functional and Architectural Growth[END_REF] and plant architecture is built with geometrical rules. For the acacia trees, we gather the different phytomers into categories, called the physiological ages, according to their morphological properties: the phytomers of the trunk, respectively the large branches and the twigs, are of physiological age 1, respectively 2 and 3. At time step n, the biomass Q(n) produced by the tree is computed with an empirical equation based on the Beer-Lambert law with extinction coefficient k:

Q(n) = E(n)µS p (1 -e -k S(n) Sp ) (1)
S(n) is the total plant leaf area at time step n, E(n) is proportional to the photosynthetically active radiation received by the plant during the time step, S p is inversely proportional to the stand density and µ is an efficiency. This amount of biomass is allocated to plant organs with a proportional model. Each organ receives an amount of biomass proportional to its sink strength multiplied by the total available biomass Q(n) divided by the plant demand D(n), that is the sum of all the sink strengths. A proportion p R of biomass is allocated to the root system.

The ratio of biomass to demand (Q(n)/D(n)) is a convenient variable to describe the feedbacks between functional mechanisms and plant development [START_REF] Mathieu | A dynamic model of plant growth with interactions between development and functional mechanisms to study plant structural plasticity related to trophic competition[END_REF]. In the model, it controls the number of new phytomers and branches. For the acacia tree, we use a simple model allowing delayed starting of buds, i.e. buds may remain dormant during τ time steps and start only if the ratio of biomass to demand is high enough (see equation 2). We assume that young buds are favored and old buds eventually die after τ time steps of dormancy. b p (t, n) denotes the number of buds of physiological age p that appeared at time step n but starts at time step t, α p , β p , γ p are parameters of the model, with γ p < 1 and x is the integer part of x. We write the equation for

t -n + 1 < τ : b p (t, n) = (α p + β p Q(t -1) D(t -1) ) * γ t-n+1 p (2)

Modelling attacks of herbivores

Herbivore attacks were modelled based on field observations and existing models [START_REF] Uriarte | A model of simultaneous evolution of competitive ability and herbivore resistance in a perennial plant[END_REF], [START_REF] Fornara | Ungulate browsing as an ecosystem process: browser-plant-soil interactions in a southern African savanna[END_REF]. Animals will primarily address the young and accessible leaves and stems. Since the studied plants are young and rather small, we suppose that the probability for a phytomer to be eaten does not depend on its height in plant architecture. This parameter could not be neglected for bigger trees. In the Green-Lab model, for each physiological age, herbivore attacks are supposed to depend on three input parameters:

• Number of terminal GU eaten at each attack. The apex of an eaten axis is assumed to be dead.

• Age at which the GU can not be eaten anymore, because herbivores attack primarily young GU .

• Time steps of attack. As we found few data on the frequency of attacks, we set it as an input parameter.

After a strong attack, the tree is assumed to be able to remoblize biomass from the root system. This behaviour is typically observed for acacia tree subject to high environmental constraints. In the model, if the ratio of available biomass to demand (Q/D) is under a given threshold parameter, a constant proportion of root biomass will be remobilized and allocated to the aerial part of the tree.

Results

To quantify changes in plant morphogenesis due to herbivore attacks, some architectural variables were analysed.

Internodes

An analysis of variance on lengths of trunk internodes showed a high effect of the individual on the internode length (P value < 0.001). It could not be linked to the intensity of the attack the tree suffered. Conversely, for young terminal branches, analysis of variance showed an influence of the herbivore attack on internode length (P value = 0.003). Internodes of trees subject to low intensity attacks are significantly higher than similar internodes in population H. Individuals not stressed by attacks seems to invest more in primary growth. This result is coherent with recent studies: trees recently protected from attacks were shown to produce greater shoots [START_REF] Belsky | Overcompensation by plants : herbivore optimisation or red herring ?[END_REF] and removing leaves on young acacia trees induces shorter GU [START_REF] Cooper | The architecture of shrubs after defoliation and the subsequent feeding behavior of browsers[END_REF]. Furthermore, data analysis show that the secondary growth is coherent with the pipe model theory but this process seems to be unaffected by the attacks (no significative differences between the populations).

Leaf masses and surfaces

At the individual scale, the linear relationship between leaf mass and leaf surface is a very good approximation (R 2 > 0.9 for our data set). There is no significant difference in specific leaf mass (SLM) between main leaves and axillary ones. But the SLM is significantly lower for trees of population L compared to population H (Fig. 3.2). An analysis of variance shows that leaf surfaces are smaller for individuals of population H with a mean value of 1.2 cm 2 to be compared to the mean value for population L of 2.4 cm 2 . This may be because the larger leaves were recently eaten on attacked trees, and only the smaller ones remained on the trees.

Parameter estimation

We use our data to estimate parameters of source and sink strength functions to model the growth of acacia trees, although we are aware of the great variability of these values inside the sampled population. Observations show that leaves stay on trees for about 20 Closed symbols:population H, open symbols:population L time steps. For the sake of simplicity, we assume that their sink strength is constant during whole growth (10 time steps). The leaf sink strength is set to 1, and the internode sink strength is 0.8, value given by calibration on real data with methods presented in [START_REF] Guo | Parameter optimization and field validation of the functional-structural model greenlab for maize[END_REF]. As we do not have data on root growth dynamics, we suppose that the root compartment receives 30% of the biomass produced at each cycle, that is the mean experimental value. Regarding the source strength, the model calibration gives µ = 30 and we assume that trees do not suffer from stand competition with other trees, hence S p = 10 m 2 . Difference in internode size between the two populations is incorporated by defining geometrical rules: for the same amount of biomass, an internode of population L will be longer but narrower compared to one of population H.

Simulation

From parameter estimation, we simulate typical plant growths for trees protected from herbivores (Fig 3 .4a) or not (Fig 3 .4b). The attack stops the growth of the terminal axis, and the one is relayed by dormant lateral buds. Attacked tree grows rather in width than in height, and has less phytomers the protected ones because, in our example, dormant buds do not suffice to fully compensate losses.

Discussion

Recent progress in structural functional plant growth models offer new possibilities for exploration in ecology [9]. This paper presents the coupling of an individual plant model with a simple model of herbivore attacks. Experimental data allow us to estimate GreenLab model parameters and reproduce tree growth, with architectural changes according to the intensity of herbivore attacks. However, due to the experimental design, the quantification of the attack intensity is not very precise. Thus measurements should be completed with data on plant growth dynamics. Likewise, definitions of the level of attack intensity should be refined, to allow us to validate this model more thoroughly.

Figure 1 .

 1 Figure 1. Relationships between mass and surfaces. Closed symbols:population H, open symbols:population L

Figure 2 .

 2 Figure 2. Simulation of acacia trees. a: tree protected from herbivore attacks, b: tree attacked three times (8,14, 20 time steps). Both trees are represented at 25 time steps.

Table 1 . Distribution of the individuals in two populations and their proportion of live apices according to their total number of apices. Population L: low intensity of at- tacks. Population H: standard intensity of attacks.

 1

	Plant number S1	S2	S3	S4	S5	S6	S7
	Proportion	-	0.46 0.26 0.15 0.22 0.18 0.50
	of live apices
	Population	L	L	H	H	H	H	L

Acknowledgement

This research is supported by the ANR project 3Worlds.


Sébastien Lemaire
email: lemaire@itbfr.org

Fabienne Maupas

Paul-Henry Cournède

Jean-Michel Allirand

Bertrand Ney

Philippe De Reffye

Analysis of the density effects on the source-sink dynamics in Sugar-Beet growth

Background and Aims This paper aims at studying the capacity of an individual-based morphogenetic crop model of sugar beet growth -GreenLab -to predict the effects of population density on growth and yield.

• Methods A field experiment was carried out for three homogeneous population densities to measure detailed plant development and growth of one cultivar. A set of in situ measurements were used to observe some model inputs (organs' expansion times and life-spans) and destructive data of dry masses were collected for model calibration.

• Keys Results and Conclusions For all density conditions, the plasticity of total biomass production and allocation patterns was accurately simulated using the sets of optimized parameters. Most of them reveal stable across different spacing, and a few parameters vary with density: specific blade mass, time of phyllochron increase, petiole sink and the characteristic surface used to compute Beer's law at individual level. These results are a first step towards developing a predictive capacity regarding crop spacing.

Introduction

The most popular crop models were built from the observation that yield is closely related to the intercepted solar radiation and therefore to the leaf area index (LAI), for example SUCROS [START_REF] Spitters | Simulation and Systems Management in Crop Protection. A simple and universal crop growth simulator: SUCROS[END_REF]. However, this key variable LAI is difficult to estimate and may need some recalibration during growth (with remote sensing techniques for example [START_REF] Launay | Assimilating remote sensing data into a crop model to improve predictive performance for spatial applications[END_REF]) or does not take into account hard stresses (in the context of low input management, in the case of pest attacks or disease). Leaf area controls the interception of radiation, and consequently the early establishment of leaf canopy in the season improves yield. Therefore, any factor reducing the speed of leaf surface expansion directly decreases the final yield. Increase in leaf area index depends on the rate at which new leaves appear and expand, on their final sizes and how long they are retained by plants [START_REF] Milford | An analysis of leaf growth in sugar-beet. I. Leaf appearance and expansion in relation to temperature under controlled conditions[END_REF][4] [START_REF] Milford | An analysis of leaf growth in sugar-beet. III. Leaf expansion in field crops[END_REF]. All these factors strongly depend on the environment (climate, irrigation, fertilization, agricultural practices including crop density…). In a context of sustainable agriculture and low input crop management, it is important to better understand plant physiology under stress conditions. [START_REF] Jaggard | Sugar Beet: Chapter 7 Agronomy[END_REF] reported that if the land area that the foliage of an individual plant covers is represented by a circle, then the maximum radius is 25,5cm. In this condition, all the land is exploited and yield is maximized. Usually, biomass yield increases asymptotically with density, but not sugar yield, which reaches a maximum for 75.000 plants/hectare. The understanding and the extrapolation of these experimental results, in given conditions, for a given cultivar, would benefit from the insights of modeling, and more particularly of individual-based models. This paper aims at studying the capacity of an individual-based morphogenetic model of plant growth, GreenLab, to predict the effects of population density on growth and yield. The GreenLab model can be seen as an intermediate between Process Based Models and Functional Structural Plant Models [START_REF] De Reffye | Plant Growth Models[END_REF]. The model does not claim to be fully mechanistic with regards to physiological and biophysical processes and fluxes involved in plant growth but a particular care is taken to follow empirically the dynamics of the carbohydrate budget, production and allocation [START_REF] Yan | A dynamic, architectural plant model simulating resourcedependent growth[END_REF]. In [START_REF] Lemaire | A morphogenetic crop model for sugar-beet (Beta vulgaris L.)[END_REF], it was shown to be well-adapted to describe sugar beet development and growth. Moreover some preliminary studies used GreenLab to analyze the effects of density on maize growth [START_REF] Ma | Parameter optimization and field validation of the functional-structural model GREENLAB for maize at different population densities[END_REF] and tomato growth [START_REF] Dong | Does the Structure-Function model GREENLAB Deal with Crop Phenotypic Plasticity Induced by Plant Spacing? A case Study on Tomato[END_REF]. GreenLab is thus a good candidate to developing a predictive capacity with respect to population density. The model is briefly recalled in Section 2, especially how the density effect is taken into account in the model [START_REF] Cournède | Computing competition for light in the Greenlab model of plant growth: a contribution to the study of the effects of density on resource acquisition and architectural development[END_REF]. A field experiment was carried out for three homogeneous population densities and is also described in Section 2. The results of these experiments and of model calibration are presented in Section 3. Finally, the predictive capacity of the model is discussed as well as the perspectives of this research.

Materials and methods

GreenLab model for sugar beet growth

The main hypothesis to compute functional growth in GreenLab is that the biomass produced by each non senescent leaf is stored in a common pool of reserves and redistributed among all organs in expansion according to their sink strengths. A sugar beet plant is composed of three types of organs: blades, petioles, and taproot. The initial seed and the green leaves are sources and petioles, blades (including the venation system) and taproot are sinks.

The time unit to compute the ecophysiological functioning (resource acquisition and allocation) is chosen to coincide with the time unit of the morphogenetic sequence based on phytomer appearance. This developmental process is driven by the phyllochron (the thermal time interval between two successive phytomers), and is called growth cycle (GC).

Therefore, the individual plant is described as a discrete dynamical system. At growth cycle n, the empirical equation of neat dry matter production n Q is derived from Monteith's equation [START_REF] Monteith | Climate and the efficiency of crop production in Britain[END_REF] at crop level and is given by: At every growth cycle, the biomass thus produced is allocated to organs individually according to their relative demands called sink strengths. The sink strength of an organ depends on its type "o" ("b" for blade, "p" for petiole and "r" for taproot), the sink variation o P is given classically in GreenLab as a function of its age j (in growth cycles). In sugar beet, the expansion time widely varies from one phytomer to another, and is determined for blades and petioles according to their topological rank k: The sink variation of an organ of type o, of rank k and chronological age j is ( [START_REF] Lemaire | A morphogenetic crop model for sugar-beet (Beta vulgaris L.)[END_REF]):

            - - = p n B p n n S S k S PAR Q exp 1
        = j T T f P j P k o o b a o k o o o , , , ) ( with 1 1 , 5 . 0 1 5 . 0 1 ) ( - -         + -         + = o o o o b o a o o b a T j T j M j f , o M chosen such that: 1 ) ( sup , = j f o o b a j and k o k o T T , max = .
A specific change in biomass allocation at canopy closing leads us to consider a variable petiole sink ( [START_REF] Lemaire | A morphogenetic crop model for sugar-beet (Beta vulgaris L.)[END_REF]):

        + = j T T f I q p j P k o o b a k p p k p o o , , , ) ( ) (
where k I denotes the competition index at growth cycle k defined by:

(1)

(2) (3) (4) k B p k k B p k k S k S S k S I                 - - - = exp 1 1
At each growth cycle n of its expansion period, an organ of age i receives a biomass increment

) , ( i n q o ∆ : n n o o D Q i P i n q ) ( ) , ( = ∆
and its accumulated biomass ) , ( i n q o is thus given by the sum of all these increments since its appearance :

( )

∑ = + - ∆ = i j o o j j i n q i n q 0 , ) , (
If we assume a constant specific blade mass (SBM), the surface area of a given leaf is directly deduced by dividing the accumulated biomass of its blade by SBM. The total green leaf area n S used in the production equation is the sum of the surface areas of all the non senescent leaves. It implies determining for all phytomers the leaf life-span, that is to say the number of growth cycles between appearance and senescence.

Field experiments

Field experiments were conducted in 2008 in the Beauce plain near Pithiviers, France N48°10'12', E2°15'7. A single commercial cultivar, Radar was sown on April 11 with 50 cm between rows and two different lengths between seed-plots,18cm and 12cm, which corresponds to sown densities of respectively 11.1 (classical density) and 16 seeds/m² (high density). The most uniform sections within a large sugar beet field were selected for the trials after plant emergence. This emergence stage (corresponding top the date when 80% of the final population is reached) occurred on March 28. Two weeks after emergence, one section of classical sown density (11.1 seeds/m²) is selected and one plant out of two is removed to get a low density. Finally the final populations were measured to be 5.4, 10.9 and 16.4 plants/m². Irrigation and fertilization were similar for all densities. Daily mean values of air temperature (°C), solar radiation (MJ/m²), relative humidity (%) as well as total daily rainfall (mm) and potential evapotranspiration or PET (mm) were obtained from French meteorological advisory services (Météo France) 5 km away from the experimental site.

Development and growth measurements

Leaf development (appearance, expansion and senescence) was measured weekly non-destructively on the same group of 15 representative adjacent plants and for the three population densities studied. Coloured plastic rings were put around the petioles of the 1 st , 5 th , 15 th , 20 th , and 25 th phytomers for an easier differentiation. Blade lengths and widths as well as petiole lengths were measured on these marked phytomers.

The phyllochron is defined as the thermal time interval that separates the emergence of successive leaves. We consider that a leaf has appeared when its length reaches 10mm. The curve giving the product of blade length time blade width as a function of the thermal time was fitted with a logistic equation (see. Lemaire et al. [START_REF] Lemaire | A morphogenetic crop model for sugar-beet (Beta vulgaris L.)[END_REF]) to determine the expansion times, k b T , and k p T , for a phytomer k. A leaf was supposed senescent when its full surface is yellow. We thus determined k s T , . The expansion time for taproot is obtained from a long period experiment (without harvesting at the usual time) conducted at ITB in 2006 and estimated to be 4527°Cdays.

Destructive biomass measurements were carried out at seven different stages during the growth period: May 20 (340°Cdays, since emergence with base temperature of 0°C), June 11 (703.4°Cdays), June 18 (804.3°Cdays), June 26 (961.8°Cdays), July 9 (1199.4°Cdays), August 12 (1863.3°Cdays) and September 17 (2466.5°Cdays). For the first three dates and the last two, 15 individual plants were selected (randomly) and the dry mass of every individual organ (blades, petioles and root storage) was measured. For the all stages, we measured the fresh matter at the level of organ compartments (total mass of blades that of petioles and taproot mass) for 30 plants; a subsampling of 3 plants was used to determine the rate of dry matter if necessary. Dry matter was obtained by drying for 48 hours at 75°C. The final stage of measurements corresponds to harvest. For this date, the (5) (7) [START_REF] Yan | A dynamic, architectural plant model simulating resourcedependent growth[END_REF] fifteen plants selected for in-situ measures were fully described (at organ level). For all plants measured at organ level, every leaf was digitalized in order to estimate its blade surface area.

Results

3.1.Model parameterization according to population density

Leaf appearance

For the three densities, the number of phytomers is plotted against thermal time in Figure 1. Two distinct intervals can be observed. The phyllochron of the first phase is constant among seasons and experimental treatments (sowing dates, N-content, plant density) as confirmed by three years of ITB experiments in various environmental conditions (data not shown).

In our experiments, the rate of leaf appearance decreases gradually with the increase of competition for light. Many hypotheses were suggested by Milford et al. in 1985 [4] to explain the curve bending: base temperature that changes when the plant gets older, photoperiodic factor and trophic competition (competition for assimilates between the developing storage root and vegetative organs). This competition may slow down the rate of leaf appearance. The density effect is quite strong on this second phase, especially regarding when it starts. Therefore, for the three densities, the growth cycle is based on the common phyllochron=30.8°C.j (see [START_REF] Lemaire | A morphogenetic crop model for sugar-beet (Beta vulgaris L.)[END_REF] for the accurate definition of the growth cycle) and the second phase is modeled by a discrete approximation of the ratio of appearance rates between the two phases -called rhythm ratio -(see also [START_REF] Lemaire | A morphogenetic crop model for sugar-beet (Beta vulgaris L.)[END_REF] for a more accurate description). The start and ratio of these changes of appearance rates vary with density see Table 1.

Density

Leaf expansions and life spans

The important works of Milford in 1985 [START_REF] Milford | An analysis of leaf growth in sugar-beet. I. Leaf appearance and expansion in relation to temperature under controlled conditions[END_REF][START_REF] Milford | An analysis of leaf growth in sugar-beet. II. Leaf appearance in field crops[END_REF][START_REF] Milford | An analysis of leaf growth in sugar-beet. III. Leaf expansion in field crops[END_REF], suggested that density has no effect on leaf expansion time and life spans. Our experiments confirm his observations. The expansion kinetics of the eleventh phytomer (corresponding to the product of its blade length and width against thermal time since emergence) is shown in Figure 2. It is interesting to note that the relative growth rate of expansion is more important in low density than high density, but that blade final sizes are reached at the same time. Such observation is done for all the phytomers.

In Figure 3 are shown the observed leaf life spans of all phytomers for the three densities. The variations between treatments can be neglected.

Specific Blade Mass

The specific blade mass is defined as the ratio of blade dry mass to blade surface area. Figure 4, plots one variable against the other for all measured data (different plants at different ages, one color for each density treatment). We observe a linear relationship between the two variables for the three densities, even though the dispersion increases for bigger leaves. Furthermore SBM decreases with population density. SBM (g/cm²) 0.0082 0.0076 0.0067 Table 2: Estimated SBM for the three studied densities

Beer-Lambert Extinction Coefficient

The Beer-Lambert Extinction coefficient B k depends on the geometric and spectral properties of the canopy. For Sugar Beet, Andrieu et al. [START_REF] Andrieu | Ground cover and leaf area index of maize and sugar beet crops[END_REF] found a relationship linking this coefficient to density, with small variations around 0.7, and deduced that it only has a small incidence on final production. Hodáñová [START_REF] Hodáñová | Structure and development of sugar beet canopy. I. Leaf angle relations[END_REF] reported its value to 0.69 in 1969. Therefore, we choose to fix

Estimation of Hidden Parameters from Experimental Data

In order to fully parameterize the model, we also need to determine the hidden parameters (summarized in Table 3).

Parameter Description µ(g/MJ)

Radiation use efficiency The data on which we calibrate the model correspond to average data on all the individual plants measured at the seven stages: average profile of organ masses when available, and average organ compartments. The calibration is done with Digiplant software [START_REF] Cournède | Structural factorization of plants to compute their functional and architectural growth[END_REF] by model inversion. However, since no specific identification tool dedicated to the study of density exists, the identification process is not fully automatic: the calibration is first done for each density independently, and then we see what parameters are not significantly different for all treatments and can be fixed to a common value. Once this set of common parameters is fixed, we run the calibration process again for each treatment. This can be done iteratively, until fixing the maximum number of parameters while keeping a comparable quality of adjustment to the data.

The parameters thus obtained for the three studied treatments are given in Table 4, and the fitting results are given in Figure 5 for the compartment masses (at 7 stages), and in Figure 6 for the blade organ masses (at every phytomer rank, and at 5 stages). The interesting point when considering the final results of parameter estimation is the strong stability of parameters across density treatment. Only three of the estimated parameters vary.

The first one is the Radiation Use Efficiency. However, we recall that the Radiation Use Efficiency thus estimated is at plant level, with a PAR given in MJ/ m². To get the Radiation Use Efficiency per m², for the closed canopy, we need to multiply µ by p S 0 and by the density, which finally gives very similar radiation use efficiencies of 3.29, 3.18 and 3.16 g/MJ for respective densities of 5.4, 10.9 and 16.4 plantes/m².

The second parameter that varies is p S 0 which is the saturation limit of p n S . Theoretically, we expected that p S 0 would be the inverse of the density, since it should correspond to the space potentially available for the plant. It is not the case, but the ratio of p S 0 to the inverse of the density seems to be (relatively) constant. Finally, the coefficient of petiole sink strength correction also varies, and increases drastically with the density. Such result is in keeping with was observed in maize [START_REF] Ma | Parameter optimization and field validation of the functional-structural model GREENLAB for maize at different population densities[END_REF], and tomato [START_REF] Dong | Does the Structure-Function model GREENLAB Deal with Crop Phenotypic Plasticity Induced by Plant Spacing? A case Study on Tomato[END_REF]: in order to overcome competition for light, a plant will tend to invest more in its petioles, to intercept more light.

4.Conclusion

This study is a first step to develop a predictive model that would allow the inference of source-sink dynamics of Sugar Beet growth at any density. Of course, there is a strong need for validation among seasons and in different environmental conditions. However, we can exhibit the base of what could be such model. It would consist in determining stable relationships linking SBM, the decrease of the rate of leaf appearance, petiole sink strength correction and p S 0 to the density. The GreenLab model could thus be applied straightforwardly with density as an input. It thus opens perspectives for applications, to optimize density, especially in the context of low input crop management, since a low density crop should also require less water, nitrogen, pesticides … However, the results obtained seem to show some hints of improvements, in terms of parameter stability and predictive capacity. Regarding the decrease of phytomer appearance rate, it does not seems directly related to density but more to some competition index, since the differences of the rate appearance (during the second phase of development) were also observed in other experiments to study Nitrogen stresses: plants enduring Nitrogen stresses had smaller foliage and more phytomers than plants without stresses.

Likewise, the coefficient of petiole sink strength correction was introduced in Lemaire et al. [START_REF] Lemaire | A morphogenetic crop model for sugar-beet (Beta vulgaris L.)[END_REF] to take into account the variation of petiole sink strength after canopy closing, by introducing this dependence to a competition index. We would have expected that this competition index, which includes density effect, would be sufficient to fully take into account the density effect. Estimation results show that it is not the case. Finally, as underlined above, we expected to find a better agreement between p S 0 and the inverse of density. Even though the ratio between the two variables seems constant for different densities, the results may indicate that we are missing an element in the individual-to-population extrapolation.

To conclude, we would like to insist on the strong stability of GreenLab model parameters in the different density conditions, few of them are sufficient to take into account the effects of competition for light.

Parameter stability is a key issue in plant growth modeling and its link to plant genetics [START_REF] Letort | Quantitative genetics and functional structural plant growth models: simulation of quantitative trait loci detection for model parameters and application to potential yield optimization[END_REF], since each genotype should be represented by a set of response parameters that are valid under a wide range of conditions [START_REF] Tardieu | Virtual plants: modelling as a tool for the genomics of tolerance to water deficit[END_REF].

  photosynthetically active radiation accumulated during growth cycle n, µ is plant Light Use Efficiency (LUE), n S is the leaf surface area of the plant at cycle n, B k is the Beer-Lambert extinction coefficient and p S is an empirical coefficient corresponding to a characteristic surface (related to the two-dimensional projection of space potentially occupied by the plant onto x-y plane). Since this surface increases with the development of foliage, and saturates at canopy closing, we propose to model α are empirical parameters estimated from experimental data.

T

  , will denote the life-span (in growth cycles) of the leaf borne by the phytomer of rank k. If the phyllochron, expansion duration ( and specific blade mass can be observed experimentally, it is not the case for the parameters: µ , will be estimated from experimental data by model inversion, as detailed by Guo et al.[START_REF] Guo | Parameter optimization and field validation for the functional-structural model GREENLAB for maize[END_REF].

Figure 1 :

 1 Figure 1: Number of phytomers as a function of the thermal time for the three density treatments

  u la te d b la d e s ize (L e n g th x w id th in m m ²) 5.4 plants/m² 10.9 plants/m² 16.4 plants/m²

Figure 2 :

 2 Figure 2: Expansion kinetics of the blade of the 11th phytomer for the 3 studied densities

Figure 3 :

 3 Figure 3: Leaf life spans (in growth cycles) according to phytomer rank for the 3 studied densities

Figure 4 :

 4 Figure 4: Blade Masses against Surface Area for all the measured leaves and for the three studied densities

Figure 5 :

 5 Figure 5: Average plants at compartment level: experimental data (dots) and simulated data (lines) for three density treatment. (a) Taproot compartment dry mass and (b) Shoot (i.e. blades and petioles dry mass) as function of thermal time since emergence

Figure 6 :

 6 Figure 6: Average plants of shoot dry masses organ level according to phytomer rank: experimental data (dots) and simulated data (lines) for three density treatments (a) 5.4 plants.m -2 (b) 10.9 plants.m -2 and (a) 16.4 plants.m -2

Table 1 :

 1 Estimated rhythm ratio and corresponding phytomer rank to the start of change of appearance rate for the three studied densities

		Rhythm	Start of 2 nd phase
	(plants/m²)	Ratio	on phytomer rank
	5.4	0.57	18 th
	10.6	0.57	21 st
	16.4	0.58	35 th

Table 3 :

 3 Summary of model hidden parameters to estimate