Dr Yannick Geiger

Thierry Achard

Dr Aline Maisse-François

Stéphane Bellemin-Laponnaz

Dr Dr T Achard

Absence of Non-Linear Effects Despite Evidence for Catalyst Aggregation

INTRODUCTION

Non-linear effects (NLEs) appear in asymmetric catalytic reactions when the product ee (eeP) does not scale linearly with the catalyst ee (eeL, as the chirality is often carried by a chiral ligand attached to a metal). [1][2][3][4] In case of a positive non-linear effect (or (+)-NLE) the eeP can be much higher than eeL and this leads to asymmetric amplification (Figure 1). On the other hand, a small reduction of the ligand's enantiopurity can also lead to a reduced eeP (negative NLE, (-)-NLE, chiral depletion). The concepts and first experimental examples of NLEs in asymmetric catalysis were reported by Kagan and co-workers in the mid-80's and early 90's. [1,2] Their work followed studies on the influence of a compound's enantiomeric purity on its physical properties by Horeau and Guetté [5] and on its chemical reactivity in stoichiometric reactions by Wynberg and Feringa in the mid-70's. [6] It should also be noted that already in 1936 (!) Langenbeck and Triem reported a NLE in the double addition of menthol to oxalyl chloride, along with the conceptual basics for asymmetric amplification and an outlook to catalytic reactions. [7,[START_REF] Ref | 7 went largely unnoticed until being highlighted by Heller and co-workers in 2000[END_REF] Nonlinear scaling in enantioselectivity is also widely discussed in relation to the origin of biological homochirality, in which (+)-NLEs as steps of chiral amplification might have played a role. [9][10][11] Figure 1. Scheme for typical product ee (eeP) vs ligand ee (eeL)-relationships. Non-linear effects (here, (-) or (+)-NLE; hyperpositive NLEs are not shown) are indicative of catalyst aggregation, which leads commonly to the assumption that the absence of NLE goes hand in hand with the absence of catalyst aggregation.

Although the concepts of non-linear effects are also valid in stoichiometric reactions and kinetic resolution processes, they have been mostly studied for catalytic reactions. This is mainly due to their usefulness in the analysis of catalytic reaction mechanisms: the presence of a NLE is generally considered to be a footprint of catalyst aggregation or catalysis by multiligand-complexes. [3,4] Establishing a NLE curve requires nothing more than having both enantiomers of the ligand, and the analysis is apparently simple and binary: if there is an NLE there is also catalyst aggregation; if not, there is no aggregation. Therefore, the use of NLEs as a probe has has been successfully applied in the study of asymmetric catalytic mechanisms, especially when other analytical tools failed. [12] Although this approach leads to correct conclusions in most cases, it may also lead to false interpretations as it consists in the simplification of a complex system. This has been recently highlighted by Fu and Kalek who reported a special case of kinetic resolution, where an intrinsic (-)-NLE was observed without any catalyst aggregation. [13] In this paper, we present a specific case complementary to that of Fu and Kalek, [13] in the form of an asymmetric catalytic system in which the catalyst does aggregate despite an observed linear eeP/eeL-relationship. We will first discuss how such a situation may arise according to the Kagan and Noyori models for NLEs. Then we will discuss experimental evidence for such a case, namely in the catalytic asymmetric addition of ZnEt2 to benzaldehyde in presence of the N-methylephedrine (NME) ligand. Since this catalytic reaction follows the model we recently developed, [14][15][16] we will use plot simulations to explore possible scenarios leading to linearity. Overall, the present contribution shows that the study of non-linear effects requires much more caution as it can lead to false interpretations in particular with respect to linearity as an indicator of non-aggregation.

RESULTS AND DISCUSSION

Aggregation without NLE: known theoretical cases. The ML2 model, developed by Kagan and co-workers in 1986, [1,2] assumes the assembly of two chiral ligands within a complex or merely a dimerization of a complex. In such system, three stereoisomers are formed (Figure 2): two homochiral dimers (RR and SS, which give rise to products with ee2 and -ee2, respectively) and one heterochiral dimer (RS, which yields racemic product). The chiral product is thus obtained from these three catalysts that operate simultaneously at a steady state. The system's overall outcome depends on the dimers' relative reactivity g and on the relative catalyst concentration β. The system will always exhibit a NLE unless the RS-dimer does not form (i.e. β = 0) or if its reactivity is equal to that of the homochiral dimers (g = 1). In these cases, the eeP vs eeL-relationship is linear despite dimers being present. In 1989, Ryoji Noyori and co-workers observed unprecedented strong (+)-NLEs in the DAIBcatalysed enantioselective addition of dialkylzinc reagents to aldehydes . [17,18] They proposed a model in which the reaction is catalysed by monomeric complexes while dimers (both homochiral and heterochiral) are inactive and in a dynamic equilibrium with the monomers (Figure 3). In this case, the NLE is governed only by the relative stabilities of the dimeric species (i.e. dimerisation constants, KHomo and KHetero), which in turn alter the R/S ratio of the catalytically active monomers (although the association of reagent, substrate and catalyst can also play a role, cf. XXX). A NLE occurs when the heterochiral dimer is either more (i.e. KHetero > 2KHomo) or less (KHetero < 2KHomo) stable than the homochiral species.. The particular case of KHetero = 2KHomo (what Kagan called "a statistical distribution of ligands" in the ML2 model, [2] cf. Supporting Information I.2.4 ) was highlighted as the "branching point" of the system: the eeP vs eeL-relationship is strictly linear here. Dimerization takes place, but in a way that it doesn't influence the original ratio of the monomeric Rand S-catalysts.

Figure 3. General scheme for the model developed by Noyori. The overall eeP is a consequence of the R/S-ratio of monomeric catalysts -the dimers are inactive -which depends on the ligand ee (eeL), the catalyst's enantioselectivity (ee1), the dimerization constants KHomo and KHetero, the total catalyst concentration [Cattot] and the catalyst-substrate-reactant association constant Kassoc. Cf.

Supplementary Figure 2 for more details and section Supporting Information II.2.5 for a comment on the adapted nomenclature.

A case study: the enantioselective addition of ZnEt2 to benzaldehyde catalysed by (-)-NME.

We are currently studying the addition of dialkylzinc reagents to benzaldehyde catalysed by N-substituted ephedrine ligands. Although these are chiral aminoalcohols similar to Noyori's DAIB, they follow a different mechanistic scheme in which monomeric as well as dimeric complexes catalyse the reaction with different enantioselectivities (or even giving different product enantiomers), thus being in equilibrium and in competition. [14][15][16] Figure 4 displays the proposed reaction network. This model stands in striking contrast to the Noyori model, which excludes dimer catalysis, and thus opens new possibilities in non-linear phenomena, the most notable being a hyperpositive NLE with N-benzylephedrine (NBE) [14,15] and an enantiodivergent NLE with the parent ligand N-methylephedrine (NME) when using ZnMe2 as the dialkylzinc reagent. [16]

Figure 4. General scheme for our monomer-dimer competition model, with equations (1), (2) and (3) that allow to simulate eeP vs eeL-curves. Both monomers and dimers are catalytically active, thus eeP depends on the ligand ee (eeL), the catalysts' enantioselectivities (ee1 and ee2, respectively) and relative kinetic constants (k2/k1, k2'/k1), the total catalyst concentration [Cattot] and the dimerization constants KHomo and KHetero. α and β (with α = [R]+[S] and β = [R][S]

) are working variables introduced to simplify the expressions. Furthermore, α links eeP and eeL through equations ( 1) and ( 2

): a valid range of αvalues must be defined first so that paired eeP-and eeL-values can be calculated then (cf. Supplementary Information II.2.2 for updated expressions of the limiting values α0 and α100).

In the light of these observations, it is worth to reconsider an earlier report from Fitzpatrick and co-workers [19] who have used NME in conjunction with ZnEt2 and found an "essentially linear" relationship between eeP and eeL. Cryoscopic molecular weight (MW) measurements showed that the NME-ZnEt complex has a MW of 770 g/mol, which the authors interpreted to be a 1:1-mixture of dimers and tetramers -trimers were excluded without further explanation -which, however, break down to monomers upon addition of one equivalent of ZnEt2. This led them conclude that NME-ZnEt does not form any aggregates under catalytic conditions, consistent with the absence of an NLE. This is rather surprising: why would switching from ZnMe2 to ZnEt2 result in a total loss of aggregation? With NBE as ligand, the homochiral dimeric complex is present as well as catalytically active with both ZnMe2 and ZnEt2, even though the dimer's contribution to the overall eeP is probably lower with ZnEt2 (as seen in a hyperpositive NLE with a weaker amplitude). [15] This inspired us to look more closely at this case and to conduct our own investigations of the NME/ZnEt2-system. Using standard reaction conditions (0.83 M benzaldehyde in toluene, 1.2 equiv. ZnEt2, 20 mol% ligand) we first could reproduce the linear eeP vs eeL-relationship (Figure 5a, blue dots) previously noted by Fitzpatrick; for comparison, we have added the data for the ZnMe2 addition for similar reaction conditions (Figure 5a, grey dots) which shows an inverted (-)-NLE. [16] When using only enantiopure NME, a reduction of the catalyst loading led to a significant increase of eeP (Figure 5b, blue dots). While the alcohol product was obtained in 68% ee when using 20% mol% of (-)-NME, 80% ee was measured in presence of only 2.5 mol% ligand. Such behaviour is typical of a monomer-dimer-competition system, as shown in Figure 4, in which a dimeric catalyst is present and gives a lower eeP than the monomer: a decrease in catalyst concentration tips the equilibrium from the dimer to the more enantioselective monomer, which leads to an increase of the overall eeP. A similar eeP increase is obtained with NBE as ligand using the same reaction conditions (Figure 5b, grey crosses) [14,15] and the NME/ZnMe2-system also shows this behaviour (Figure 5b, grey dots). [16] This points to the presence, and also to the catalytic activity, of homochiral NME-ZnEt dimers -even though no NLE is observed for this system in Figure 5a. [16] ) or (-)-NBE (ZnEt2, grey crosses). [14,15] Each point is the mean of three independent experiments; the vertical bars depict standard deviations. The linear (a, blue dotted lines) and 2 nd -order polynomial fits (a and b, grey dotted or dashed lines) serve as visual guidelines.

Next, we investigated the behaviour of the complexes in solution using 1 H and DOSY NMR spectroscopy. The isolated enantiopure (-)-NME-ZnEt complex in toluene-d8 (Figure 6a) shows a clear trimodal pattern which is readily seen in the -CH-O-region of the spectra (d = 4.8 to 5.2 ppm). The integrals of the three doublets are identical, even if the concentration is decreased from 185 to 6 mM (cf. Supporting Information). This indicates that the three sets of signals do not arise from different species in equilibrium (i.e. monomer and dimer) but from a single molecule containing several non-equivalent NME-ZnEt-subunits -in this case a trimeric aggregate, for which ample precedent exists with other aminoalcohols and organozincs. [20][21][22][23][24][25] This confirms Fitzpatrick's MW measurements, even though they excluded the possibility of trimer formation. Our own measurements using 1 H DOSY NMR and the ECCmethodology developed by Stalke and co-workers [26][27][28] gave a mean molecular weight (mMW) of 732 g/mol at a concentration of 185 mM, which corresponds to a degree of aggregation (Ad) of 2.7. This is consistent with Fitzpatrick's cryoscopic molecular weight measurements, even though they did not interpret the data as resulting from the trimerization of the complex. At 6 mM, 619 g/mol and Ad = 2.3 were obtained, indicating that the trimer is in equilibrium with smaller aggregates. However, addition of 5 equiv. of ZnEt2 didn't disrupt the trimer: the trimodal pattern persists in the 1 H NMR spectrum at 185 mM (Figure 6b) and the mean MW decreased only slightly (701 g/mol, Ad = 2.6)). This is in contradiction to Fitzpatrick's conclusions. On the other hand, adding benzaldehyde instead of ZnEt2 gave rise to several new signals (Figure 6c) and the DOSY NMR spectrum became more disparate, with sets of signals at ca. 665 (Ad = 2.4) and ca. 753 g/mol (Ad = 2.8). We subsequently performed catalytic runs (i. e. NME-ZnEt in presence of both ZnEt2 and benzaldehyde) monitored by 1 H NMR. Figure 7 shows the spectra (a) before, (b) 5 min after and (c) 6 h after the beginning of the reaction. Upon addition of the benzaldehyde substrate the CH-O region dramatically changes: the trimodal pattern broadens and almost totally disappears. At the end of the catalytic reaction the trimer signal becomes more clearly defined (ca. 20% of its original intensity). This is a strong indication for the trimer being only a resting species which dissociates to the actual catalysts in presence of substrate and reagent, and which partly re-forms once the latter are consumed. At this stage we could not assign the various broad resonances to possible catalytic species, but the complexity of the spectrum renders the hypothesis unlikely that a single monomeric complex is present. A DOSY analysis after completion of the reaction (the reaction is too fast to be analysed during the conversion) resulted in a complex spectrum, leading to average MW values between 679 and 805 g/mol (Ad = 2.5 to 3.0, cf. Figure 7c).

Figure 7. 1 H NMR spectrum of the catalytic addition of ZnEt2 to benzaldehyde catalysed by NME-ZnEt, a) before addition of benzaldehyde, b) after 5 minutes reaction time, c) at the end of the reaction (ca. 6h reaction time), with the MW values obtained for each peak from DOSY NMR. The part of the spectrum depicted here shows the O-CH region of the catalyst as well as of the reaction product (4.6-4.7 ppm). The vertical zoom of all three spectra was normalised with an internal reference (mesitylene).

Modelling the absence of NLEs

Altogether, the variation of the catalyst loading and the NMR data show that NME-ZnEt complex does aggregate under reaction conditions, despite the eeP vs eeL-relationship being linear and in contrast to conclusions drawn in earlier reports. Therefore, the latter must be the result of a "serendipidous" coincidence: the parameters governing the reaction are related in a way that the curve does not deviate from linearity despite catalyst aggregation -very similar to the examples shown for the Kagan and Noyori models vide supra. This drove us to investigate such particular parameter settings using our own model for monomer-dimer competition reactions, which has allowed us to simulate hyperpositive NLE curves, eeP vs [Cattot] curves and kinetic profiles. [16] The following simulations have been calculated using equations ( 1)-(3), cf. Figure 4. It should be noted that this set of equations does not allow to set KHetero = 2KHomo because the denominator in eq. ( 3) would become equal to 0. This issue can be easily circumvented but, for the sake of completeness, we have developed a new set of expressions which includes the case of KHetero = 2KHomo ("Model II.b", cf. section II.2.3 of the Supplementary Information). Also, we abstained from extending the model to trimers (either catalytically active or not) as it would significantly complexify the discussion. Since the observed NME-ZnEt-trimer probably is only a transient species, the present monomer-dimer competition model should be sufficient for the purpose of this study.

The present theoretical study shows that there are various scenarios which allow the emergence of a linear eeP vs eeL-relationship, which depend on the contribution of the different catalysts to the overall reaction outcome.

The contribution depends on the different catalyst concentrations ([R], [S], [RR], [SS] and [RS]

) and their respective kinetic constants (k1, k2 and k2'). The catalyst concentrations are considered to be constant over time and depend on the interplay of the dimerization constants (KHomo, KHetero) and the overall catalyst concentration [Cattot]. The monomer's and homochiral dimer's enantioselectivities (ee1 and ee2) have been set to 100% and -100%, respectively. Globally, we can distinguish between two different extreme cases: one in which only the monomers and other in which only the dimers are catalytically active. In those extreme cases the system behaves either like the Noyori or Kagan's ML2 model, respectively, and the eeP vs eeL-plots become linear for the parameters as set above.

For a negligible activity of both homo-and heterochiral dimers (k2/k1 = 0.01 and k2'/ k1 = 0.01), the plot becomes linear if KHetero/KHomo = 2 -exactly as in the Noyori model -regardless of the absolute value of KHomo, of KHetero or of the kinetic constants. This changes when k2/k1 is increased: the necessary KHetero/KHomo ratio deviates from 2 and becomes dependent on the absolute values of KHomo, KHetero and [Cattot]. Figure 8a shows a variation of k1/k2 in a system with only weak dimerization (KHomo = 10, [Cattot] = 0.01). The KHetero/KHomo ratio necessary to keep the plot linear decreases progressively to 0 as k2/k1 is increased, the plot's slope also decreases. Interestingly, once it has crossed the eeL axis and gets a negative slope, no linear behaviour is possible anymore; the dashed lines in Figure 8a show plots which gets as close to linear as possible (cf. Supporting Information section I.2.1 for the fitting procedure used here). For k2/k1 = 20 one obtains a seemingly linear plot when setting KHetero/KHomo to 14.1 (a close zoom shows that it is slightly S-shaped) but for higher k2/k1-values the plots become distorted. A similar behaviour is found when working in conditions with high levels of catalyst aggregation (i.e. KHomo = 1000 and [Cattot] = 0.16, Figure 8b). For k2/k1 = 0.01 the optimal KHetero/KHomo for a linear plot is again 2. However, here the system is much more sensitive to changes in k2/k1. At k2/k1 = 0.2 the plot's slope is negative and, again, only seemingly linear. At higher values of k2/k1 (0.4 to 10, dotted lines in Figure 8b), the system's behaviour differs from Figure 8a: close to linear plots are obtained by setting KHetero = 0.

The precedent conclusions prompted us to study how variation of k2' (i.e. racemic contribution) may drive system towards linearity. Figure 9a shows the exact same simulated curves as in Figure 8a (dashed lines in both figures) with k2' = 0.01. Their parameter sets were taken to calculate another set of curves (Figure 9a, full lines) with increased k2'-values in order to obtain more linear plots. Increasing k2'/k1 indeed flattens the curves and reduces the difference (i. e. the sum of squares, cf. Supporting Information I.2.1) from linearity, but only to a certain extent. This is due to the high KHetero/KHomo-ratio which twists the curve into an "S" shape (and makes it hyperpositive [16] ) and which cannot be fully compensated by tweaking k2'/k1. However, for KHetero/KHomo < 100 the S-shape is negligible and linear plots are obtained for any k2/k1-ratio if the correct k2' is chosen (Figure 9b), regardless of the plot being in the upper or lower half of the graph. Globally, we can say that especially when the meso dimer is active, a multitude of parameter configurations exists where the eeL vs eeP-plot becomes linear although catalyst dimers are present (and possibly catalytically active) in the reaction mixture.

Conclusion

The present study shows that probing for non-linear effects alone can lead to false interpretations considering the aggregation state of the catalyst in a catalytic reaction: although the absence of a NLE is commonly interpreted as the catalyst not aggregating, we show here an example where this is not the case. The NME-catalysed addition of ZnEt2 to benzaldehyde displays a linear eeP vs eeL correlation although catalyst aggregates are present and, moreover, catalytically active as seen in the eeP vs catalyst loading screening and the 1 H and DOSY NMR studies. Kagan and Noyori already had pinpointed such exceptions to the rule in their respective models when applying a particular set of parameters. Our own model confirms these exceptions and shows that, in the presence of both monomeric and dimeric entities, a plethora of different parameter combinations can lead to linearity. What we can derive from this study is that probing for non-linear effects remains a useful tool to find evidence for catalyst aggregation, but it should be complemented with other analytical methods in order to avoid erroneous conclusions. In this study, 1 H and DOSY NMR proved to be such additional and valuable tools, and eeP vs catalyst loading-plots have now become good indicators for the possible catalytic activity of both monomeric and aggregated catalyststhus, they also provide evidence about the catalyst's degree of aggregation. Kinetic studies have also been used to complement NLE probing: relationships such as rate vs eeL [29] and catalyst order vs catalyst concentration [14,16,30] may show features that are indicative of certain models.

Table of Contents Graphic

One might think that the absence of a non-linear effect (NLE) might necessarily mean that only monomeric complexes catalyse a studied enantioselective reaction. Here we show a counterexample using an ephedrine-derived ligand in asymmetric alkylzinc additions. The system is catalysed by both monomers and dimers, although no NLE is apparent.

Figure 2 .

 2 Figure 2. General scheme for the ML2 model developed by Kagan. The overall eeP will depend on ligand ee (eeL), the enantioselectivity of the dimer catalyst (ee2), the catalysts' concentrations [RR], [SS] and[RS] (related through the constantβ) and their relative reactivity g. Note: the nomenclature of all constants in this paper has been unified for the sake of clarity and thus differs from Kagan's original notation, cf. Supplementary Figure1for more details.

Figure 5 .

 5 Figure 5. a) eeP vs eeL-relationship in the (-)-NME-catalysed enantioselective addition of ZnEt2 (blue dots) or ZnMe2 (grey dots) to benzaldehyde, b) catalyst loading screening of the same reaction using enantiopure (-)-NME (ZnEt2: blue dots, ZnMe2: grey dots[16] ) or (-)-NBE (ZnEt2, grey crosses).[14,15] Each point is the mean of three independent experiments; the vertical bars depict standard deviations. The linear (a, blue dotted lines) and 2 nd -order polynomial fits (a and b, grey dotted or dashed lines) serve as visual guidelines.

Figure 6 .

 6 Figure 6. 1 H NMR spectra of enantiopure (-)-NME-ZnEt (a), with additional 5 equiv. ZnEt2 (b) or 5 equiv. benzaldehyde (c). All analyses were performed in anhydrous toluene-d8; the complex and ZnEt2/Benzaldehyde concentrations (0.17 and 0.83 M, respectively) are the same as used for the catalytic runs in Figure 5. The vertical zoom of all three spectra was normalised on the mesitylene reference peak. The integrated peaks correspond to the O-CH (δ = 4.8 to 5.2 ppm), N-CH (3.5 to 2.5 ppm) and N-CH-Me signals (1.0 to 0.6 ppm) and illustrate well the trimodal pattern of the spectrum.

Figure 8 .

 8 Figure 8. Screening of the necessary KHetero/KHomo-ratio to obtain linear eeP vs eeL-plots at different k2/k1 ratios, a) in conditions of low catalyst aggregation (low KHomo, low [Cattot]) and b) high catalyst aggregation (high KHomo, high [Cattot]). The non-full lines (dashed and coloured in a), dotted in b)) represent cases were linearity cannot be achieved, they show instead the curve closest to linearity which could be obtained. All curves have been obtained by setting the indicated parameters (k1 = 100, k2' = 1, ee1 = 100%, ee2 = -100%) and by fitting the curve, through variation of KHetero, to a simulated linear plot. The full procedure and data are found in the Supporting Information.

Figure 9 .

 9 Figure 9. a) Non-linear curves from Figure 8a (dashed lines ; with KHomo = 10, [Cattot] = 0.01, k1 = 100, k2' = 1, ee1 = 100%, ee2 = -100%) and the corresponding modified curves with adapted k2' (full lines) in order to minimise deviation from linearity, b) k2/k1-screening of the system by fitting to linearity through variation of k2' (k2 = 100, KHetero = 400, KHomo = 100, ee1 = 100%, ee2 = -100%).

Acknowledgments

We acknowledge the CNRS and the Ministère de l'Enseignement Supérieur et de la Recherche (MESR) for a PhD grant to Y.G. This work was funded by the French National Research Agency (ANR) through the Programme d'Investissement d'Avenir under contract ANR-11-LABX-0058_NIE within the Investissement d'Avenir program ANR-10-IDEX-0002-02.

Data availability

Synthetic and mathematical procedures, as well as additional data, chromatograms and NMR spectra ( 1 H and DOSY) can be found in the SI.docx. Detailed experimental data (reactant quantities, reaction conditions, raw and treated results for all catalytic runs, including reactions with NBE) and calculated data for all simulated graphs in this paper can be found in the respective Microsoft Excel files (SI_Experimental_Data.xlsx, SI_Simulation_Data.xlsx).


Ahmed Attia

Yann Nouvellon

Santiago Cuadra

Osvaldo Cabral

Jean-Paul Laclau

Joannès Guillemot

Otavio Campoe

José-Luiz Stape

Marcelo Galdos

Rubens Lamparelli

Guerric Le Maire

Modelling carbon and water balance of Eucalyptus plantations at regional scale: effect of climate, soil and genotypes

Keywords: Eucalyptus plantations, ecophysiological model, G'DAY, optimization, productivity

Carbon and water budgets of forest plantations are spatially and temporally variable and hardly empirically predictable. We applied G'DAY, a process-based ecophysiological model, to simulate carbon and water budgets and stem biomass production of Eucalyptus plantations in São Paulo State, Brazil. Our main objective was to assess the drivers of spatial variability in plantation production at regional scale. We followed a multi-site calibration approach: the model was first parameterized using a detailed experimental dataset. Then a subset of the parameters were re-calibrated on two independent experimental datasets. An additional genotype-specific calibration of a subset of parameters was performed. Model predictions of key carbon-related variables (e.g., gross primary production, leaf area index and stem biomass) and key water-related variables (e.g., plant available water and evapotranspiration) agreed closely with measurements. Application of the model across ca. 27,500 ha of forests planted with different genotypes of Eucalyptus indicated that the model was able to capture 89% of stem biomass variability measured at different ages. Several factors controlling Eucalyptus production variability in time and space were grouped in three categories: soil, climate, and the planted genotype. Modelling analysis showed that calibrating the model for genotypic differences was critical for stem biomass prediction at regional scale, but that taking into account climate and soil variability significantly improved the results. We conclude that application of process-based models at regional scale can be used for accurate predictions of Eucalyptus production, provided that an accurate calibration of the model for key genotype-specific parameters is conducted.

Introduction

Forest vegetation plays a major role in determining the state of the global climate system and carbon cycle, both of which are undergoing significant anthropogenic perturbations.

Among forest vegetation types, Eucalyptus is the most widely planted tropical hardwood genus, covering more than 20 million ha worldwide [START_REF] Albaugh | Eucalyptus and Water Use in South Africa[END_REF]. In Brazil, forest plantations cover 7.6 million ha, of which 72% are planted with high-productivity Eucalyptus clones (average annual increment of 40 m 3 ha -1 yr -1 of roundwood, ranging from 25 to 60 m 3 ha -1 yr -1 [START_REF] Gonçalves | Integrating genetic and silvicultural strategies to minimize abiotic and biotic constraints in Brazilian eucalypt plantations[END_REF]). Eucalyptus plantations in Brazil are generally managed in 6-7 years rotation, with canopy closure occurring within 2-3 years after planting. The relatively low susceptibility to pests and diseases, rapid growth and high productivity, adaptability to varying soil and climate, and adequate fiber quality for the industry explain the expansion of commercial Eucalyptus tropical plantations worldwide.

Water deficit, nutrient deficiency, soil type and compaction are the main drivers of Eucalyptus plantation functioning in Southern Brazil [START_REF] Ryan | Factors controlling Eucalyptus productivity: How water availability and stand structure alter production and carbon allocation[END_REF][START_REF] Stape | The Brazil Eucalyptus Potential Productivity Project: Influence of water, nutrients and stand uniformity on wood production[END_REF].

However, predicting how forests grow in response to soil and climate constraints, and determining their carbon storage capacities remains a key challenge for modelers. Forest productivity is driven by complex interactions and feedbacks among biological mechanisms. Furthermore, spatial variability in resources supplies, management, and characteristics of the genetic tree material critically influence forest productivity.

Ecophysiological process-based models that simulate water and carbon fluxes in forests proved to be useful tools to formalize biophysical hypotheses on forest functioning and to test for the importance of environmental drivers on productivity. Over the last two decades, a range of process-based models was developed, varying in resolution, complexity, generality, and applicability [START_REF] Mäkelä | Process-based models for forest ecosystem management: current state of the art and challenges for practical implementation[END_REF][START_REF] Battaglia | CABALA: a linked carbon, water and nitrogen model of forest growth for silvicultural decision support[END_REF]Corbeels et al., 2005a;[START_REF] Dufrêne | Modelling carbon and water cycles in a beech forest: Part I: Model description and uncertainty analysis on modelled NEE[END_REF][START_REF] Marsden | Modifying the G'DAY process-based model to simulate the spatial variability of Eucalyptus plantation growth on deep tropical soils[END_REF]. Some of these models were developed for research purposes to understand and quantify carbon and water cycling at fine time-scale [START_REF] Dufrêne | Modelling carbon and water cycles in a beech forest: Part I: Model description and uncertainty analysis on modelled NEE[END_REF], others were constructed as management tools, in close collaboration with intended end-users [START_REF] Landsberg | A generalised model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning[END_REF][START_REF] Sands | Application of process-based models to forest management: experience with PROMOD, a simple plantation productivity model[END_REF]Almeida et al., 2004b).

Importantly, wood productivity relates on both the net amount of carbon (C) sequestered by trees (net C balance) and the way this C is allocated among tree organs.

The modelling of Eucalyptus plantations at regional scale is of critical economic importance and has been the focus of a growing body of studies (e.g., Almeida et al. (2004a). [START_REF] Marsden | Modifying the G'DAY process-based model to simulate the spatial variability of Eucalyptus plantation growth on deep tropical soils[END_REF] modified the Generic Decomposition And Yield model (G'DAY) to simulate the productivity of 16 Eucalyptus plantation stands of the same genotype in São Paulo region, Brazil. It was shown that soil water holding capacity explains a large part of spatial variability in tree growth rate and biomass productivity in this area. This study also highlighted two important limitations that currently hinder regional modelling using G'DAY.

First, the impact of water availability on the C allocation toward root growth, leaf growth and litterfall production is not correctly represented in the model. Because leaf and root areas are critical C and water exchange surfaces, further improvements of the G'DAY C allocation scheme is required for accurate simulations of carbon and water fluxes and productivity along soil water gradients. Second, the diversity of species and hybrid Eucalyptus genotypes used in Brazil's plantations was shown to be challenging for large-scale modelling [START_REF] Almeida | Mapping the effect of spatial and temporal variation in climate and soils on Eucalyptus plantation production with 3-PG, a process-based growth model[END_REF]. Eucalyptus materials may indeed strongly differ in the control of several key processes, such as photosynthesis and light use efficiency [START_REF] Warrier | Gas exchange characteristics in Eucalyptus clones[END_REF]le Maire et al., accepted) and carbon allocation [START_REF] Ngugi | Dry matter production and allocation in Eucalyptus cloeziana and Eucalyptus argophloia seedlings in response to soil water deficits[END_REF]. However, a detailed parameterization of process-based models for each genotype from field measurements is today out of reach, because physiological and C partitioning measurements along plantation rotation are missing for most of the numerous cultivated genotypes. Previous studies therefore generally used a unique model parameterization over a whole region, either by using generic parameterization at the level of plant functional type (e.g., Almeida et al., 2004a), or by calibrating the model based on some genotype-specific measurements (Almeida et al., 2004a;[START_REF] Gonzalez-Benecke | Regional validation and improved parameterization of the 3-PG model for Pinus taeda stands[END_REF]. Calibration can be conducted by using independent measurements to parameterize the equations of the model, or by inverting the model to constrain the values of these parameters [START_REF] Gonzalez-Benecke | Regional validation and improved parameterization of the 3-PG model for Pinus taeda stands[END_REF]. This last option is considered more reliable for application purposes, since the model simulations are constrained within the range of a plausible domain, even though equifinality of the simulations (i.e., the possibility of 'getting the right answers for the wrong reason') may lead to a wrong interpretation of the results. Although using a unique parameterization for regional modelling proved useful to capture the response of productivity to large environmental gradient, it creates an additional uncertainty that need to be quantified if the simulations are to be used for management purposes.

The overall objectives of the present study was 1) to predict temporal and spatial variations of stem biomass production in Eucalyptus plantations in the São Paulo state of Brazil using the G'DAY model 2) to use the validated model to test for the effect of genotype, soil and climate on plantation wood productivity at regional scale. Efforts were done to improve the ability of the model in capturing the responses of forest functioning and C allocation to water stress. We used a multi-site calibration approach to provide an end-toend calibration and application scheme from local to regional scale. G'DAY was first optimized using detailed measurements from one experimental site. A subset of "genotypespecific" parameters was then calibrated at two other sites where physiological data were also available, and presenting differences in genotype and location. Finally a calibration and application of the model at large spatio-temporal scale was performed using a network of 1472 stands. Constrained simulations were then conducted using the final model to highlight the drivers of productivity at regional scale.

Materials and Methods

The four Eucalyptus datasets

The site of the first experimental dataset (DATASET 1) was a commercial plantation of Eucalyptus grandis, located at 22°58'04''S, 48°43'40''W, 750 m.a.s.l, planted in November 2009 at a 3x2 m spacing, and monitored continuously since then as part of the EUCFLUX project (http://www.ipef.br/eucflux/en/, [START_REF] Nouvellon | Water and energy fluxes above an Eucalyptus plantation in Brazil: environmental control and comparison with two eucalypt plantations in Congo[END_REF][START_REF] Nouvellon | Full-rotation carbon, water and energy fluxes in a tropical eucalypt plantation[END_REF]). The average annual precipitation was 1540 mm from 2008 to 2016, with an average temperature of 19.3°C, and a wet hot summer from October to May and dry cold winter from June to September. A detailed description of the site is given in [START_REF] Campoe | Stem production, light absorption and light use efficiency between dominant and non-dominant trees of Eucalyptus grandis across a productivity gradient in Brazil[END_REF] and [START_REF] Christina | Sensitivity and uncertainty analysis of the carbon and water fluxes at the tree scale in Eucalyptus plantations using a metamodeling approach[END_REF][START_REF] Christina | Importance of deep water uptake in tropical eucalypt forest[END_REF]Measurements included volumetric soil moisture content, tree height, leaf area index (LAI), H2O and CO2 gas exchanges between the ecosystem and the atmosphere, leaf, bark and branch litterfall, biomass of all tree compartments, leaf photosynthesis. Daily net ecosystem exchange (NEE) and evapotranspiration (ET) were obtained from a flux-tower, using the Eddy-covariance method [START_REF] Christina | Importance of deep water uptake in tropical eucalypt forest[END_REF][START_REF] Nouvellon | Full-rotation carbon, water and energy fluxes in a tropical eucalypt plantation[END_REF][START_REF] Vezy | Measuring and modelling energy partitioning in canopies of varying complexity using MAESPA model[END_REF]. Gross primary productivity (GPP) was estimated from NEE and meteorological data using the standard computation from [START_REF] Reichstein | On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm[END_REF]. Soil water content was measured using calibrated CS616 probes (Campbell Scientific Inc., Logan, UT, USA) inserted at 0.15, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 m soil depths. Soil is sandy loam with an average maximum plant available water (PAWmax) of 92 mm m -1 [START_REF] Christina | Importance of deep water uptake in tropical eucalypt forest[END_REF]. Leaf area index was measured by combining destructive measurements and dendrometric inventories, and leaf and bark litterfall was collected monthly using 48 litterfall traps of 0.52 by 0.52 m placed in the field following a voronoï scheme to sample every distance to trees. Branch litterfall was collected monthly on four 6.6 m 2 area. Litter samples were dried, weighted, and values of dry mass by unit area were computed (gDM/m 2 ).

The site of the second experimental dataset (DATASET 2), located at Santa Rita municipality (21°35' 48'' S, 47°36' 0'' W, 761 m.a.s.l.) was planted on April 2004 with a clone of E. grandis × E. urophylla, at a 3 x 2.4 m spacing. The average annual precipitation is 1500 mm, and mean temperature was 22.1°C. Data collection started two years after planting on March 2006until March 2008. Collected data included measurements of soil water content, LAI, biomass of leaves and stem, fully described in [START_REF] Cabral | The energy and water balance of a Eucalyptus plantation in southeast Brazil[END_REF] and [START_REF] Cabral | Fluxes of CO2 above a plantation of Eucalyptus in southeast Brazil[END_REF]. Soil moisture content was measured at 0.1, 0.3, 0.5, 0.7, and 1 m soil depth. The soil was sandy, with PAWmax of 50 mm m -1 [START_REF] Cabral | Fluxes of CO2 above a plantation of Eucalyptus in southeast Brazil[END_REF]. Daily LAI values was computed from the reflection coefficient of photosynthetically active radiation with a linear relationship previously calibrated against destructive LAI measurements. NEE and ET were measured through the Eddy covariance method by using a flux tower at a height of 27 m, together with meteorological variables, including photosynthetically active radiation measurements [START_REF] Cabral | The energy and water balance of a Eucalyptus plantation in southeast Brazil[END_REF][START_REF] Cabral | Fluxes of CO2 above a plantation of Eucalyptus in southeast Brazil[END_REF].

The third experimental stand (DATASET 3) was located at Itatinga (23°02'28'' S, 48°37'33 W, 850 m.a.s.l) and planted with Eucalyptus grandis seedlings on April 2004. Data collection started at planting and continued until the end of the rotation, which included soil water, biomass of leaves, stem, roots, leaf area index, leaf, branch and bark litterfall. Soil water content was measured using TDR probes installed at various soil depth down to 10 m depth [START_REF] Laclau | Organic residue mass at planting is an excellent predictor of tree growth in Eucalyptus plantations established on a sandy tropical soil[END_REF]. The soil was sandy loam, with PAWmax of 110 mm m -1 [START_REF] Maquere | Dynamics of Mineral Elements under a Fast-growing Eucalyptus Plantation in Brazil[END_REF][START_REF] Marsden | Modifying the G'DAY process-based model to simulate the spatial variability of Eucalyptus plantation growth on deep tropical soils[END_REF]. This site was used in a previous application of G'DAY model [START_REF] Marsden | Modifying the G'DAY process-based model to simulate the spatial variability of Eucalyptus plantation growth on deep tropical soils[END_REF].

The fourth dataset (DATASET 4) was made of measurements of 1472 Eucalyptus species stand polygons comprised in an area ranging from 22°33' S to 20°50' S and 48°14' W to 46°49'W (i.e., spread on an 183 × 151 km area). The 1472 polygons have an average area of 29 ha, and were planted with various genotypes of Eucalyptus (there is almost no coppicing practice in this dataset). Soil was highly variable in sand and clay content and therefore showed a large range of PAWmax (see section 2.2). Plant height and stem biomass were obtained from field inventories performed at two to three inventory dates (ages) for each polygon, between 2000 and 2012. The dataset and biomass calculations are presented in le Maire et al. (2011a), [START_REF] Baghdadi | Testing different methods of forest height and aboveground biomass estimations from ICESat/GLAS data in eucalyptus plantations in Brazil[END_REF] and [START_REF] Baghdadi | Evaluation of ALOS/PALSAR L-Band Data for the Estimation of Eucalyptus Plantations Aboveground Biomass in Brazil. Selected Topics in Applied Earth Observations and Remote Sensing[END_REF].

2.2. Regional scale soil and meteorological datasets Application of the model at regional scale required reliable sources of weather and soil data. We obtained the gridded weather data from the open-access dataset for daily meteorological variables in Brazil (1980Brazil ( -2013) ) [START_REF] Xavier | Daily gridded meteorological variables in Brazil (1980-2013)[END_REF]. This dataset provides high resolution grids (0.25° by 0.25°) of daily precipitation, evapotranspiration, maximum and minimum temperature, solar radiation, relative humidity, and wind speed developed by the CLIMA research team using ground-based weather stations in Brazil, operated by federal (INMET, ANA) and state (DAEE for São Paulo) agencies. This gridded weather dataset was compared with the three DATASET 1, 2 and 3 for solar radiation, maximum and minimum temperature, and precipitation. Gridded solar radiation (Supplementary Figs. 1A, 1B, and1C) and maximum and minimum temperature recorded at the three sites (Supplementary Figs 1D, 1E, and 1F and Figs. 2A, 2B, and2C) matched very well the local measurement data. For precipitations, the errors were larger, but the order of magnitude was respected (Supplementary Figs. 3 and4). Overall, the gridded weather data was considered a reliable source for the regional simulations.

Several sources of soil data were tested to determine the dataset that best captured the variability of soil texture, and possibly the water retention properties subsequently used in the model. The digitized 1:250.000 scale soil map of Sao Paulo state from [START_REF] Rossi | Mapa pedológico do Estado de São Paulo: revisado e ampliado[END_REF] was the most detailed one and included a soil texture attribute (Fig. 1). We overlay this map with about 600 soil profiles measured at various locations across São Paulo state obtained from the forest company, the Brazilian Agricultural Research Corporation (Embrapa) and the Geological Service of Brazil (CPRM). These soil profiles generally were to 2 m deep and included soil texture properties, soil organic matter content and bulk density. This information was used within Tomasella's pedotransfer function (level 3) to estimate the water retention parameters [START_REF] Tomasella | Pedotransfer functions for the estimation of soil water retention in Brazilian soils[END_REF], which enables the computation of the maximum plant available water content (PAWmax) on each soil profile, with a per meter depth unit (mmwater msoil -1 ). The average PAWmax value per soil textural class extracted from Sao Paulo soil map were computed, together with their standard deviation, and further extrapolated to the entire map (Supplementary Fig. 5).

G'DAY model

In the present study, the G'DAY model [START_REF] Comins | Long-Term Response of Nutrient-Limited Forests to CO"2 Enrichment; Equilibrium Behavior of Plant-Soil Models[END_REF]Corbeels et al., 2005a;Corbeels et al., 2005b;[START_REF] Marsden | Modifying the G'DAY process-based model to simulate the spatial variability of Eucalyptus plantation growth on deep tropical soils[END_REF] was used to simulate water and carbon budgets of Eucalyptus plantations at a number of experimental sites and commercial stands in the state of São Paulo, Brazil. The G'DAY model is an ecophysiological process-based model, functioning at a daily time-step, which uses minimum and maximum air temperatures, precipitation, vapor pressure deficit and solar radiation as daily weather inputs. G'DAY simulates the water and carbon fluxes between the environment and a number of soil and tree biomass pools. The ecophysiological processes in the model are represented by several submodels of plant production, soil water balance, and litter decomposition. [START_REF] Marsden | Modifying the G'DAY process-based model to simulate the spatial variability of Eucalyptus plantation growth on deep tropical soils[END_REF] modified the soil water balance in order to reflect the dynamics of moisture content in three layers: litter, top soil, and rooting zone. The maximum PAW in the litter layer was modified to increase linearly as the C mass of the forest floor increased, whereas in the rooting zone the total maximum PAW increased with rooting depth during the rotation. Here, we used the version presented in details in [START_REF] Marsden | Modifying the G'DAY process-based model to simulate the spatial variability of Eucalyptus plantation growth on deep tropical soils[END_REF], but added a new modification of the C allocation scheme to better consider the impact of environmental constraints such as water availability on the partitioning of C between shoots and roots and to enhance the model capability of capturing seasonal leaf area index (LAI) variations.

Allocation fraction of the net primary production (NPP) to the different organs, is obtained following a "goal-seeking" scheme, where a constant allocation value was modulated in function of a target compartment biomass (see Corbeels et al. (2005a) and [START_REF] Marsden | Modifying the G'DAY process-based model to simulate the spatial variability of Eucalyptus plantation growth on deep tropical soils[END_REF]), and mortality is computed as turnover rates. We slightly modified this C allocation scheme as following:

i)

The proportion of carbon allocated to fine roots compared to foliage production was based on an assumption of a higher allocation to root when the soil is dry [START_REF] Landsberg | Physiological ecology of forest production: principles, processes and models[END_REF];

ii)

The leaf turnover was set to increase when the soil was dry (severe and prolonged stress episodes) and was positively correlated to the production of new leaves, as observed in litterfall measurements [START_REF] Pook | Long-term variation of litter fall, canopy leaf area and flowering in a Eucalyptus maculata forest on the south coast of New South Wales[END_REF];

iii)

The "target" value of leaf area, which came from a height-dependent target ratio between leaf area and sapwood cross-sectional area (Corbeels et al., 2005a;Corbeels et al., 2005b;[START_REF] Marsden | Modifying the G'DAY process-based model to simulate the spatial variability of Eucalyptus plantation growth on deep tropical soils[END_REF], was previously reached by modulating the carbon allocation to leaves. Within the new allocation scheme, we kept the same principle but modulated the leaf turnover ratio in order to match this algometrical constraint, while C allocation to new leaves was modulated by the equation described in i)

Use of water from very deep layers was shown to sustain transpiration of Eucalyptus plantation during extreme drought events [START_REF] Christina | Importance of deep water uptake in tropical eucalypt forest[END_REF], permitted by the fast soil exploration of Eucalyptus trees along the rotation and the ability for the root to access soil water table . A minimum access to water (i.e, a minimum value for soil PAW), was therefore included in the model in order to avoid a total transpiration break and reproduce the effect of deep soil water use on tree functioning observed in the field. Such minimum soil PAW was however very low, and correspond to high water stress effects on trees (reduced transpiration and photosynthesis).

Model Calibration

Calibration on the detailed DATASET 1

The G'Day model has more than 200 parameters, among which physical constants, parameters linked to the nitrogen cycle that were not considered in the present study, and a subset of plant and soil-specific parameters. The entire set of these 57 plant and soil-specific parameters, described in Supplementary Table 1, were optimized against DATASET 1 using the Hooke-Jeeves algorithm for derivative-free optimization in the R package dfoptim [START_REF] Varadhan | dfoptim: Derivative-Free Optimization[END_REF]. This algorithm minimizes the residual sum of squares between observed and predicted data while allowing bond constraints on the parameters. These bonds, given in Supplementary Table 1, were set based on the literature and previous work on the model [START_REF] Marsden | Modifying the G'DAY process-based model to simulate the spatial variability of Eucalyptus plantation growth on deep tropical soils[END_REF]. The calibration was done using all the following measurements, after normalization of their residual sum of squares: plant height (m), LAI, C mass of leaves (kg C m -2 ), leaves, branches, and bark litter fall (kg C m -2 ), stem biomass (kg C m -2 ), PAW (mm), evapotranspiration (ET, mm day -1 ), net ecosystem exchange (NEE, g C m -2 day -1 ), and gross primary photosynthesis (GPP, g C m -2 day -1 ).

Calibration on the DATASET 2 and 3

Only a subset of the parameters were re-calibrated (Supplementary Table 1) using DATASET Table 1). When little information was available on literature, larger bounds of 30% was given. On these 2 datasets, local climate information from meteorological station and fixed soil parameters were used. The recalibration was done independently for DATASET 2 and 3.

For DATASET 2, the calibration was done on LAI, C mass of stem, NEE, and ET, and for DATASET 3 on LAI, C mass of leaves and stem, cumulative leaves, branches, bark litterfall, and PAW.

Calibration on DATASET 4

Genotypes having fewer than 30 stands on DATASET 4, and therefore with little representativity for model calibration, were not used. Eleven genotypes were kept, which included many common genotypes used in fast growing Eucalyptus plantation in this region of Brazil at that time, including E. grandis and E urophylla x grandis materials. The dataset was then partitioned into two equal halves for each genotype, one for calibrating the parameters of the model and the other for evaluating the model (2*736 stands). Among the 57 parameters of G'DAY calibrated on DATASET 1, only 22 were calibrated for each genotype (Supplementary Table 1) using DATASET 4, with the same bounds defined for DATASET 2 and 3. Indeed, 21 parameters re-calibrated on DATASET 2 and 3 out of the 43 genotypespecific parameters were found to be hardly variable between the 3 first datasets, or to have no influence on simulated biomass. Values of the parameters obtained on DATASET 1 were used as generic values for parameters not recalibrated. For DATASET 4, the calibration was performed on two measured variables: plant height (m) and stem biomass (kg C m -2 ) measured at two to three inventory dates for each polygon, but on the high range of soil and climate conditions found in this area. In these simulations, the planting date was prescribed as its observed value recorded in DATASET 4. The climate and soil information from the gridded inputs (see section 2.2) was used for each polygon. Model performance was evaluated by R 2 , RMSE, relative RMSE, and Nash-Sutcliffe efficiency (NSE) statistics using modeval function in the R package sirad [START_REF] Bojanowski | Package 'sirad[END_REF]. Because some stands were measured two or three times during their growth, some of the observations within the evaluation set were not statistically independent. To avoid this eventual issue, the statistics described above were computed on a subset where only one measurement per polygon was randomly selected. This random selection was repeated a thousand time, to generate a distribution of the model performance statistics, from which the average and standard deviation was computed.

2.5. Analysis of the model results at regional scale Simulation of biomass at regional scale by G'DAY model was compared to a subset of 736 polygons set apart during the parameter calibration step. We quantified the importance of spatial variation of climate, soil properties, and genotype specific parameters sets, and their combinations, for the model accuracy. Eight different simulations were run on all 736 validation polygons (Table 2). The first set was a "base simulation", representing a scenario where no spatial variation of climate, soil or genotype was taken into account: one grid point of the meteorological gridded dataset, chosen at the center of DATASET 4, was selected for all polygons. Similarly, one value of PAWmax was used, based on the most frequent value in the total area. In base simulation, the set of parameters obtained on DATASET 1 after optimization was used for all stands, to eliminate the spatial variability due to change in genotypes (Table 1). The second, third, and fourth sets of simulations were performed with spatial variation of climate only, genotype specific parameters only, and soil data only, respectively. The fifth, sixth, and seventh sets of simulations included combinations of two of these spatial variability. The last set of simulations had all drivers varying spatially, which correspond to the more precise parameterization and input data in this study. The abilities of the eight sets of simulations to predicting stem biomass were compared using Tukey's HSD.

The absolute difference between stem biomass predictions and observations for each individual point (inventory date) in the 736 polygons were computed. The dependency of these absolute differences (a measure of the model performance) to Climate, Genotype and Soil was quantified using one way ANOVA. Climate, Genotype and Soil were considered as binary dummy variables, with 0 value when no variation was accounted for, and 1 when spatial variation of that characteristics was accounted for. All interactions among predictors were included. This ANOVA therefore allowed to quantify the importance of the spatial variability in climate, genotype and soil for the model performance in simulating wood productivity at regional scale. After calibration, the simulated LAI on DATASET 1 were fairly similar to observed LAI, except little underestimation in the third and fourth years of the rotation when observed LAI was measured at approximately 5.7 m²leaf/m²soil while simulated LAI was little less than 5 m²leaf/m²soil (Table 1 and Fig. 2B). The LAI curve followed the typical trend observed in commercial rotation of Eucalyptus trees, increasing to a maximum value three years after planting and then gradually decrease until harvest, with seasonal variations. Similar observations were reported for the prediction of C mass of leaves when the model slightly underestimated this C variable in the second and third years of the rotation, but matched the observed data in the following years (Table 1 and Fig. 2C).

Results

3

Leaf, branch, and bark litterfalls were also correctly predicted by the model for cumulative sum, with R 2 of 0.99 and NSE close to 1 (Table 1 and Fig. 2D). Plant height was correctly predicted throughout the rotation reaching maximum height of 28 m by the end of the rotation in 2016 (Fig. 2A).The model also simulated correctly the final stem biomass that reached 8 kg C m -2 at the end of the rotation (Table 1 and Fig. 2E).

The calibrated model correctly simulated PAW throughout the rotation, except for a little underestimation in early 2016 (Table 1). Figure 2F shows observed and predicted PAW and the green line in the graph represents PAWmax as used in the model. This value do not corresponds to the real PAWmax experienced by trees (that should be integrated down to the real depth of roots) but reflects a "functional" PAWmax. Indeed, PAW is used in many processes of G'DAY, and increasing PAWmax for very deep roots would results in very low fraction of available water during most of the rotation, even on rainy seasons. As explained in 2.3, we however included a minimum PAW in the model, which represents access to deep soil water. Such access is important for sustaining transpiration throughout the year, but represents a small amount of transpired water as shown in [START_REF] Christina | Importance of deep water uptake in tropical eucalypt forest[END_REF] and confirmed in these simulations.

Simulation of ET correctly matched the observed data throughout the rotation, except a small underestimation in the second half of 2013 and early 2014 (Fig. 3A). The goodness of fit statistics of that variable indicated acceptable level of agreement between observed and simulated data (Table 1). The model reasonably simulated the NEE with NSE of 0.75 and R 2 of 0.75 (Fig. 3B), but most importantly showed good seasonality. For GPP, simulated data followed a temporal dynamic similar to observed data, which show peaks every summer season along the rotation (Fig. 3C). There was an underestimation of GPP by the model in 2012-2013 which could be attributed to the underestimation of LAI during this period. Other research is ongoing to provide more confidence on these estimates. Nevertheless, observed and simulated daily GPP had R 2 of 0.92 with NSE of 0.84 and RRMSE of 16.5% (Table 1), and well simulated seasonality.

A subset of parameters was re-calibrated using DATASET 2 and 3. The model reasonably simulated LAI and stem biomass of DATASET 2 with average LAI of 3 m²leaf/m²soil and C mass of stem of 3 kg C m -2 in the beginning of 2008 (Figs. 4A and4B).

Note that LAI at this site was not measured using destructive sampling but was an estimation based on reflected light, which may result in higher uncertainty. This site has a sandy soil with average PAWmax of 50 mm m -1 of soil which could explain the lower values of LAI and biomass yield at this site compared to DATASET 1. In addition, low precipitation in March to November of 2006 and 2007 resulted in reduced LAI (Supplementary Fig. 3D andE). The model was able to capture the seasonal variations of NEE (Fig. 4C) and ET (Fig. 4D), except for underestimation of NEE at the beginning of 2007. The high ET values in summer 2007 were most likely due to the high precipitation (about 1000 mm) observed at that time (Supplementary Fig. 3D). The model predicted ET with R 2 of 0.76 and RMSE of 1.08 (mm day -1 ) and NEE with R 2 of 0.28 and RMSE of 1.7 (g C m -2 day -1 ) (Table 1). To reach this result, only some of the re-calibrated parameters were significantly changed such as those controlling the relationship between stem biomass and height (Ht0 and Htpower), turnover rate (Bfall), stomatal conductance (Fs1 and Fs2) and photosynthesis (Jref) (Supplementary Table 1).

Results at DATASET 3 included the model simulations of LAI, carbon mass of leaves and stem, litterfall, and PAW (Fig. 5). Overall, after calibration, the model correctly simulated LAI and C mass of leaves throughout the 6-year rotation (Table 1), except for a little underestimation of LAI in 2006 and overestimation of C mass of leaves in 2007 (Fig. 5A andB). Other plant compartments of leaves, branches, and bark litterfall were also well simulated with R 2 ranging from 0.99 to 0.80 and NSE close to 1 (Fig. 5C) as well as C mass of stem with RMSE of 0.33 (kg C m -2 ), RRMSE of 14.7% and NSE close to 1 (Fig. 5D).

Simulated and observed PAW were comparable with R 2 of 0.71 and NSE of 0.69, apart from some underestimation in the first half of 2006 which could be related to the underestimation of LAI in this period (Figs. 5A andG). Similarly to DATASET 2, to reach these results, only some of the re-calibrated parameters were significantly changed such as parameters of growth efficiency (GrowthEff), relationship between stomatal conductance and VPD (gmax and Fs2), relationship between stem biomass and height (Ht0 and Htpower), and turnover rates (Bdecay and Bfall) (Supplementary Table 1).

Finally, based on the results from DATASET 1, 2 and 3, we concluded that: 1) the model is flexible and generic enough to represent the main processes controlling carbon and water balance of these eucalypt plantations; 2) the use of large dataset including many field measured variables, and the choice of adequate bounds for parameters based on field measurements and literature, allow to constrain enough the model in a calibration procedure;

3) only a small subset of 22 parameters is necessary to be modified in the model to be able to simulate very different plantations, in terms of climate conditions, soils and genotypes. The procedure was therefore considered suitable, and was extrapolated to DATASET 4, which have large range of spatial and genotype variation. In that case, the model was calibrated at genotype-scale on biomass and stand height measurements.

3.1. Model calibration and validation at regional scale (DATASET 4)

The large scale spatio-temporal DATASET 4 included measurements of plant height and stem biomass on commercial plantations. Calibrated values from DATASET 1 were used for the regional simulation along with genotype specific parameters that were calibrated for each genotype in the regional simulations (Supplementary Table 1). Plant height was fairly well simulated with some differences in the level of precision among genotypes. For instance, prediction of plant height on independent stands of genotypes C, D, F, and K had greater correlation with observed data and lower RMSE values than genotypes B, H, and J where the level of correlation was lower and the RMSE was greater (Supplementary Figs 6,7,[START_REF] Le Maire | Light absorption, light use efficiency and productivity of 16 contrasted genotypes of several Eucalyptus species along a 6-year rotation in Brazil[END_REF].

Across all genotypes, simulated and observed height were correlated by 93% with RMSE of 1.83 m and RRMSE of 10.47% and NSE 0.86 (Fig. 6).

Similar results was obtained for biomass yield as the model performances were higher for genotypes C, F, and G with R 2 of 0.91-0.93 and RMSE of 0.63-0.73 kg C m -2 , while for some other genotypes like B the R 2 was ca. 0.80 and RMSE was 0.84 kg C m -2 (Supplementary Figs. 6, 7, and 8). Nonetheless, the overall correlation between observed and simulated biomass was 89% and RMSE was 0.75 kg C m -2 (Fig. 6, All genotypes), on stands independent from model calibration.

Impact of climate, genotype specific parameters, soil data, and their combinations on model error estimation

There was significant differences among simulation scenarios for stem biomass prediction (Figure 7 and Table 3). Prediction of stem biomass and plant height without taking into account the spatial variation of climate, soils and genotypes had the greatest error.

Improvement of the model to simulate stemwood biomass was drastically improved when including differences of Genotype, then Climate and then Soil. Effect of Genotype was even further improved when climate or soil was considered together (Figure 7). The results of the ANOVA similarly showed significant individual and two-way interaction effects of climate, soil, and genotype on the accuracy of stem biomass prediction (Table 3). Fisher values confirmed the dominant effect of Genotype calibration, followed by Climate and Soil.

Interactions reached a high level, in particular the Soil x Genotype interactions.

Discussion

Parameterization and evaluation of the model

We used a slightly modified version of the G'DAY allocation model based on the "balanced growth hypothesis" that postulates that plants invest more carbon to roots when the limiting factor for growth is water or nutrients, and reduce the foliage size as drought avoiding mechanism [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF]. In addition, the influence of water stress on litterfall production and LAI dynamics was also taken into consideration, as observed on long term measurements [START_REF] Pook | Long-term variation of litter fall, canopy leaf area and flowering in a Eucalyptus maculata forest on the south coast of New South Wales[END_REF]. Evaluation of the model showed that the incorporation of the modified carbon scheme enhanced the model's ability to capture the decline of LAI as a result of soil drought, and consequently improved the simulation of leaf biomass.

In the DATASET 1 calibration, the final values of some parameters were restrained by the lower or upper bounds of prior distributions (Supplementary Table 1). In the present work the objective was mainly to get plausible prior probability distribution of parameters to be further included on larger scale applications. These parameter values, even when constrained on bound limits, are still plausible and gave good results at this site with respect to carbon and water measurements. On other datasets, parameters also were constrained by bounds, but this was necessary because of the fewer measurements available for calibration.

The overall picture and the good results obtained on independent data on DATASET 4

showed that these limits were sufficient for a first trial. More spatial information on vegetation characteristics, such as variables estimated from remote sensing (le Maire et al., 2011b;[START_REF] Baghdadi | Testing different methods of forest height and aboveground biomass estimations from ICESat/GLAS data in eucalyptus plantations in Brazil[END_REF][START_REF] Baghdadi | Evaluation of ALOS/PALSAR L-Band Data for the Estimation of Eucalyptus Plantations Aboveground Biomass in Brazil. Selected Topics in Applied Earth Observations and Remote Sensing[END_REF], and to more datasets could allow to improve data assimilation, with a refinement of the distribution of calibrated parameters probability distribution.

Widely distributed Eucalyptus species show high levels of genetic diversity in light absorption, gross primary production and differences in C allocation which determine together genotype productivity [START_REF] Aspinwall | Photosynthesis and carbon allocation are both important predictors of genotype productivity responses to elevated CO 2 in Eucalyptus camaldulensis[END_REF]le Maire et al., accepted). In the current study, we calibrated 22 parameters to address differences among genotype.

Parameters related to leaf and branch mortality and turnover rate, C allocation, and stomatal conductance in response to climatic factors showed the greatest variations among genotypes (Supplementary Table 1). Previous research on similar Brazilian Eucalyptus plantations has reported variations in photosynthetic capacity between E. grandis × urophylla hybrids in comparison with E. grandis impacting productivity (Almeida et al., 2004a). Leaf photosynthetic parameters of Eucalyptus plantations can also vary considerably within clones of the same species [START_REF] Shem | Gas exchange responses of Eucalyptus, C. africana and G. robusta to varying soil moisture content in semi-arid (Thika) Kenya[END_REF][START_REF] Warrier | Gas exchange characteristics in Eucalyptus clones[END_REF]le Maire et al., accepted). Other reports suggested that C allocation to roots and other C sinks may be dominant drivers of genotypic variation in productivity in response to environmental factors [START_REF] Resco De Dios | Intraspecific variation in juvenile tree growth under elevated CO2 alone and with O3: a meta-analysis[END_REF]. This results in large differences in light absorption and light use efficiencies and growth along rotations (le Maire et al., accepted).

Drivers of modelled spatial variability of stand biomass

Several factors influence biomass production in Eucalyptus plantations among which the climate, soil type, genetic material, and their interactions. One of the main advantages of process-based models, following proper parameterization, is the possibility to identify, quantify and disentangle the influence of spatial variation in soil, climate, and genotypes on productivity. Among these factors, genotype specific parameters were the most important for the accuracy of stem biomass predictions (Fig. 7 and Table 3). Genotype-specific calibration of ecophysiological models is critical for accurate characterization of stem biomass production among stands planted with different hybrids of Eucalyptus. In the present study, genotype-specific parameters were obtained by constrained optimization on only two variables, trunk biomass and height. In the future, other measurements could help improving the estimation of these parameters, such as data issued from remote sensing (le Maire, 2018).

The spatial variability of climate was the second factor in importance for prediction accuracy of stem biomass (Fig. 7 and Table 3). G'DAY was sensitive to varying climate scenarios: a reduced error in stem biomass and plant height predictions of about 20% was obtained when taking into account the local climate compared with holding the climate scenario constant for all stands (from a gridded dataset). While this result is expected on large climate gradient ("regional scale"), we showed here that even on reduced area ("regional scale") of 183 × 151 km the climate has a major importance. Rainfall and drought occurrence was the most important climatic factor constraining the growth of Eucalyptus in this area, as observed in other studies [START_REF] Mummery | Significance of rainfall distribution in predicting eucalypt plantation growth, management options, and risk assessment using the process-based model CABALA[END_REF][START_REF] Stape | Eucalyptus production and the supply, use and efficiency of use of water, light and nitrogen across a geographic gradient in Brazil[END_REF][START_REF] Whitehead | Physiological regulation of productivity and water use in Eucalyptus: a review[END_REF]. In addition, research conducted in tropical Brazilian forest, near the Atlantic Coast, reported 70 to 110% annual variations in production in response to soil and air humidity [START_REF] Almeida | Mapping the effect of spatial and temporal variation in climate and soils on Eucalyptus plantation production with 3-PG, a process-based growth model[END_REF]. Accuracy of the input rainfall dataset could be ameliorated using denser networks of pluviometers.

Among the three drivers of modelled spatial variability of stand biomass, varying the soil type parameter was the least important driver, reducing the error in prediction by only 15% (Fig. 7 and Table 3). Nevertheless, the soil x climate or soil x genotype interaction effects reduced the error in stem biomass prediction in a larger way than climate alone or genotype alone, indicating an improvement of the simulations when considering variations in soil water holding capacity, in the studied range. A previous study aiming at simulating the spatial variability of Eucalyptus plantation in São Paulo state using the G'DAY model reported significant improvement in stem biomass prediction when using stand-specific PAW values compared with holding it constant for all stands [START_REF] Marsden | Modifying the G'DAY process-based model to simulate the spatial variability of Eucalyptus plantation growth on deep tropical soils[END_REF]. The authors attributed the higher performances of simulations using stand-specific PAW values to the linkage of maximum PAW with tree height. The current version of the model has shown improvements in simulation of stem biomass using stand-specific PAW values as well; however, it did not show much improvements when comparing all drivers variable (scenario 8) with varying climate and genotype (scenario 5) (Figure 7). This indicates that further improvements in the model responses to the climate x genotype x soil interactions is needed. Improvement should also be done on spatial parameterization of the maximum plant available water (PAWmax).

In the current study, this gridded map has several limitations: 1) the texture class associated to each soil type of the soil map was not always precise, in particular for some complex soils;

2) associating a unique soil textural class in function of soil types do not allow to represent eventual variability within a soil type; 3) the association of a texture class to a single PAWmax value does not reflect the variability visible on Supplementary Figure 5; 4) differences in soil depth is not taken into account: here all stands are supposed to have access to the first ~3 meters of soil, and keep having minimal water access along their rotation.

However, even with these limitations, the use of these gridded estimates of soil PAWmax did improved the results compared to using a single value at all location.

Conclusions

A modified version of the G'DAY model was able to simulate seasonal variations in growth and the exchange of key C and H2O variables between the ecosystem and the atmosphere along complete commercial rotations of Eucalyptus plantations. Application of the model at the regional scale showed reasonable level of accuracy in the simulation of stem biomass and plant height. The main drivers of spatial variability in simulated stem productivity was the genetic differences among genotypes followed by climatic and soil variables. This work will benefit in the future from other data sources, such as remotesensing, allowing to further constrain the mechanisms embodied in ecophysiological models at various temporal and spatial scales.

Acknowledgements

This project was funded by the FAPESP-Microsoft Research the project "Characterizing and Predicting Biomass Production in Sugarcane and Eucalyptus Plantations in Brazil" (SEMP) of the Sao Eight different simulations were run on all 736 validation polygons (see description of the scenarios in Table 2). Scenarios with the same letter are not significantly different.

2

  and 3. A subset of 43 parameters among the 57 was selected based on expert knowledge on their possible variations among genotypes: parameters which are known to be stables among genotypes, or showing very low sensitivity of stem biomass evolution through time, were kept at the values obtained after calibration on DATASET 1. The 43 parameters are mainly linked to carbon assimilation and allocation, allometry, wood density properties, water balance submodels. Parameters obtained after calibration on DATASET 1 were used as initial values of those re-calibrated parameters. The parameters re-calibrated locally were constrained by new bounds values, based on the calibrated parameters obtained on DATASET 1. Indeed, these parameters are not expected to vary much among genotypes or with management practices, and therefore the values obtained on DATASET 1 were considered reliable priors. Bounds of +/-10 %, 20% or 30 % around this value was chosen in function of the current knowledge on physiological differences among clones (Supplementary

Figure 1 .

 1 Figure 1. Soil map of São Paulo state[START_REF] Rossi | Mapa pedológico do Estado de São Paulo: revisado e ampliado[END_REF] showing the location of the four datasets.

Figure 2 .

 2 Figure 2. Comparison of measured (OBS, black line/symbol) and simulated (SIM, red

Figure 3 .

 3 Figure 3. Comparison of measured and simulated variables on DATASET 1: (A)evapotranspiration (mm day -1 ), (B) net ecosystem exchange (g C m -2 day -1 ), and gross primary production (g C m -2 day -1 ).

Figure 4 .

 4 Figure 4. Comparison of measured and simulated variables on DATASET 2: (A) leaf area index (m 2 m -2 ), (B) C mass of stem (kg C m -2 ), (C) net ecosystem exchange (g C m -2 day -1 ), and (D) evapotranspiration (mm day -1 ).

Figure 5 .

 5 Figure 5. Comparison of measured and simulated variables on DATASET 3: (A) leaf area index (m 2 m -2 ), (B) C mass of leaves (kg C m -2 ), (C) cumulative leave, branch, and bark litterfall (kg C m -2 ), and (D) C mass of stem (kg C m -2 ), and plant available water (mm). The

Figure 6 .

 6 Figure 6. Measured vs. simulated plant height (m) and stem biomass (g C m -2 ) of all

Figure 7 .

 7 Figure 7. Quantifying the effects of including the spatial variation of climate, soil properties,

Table 1 .

 1 Model performance statistics after model calibration on three experimental stands (DATASET 1, 2, and 3) : coefficient of determination (R 2 ), root mean square error (RMSE), relative RMSE (RRMSE), and Nash-Sutcliffe efficiency (NSE).

	Variables	Corresponding Figure	R 2	RMSE	RRMSE	NSE
	DATASET 1
	Plant height	Fig. 2A	0.98	1.03 (m)	11.4%	0.98
	Leaf area index	Fig. 2B	0.94	0.42	13.1%	0.93
	C mass of leaves	Fig. 2C	0.85	0.03 (kg C m -2 )	21.7%	0.82
	Leaves litterfall	Fig. 2D	0.99	0.03 (kg C m -2 )	6.2%	0.99
	Branches litterfall	Fig. 2D	0.99	0.02 (kg C m -2 )	9.5%	0.98
	Bark litterfall	Fig. 2D	0.97	0.004 (kg C m -2 )	13.3%	0.96
	C mass of stem	Fig. 2E	0.99	0.33 (kg C m -2 )	10.2%	0.98
	Plant available water	Fig. 2F	0.86	38.6 (mm)	37.9%	0.81
	Evapotranspiration	Fig. 3A	0.71	0.91 (mm day -1 )	24.1%	0.66
	Net ecosystem exchange	Fig. 3B	0.75	1.64 (g C m -2 day -1 )	57.9%	0.74
	Gross primary production	Fig. 3C	0.85	1.5 (g C m -2 day -1 )	16.5%	0.84
	DATASET 2
	Leaf area index	Fig. 4A	0.31	0.28	9.8%	0.29
	C mass of stem	Fig. 4B	0.99	0.25 (kg C m -2 )	13.4%	0.96
	Net ecosystem exchange	Fig. 4C	0.28	1.7 (g C m -2 day -1 )	55.2%	0.24
	Evapotranspiration	Fig. 4D	0.76	1.08 (mm day -1 )	32.9%	0.72
	DATASET 3
	Leaf area index	Fig. 5A	0.95	0.23	7.6%	0.95
	C mass of leaves	Fig. 5B	0.93	0.02 (kg C m -2 )	16.9%	0.85
	Leaves litterfall	Fig. 5C	0.97	0.09 (kg C m -2 )	15.5%	0.92
	Branches litterfall	Fig. 5C	0.99	0.03 (kg C m -2 )	29.8%	0.94
	Bark litterfall	Fig. 5C	0.80	0.01 (kg C m -2 )	87.5%	0.73
	C mass of stem	Fig. 5D	0.99	0.33 (kg C m -2 )	14.7%	0.96
	Plant available water	Fig. 5E	0.71	46.6 (mm)	33.4%	0.69

Paulo Research Foundation (FAPESP-Microsoft Research n° 2014/50715-9). Ahmed Attia was supported by a post-doctoral grant 2017/00886-0 from Sao Paulo Research Foundation (FAPESP). The contribution of MG was supported by grant 2012/06933-6 from Sao Paulo Research Foundation (FAPESP). This project was funded by the EUCFLUX project (cooperative program with Arcelor Mittal, Cenibra, Bahia Specialty Cellulose, Duratex, Fibria, International Paper, Klabin, Suzano, and Vallourec Florestal), coordinated by the Forestry Science and research Institute -IPEF (https://www.ipef.br/). The experiment also partially benefited of fundings from Agence Nationale de la Recherche (MACACC project ANR-13-AGRO-0005, Viabilité et Adaptation des Ecosystèmes Productifs, Territoires et Ressources face aux Changements Globaux AGROBIOSPHERE 2013 program). The experimental site belongs to the SOERE F-ORE-T, which is supported annually by Ecofor, Allenvi and the French National Research Infrastructure ANAEE-F (http://www.anaeefrance.fr/fr/). We acknowledge support from the IN-SYLVA network. We are grateful to the staff at the Itatinga Experimental Station, in particular Rildo Moreira e Moreira (Esalq, USP) and Eder Araujo da Silva (http://www.floragroapoio.com.br) for their technical support.

Table 2:

The height different simulations scenarios at regional scale, taking into account or not the spatial variation in climate, Eucalyptus genotype and soil. When no climate variation were input, the meteorology from the central grid point was used. When no soil variation were input, the more frequent soil type was used. When no genotype variation were input, the genotype from DATASET1 was used.