Théophile Lohier

Franck Jabot
email: franck.jabot@irstea.fr

Driss Meziane

Bill Shipley

Peter B Reich

Guillaume Deffuant

Explaining ontogenetic shifts in root-shoot scaling with transient dynamics

Keywords: allometry, functional trait, grassland, model selection, optimal partitioning theory, transient dynamics

Backgrounds and Aims Simple models of herbaceous plant growth based on optimal partitioning theory predict, at steady state, an isometric relationship between shoot and root biomass during plant ontogeny: i.e. a constant root-shoot ratio. This prediction has received mixed empirical support, suggesting either that optimal partitioning is too coarse an assumption to model plant biomass allocation, or that additional processes need to be modelled to account for empirical findings within the optimal partitioning framework. Here simulations are used to compare quantitatively two potential explanations for observed nonisometric relationships, namely nutrient limitation during the experiments and initial developmental constraints.

• Methods A simple plant growth model was built, based on optimal partitioning theory combined with empirically measured plant functional traits. We assessed its ability to reproduce plant relative growth rate and final root weight ratio. Predicted root-shoot ratios during plant ontogeny were compared to experimental observations. The effects of nutrient limitation and initial developmental constraints on root-shoot ratios were then tested.

• Key Results The model was found to accurately reproduce overall plant growth patterns, but failed, in its simplest form, at explaining non-isometric growth trajectories. Both nutrient limitation and ontogenetic developmental constraints were further shown to cause transient dynamics resulting in a deviation from isometry. Nitrogen limitation alone was not sufficient to explain observed trajectories of most plant species. The inclusion of initial developmental constraints (fixed non-optimal initial root-shoot ratios) enabled the reproduction of the observed trajectories and were consistent with observed initial root-shot ratios.

• Conclusions This study highlights that considering transient dynamics enabled us to reconcile theoretical predictions based on optimal partitioning with empirically measured ontogenetic root-shoot allometries. The transient dynamics cannot be solely explained by nutrient limitation during the experiments, pointing to a likely role of initial developmental constraints in the observed non-isometric growth trajectories.

INTRODUCTION

Plants both fix atmospheric carbon in their leaves by photosynthesis and capture soil water and nutrients by their roots. The way these basic resources are allocated among different plant organs during plant growth is of utmost importance to understanding such basic ecological phenomena as competition between plants [START_REF] Tilman | Plant strategies and the dynamics and structure of plant communities[END_REF][START_REF] Grime | Plant strategies, vegetation processes, and ecosystem properties[END_REF], the global carbon cycle and the consequences of rising atmospheric carbon dioxide [START_REF] Hungate | On the fate of carbon in grasslands under carbon dioxide enrichment[END_REF][START_REF] Mokany | Critical analysis of root: shoot ratios in terrestrial biomes[END_REF][START_REF] Bonan | Forests and climate change: forcings, feedbacks, and the climate benefits of forests[END_REF]. Allocation of biomass between roots and shoots in plants has received much attention [START_REF] Poorter | Biomass allocation to leaves, stems and roots: meta-analyses of interspecific variation and environmental control[END_REF]. This allocation strongly depends on environmental conditions [START_REF] Chapin | The mineral nutrition of wild plants[END_REF][START_REF] Poorter | Biomass allocation to leaves, stems and roots: meta-analyses of interspecific variation and environmental control[END_REF], so that it constitutes a major difficulty for plant growth modelling [START_REF] Thornley | Shoot: root allocation with respect to c, n and p: an investigation and comparison of resistance and teleonomic models[END_REF][START_REF] Roux | Carbon-based models of individual tree growth: a critical appraisal[END_REF].

Many approaches have been used to model allocation in plants [START_REF] Génard | Carbon allocation in fruit trees: from theory to modelling[END_REF][START_REF] Franklin | Modeling carbon allocation in trees: a search for principles[END_REF]. Allometric relationships [START_REF] Niklas | Plant Allometry The scaling of Form and Process[END_REF][START_REF] West | The fourth dimension of life: Fractal geometry and allometric scaling of organisms[END_REF] can be used to constrain the growth of different plant parts so that allometric equations are always satisfied during plant ontogeny [START_REF] Taubert | A review of grassland models in the biofuel context[END_REF]. Although this approach is conceptually and technically simple, it requires empirical measurements of allometric coefficients in multiple environmental conditions. As a result, this approach cannot predict how allometric relationships are likely to vary depending on the environmental conditions encountered by the plant [START_REF] Bloom | Resource limitation in plants -an economic analogy[END_REF]. A second approach is to represent the capture of basic resources by the plant and their transport across different plant organs [START_REF] Thornley | Grassland Dynamics: An Ecosystem Simulation Model[END_REF]. This approach aims at being mechanistic, but it requires the quantification of mechanistic properties such as the resistance to nutrient flow in the plant, as well as processes of internal regulation. A third type of approach relies on various optimization principles. Studies of this type generally consider that allocation in the plant aims at maximizing some criterion, used as a proxy for plant fitness, such as plant relative growth rate [START_REF] Charles-Edwards | Models for mesophyll cell arrangement in leaves of ryegrass (lolium perenne l.)[END_REF][START_REF] Reynolds | Modeling whole-plant allocation in relation to carbon and nitrogen supply: Coordination versus optimization[END_REF]. They sometimes also make use of game-theoretic or adaptive dynamics methods to take into account the ecological and evolutionary impacts of competition between plants [START_REF] Franklin | Modeling carbon allocation in trees: a search for principles[END_REF][START_REF] Mcnickle | Game theory and plant ecology[END_REF].

The idea that plants may allocate their assimilates among organs so as to balance their root activity of water and nutrient uptake and their shoot activity of photosynthesis dates back at least to the work of [START_REF] Brouwer | Nutritive influences on the distribution of dry matter in the plant[END_REF]. Since then, this hypothesis has been variously called "optimal partitioning", "functional equilibrium" or the "balanced growth hypothesis" [START_REF] Poorter | Biomass allocation to leaves, stems and roots: meta-analyses of interspecific variation and environmental control[END_REF]. This simple idea has proven powerful at qualitatively explaining how environmental conditions and perturbations affect patterns of root-shoot allocation [START_REF] Iwasa | Shoot/root balance of plants: Optimal growth of a system with many vegetative organs[END_REF][START_REF] Bloom | Resource limitation in plants -an economic analogy[END_REF][START_REF] Poorter | Biomass allocation to leaves, stems and roots: meta-analyses of interspecific variation and environmental control[END_REF]. Optimal partitioning is thus used in a large number of plant growth models [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF][START_REF] Franklin | Modeling carbon allocation in trees: a search for principles[END_REF]. Some simple models of plant growth that assume optimal partitioning predict an isometric relationship between shoot and root biomass during the exponential phase of growth in nonlimiting conditions; i.e. root and shoot biomass remain proportional [START_REF]Shoot and root activities during steady-state plant growth[END_REF][START_REF] Robinson | Limits to nutrient inflow rates in roots and root systems[END_REF][START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF]. This prediction has received mixed empirical support. According to the meta-analysis of [START_REF] Poorter | Biomass allocation to leaves, stems and roots: meta-analyses of interspecific variation and environmental control[END_REF], ontogenetic shifts in rootshoot ratios are variable across experiments performed so far. For instance, [START_REF] Mcconnaughay | Biomass allocation in plant: ontogeny or optimality? A test along three resource gradients[END_REF] explored the impact of resource gradients on three annual species and found that the root-shoot ratio decreases during plant development. [START_REF] Müller | The effect of nutrient availability on biomass allocation patterns in 27 species of herbaceous plants[END_REF] studied allocation patterns of 27 herbaceous plant species and also found a decreasing rootshoot ratio for 14 species. In contrast, [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF] studied 22 herbaceous plant species during 35 days and found a preferential allocation to roots during plant ontogeny in general, although deviations from isometry were weak in most cases. Arredondo et al.

(1998) also found an increase in root-shoot ratio during plant ontogeny. The variable rootshoot ratios evidenced in these studies question the validity of the optimal partitioning hypothesis.

However, rather than a fundamental flaw in the assumption of optimal partitioning, the discrepancies between data and model predictions could also be due to the requirement for additional model assumptions beyond optimal partitioning. For instance, Reynolds andThornley (1982), Johnson (1985) and Johnson and Thornley (1987) made the point that optimal partitioning implied that plants should be equally limited by shoot and root activities, which do not need to be constant over time, but rather depend on dynamical environmental conditions and potential disturbances. In this vein, [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF] argued that non-isometric relationships may be explained by a progressive nitrogen limitation of plant growth during their experiments or by a decrease of intrinsic root uptake capacity with their age. These two features are susceptible to cause a transient phase of preferential allocation to roots (Ingestad and Ågren 1991). An alternative explanation for non-isometric trajectories might be that root-shoot partitioning is ontogenetically constrained, especially during the early stages of growth [START_REF] Gedroc | Plasticity in root/shoot partitioning: optimal, ontogenetic, or both ?[END_REF]. If allocation is ontogenetically constrained, the shoot-root ratio is likely to differ from the ratio predicted by optimal partitioning. Importantly, even if developmental constraints cease early during plant ontogeny, they are likely to have persistent effects on plant growth trajectories during a transient phase of root-shoot ratio adjustment by the plant. This will be tested here by considering that the initial shoot-root ratio may differ from the one predicted by optimal partitioning, but that subsequent dynamics is controlled by optimal partitioning equations. Transient dynamics, either due to nutrient limitation or to some initial developmental constraints, could potentially explain the discrepancy between steady-state predictions based on optimal partitioning theory and experimental findings. This study aims at quantitatively testing these potential explanations.

The study is structured in three main parts. First, a simple plant growth model is built, which is based on the optimal partitioning hypothesis and on plant functional traits that can be empirically measured in practise. This model represents the basic processes of photosynthesis, nutrient uptake, and root-shoot carbon and nitrogen allocation. Second, we investigate the ability of such a simple model to reproduce patterns of relative growth rate (RGR) and final root weight ratio (RWR) experimentally measured for 25 species by [START_REF] Reich | Variation in growth rate and ecophysiology among 34 grassland and savannah species under contrasting n supply: a test of functional group differences[END_REF]. This particular study was chosen because most plant traits used in the model were measured during the experiments. The other model parameters, which were not experimentally measured, are estimated so as to maximize the model fit to the two growth indicators (RGR and RWR).

Model goodness of fit is then assessed, as well as the realism of fitted parameter values. This part of the study served to determine if the simple model considered is a realistic approximation of plant growth dynamics. Third, armed with this simple but realistic plant growth model, the experiments of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF] are re-analyzed. In these experiments, root and shoot biomass trajectories of 22 plant species were measured experimentally in varying environmental conditions. The simple model is shown not to be able to explain the observed non-isometric root-shoot biomass relationships when model parameters are constrained so that the model accurately fits overall plant growth data. We then test whether adding nitrogen limitation and a decrease of root uptake capacity with root age may lead to the observed non-isometric relationships. We finally explore whether considering ontogenetic constraints through variations in initial shoot-root ratio may improve the model fit to data. Fitted initial shoot-root ratios were finally compared to observed ones by re-analyzing the data of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF].

MODEL DESCRIPTION

The model simulates the growth of herbaceous species in non-limiting conditions of water supply. In this model, a plant is described by its total biomass B(t) at time t. This biomass is divided into above-(B s ) and below-(B r ) ground biomass, so that:

B (t) = B s (t) + B r (t) (1)
Four growth processes are modelled: i) shoot photosynthesis, ii) nitrogen uptake by roots, iii) nitrogen allocation among roots and shoots, and iv) carbon allocation among roots and shoots.

Leaves and stems are not distinguished in the shoot component for two reasons. First, since both leaves and stems contribute to photosynthesis in herbaceous plants [START_REF] Nilsen | Plant Stems: Physiology and Functional Morphology[END_REF], pooling these two plant components makes sense functionally. Second, distinguishing these two plant components would increase model complexity by adding leaf-and stem-specific activity rates (photosynthesis and respiration) and two additional leaf-stem allocation rules for carbon and nitrogen, while these processes are poorly documented. Certainly, increasing the complexity of the model with both stem and leaves components could be easily achieved for cases in which additional information is available.

Plant Development

A simple difference equation is used with one time step representing one hour. The plant biomass at t+1 is given by:

B (t + 1) = B (t) + ∆B (t) = B (t) + P net (t) (2)
where P net (t) is the net primary production at time t.

The increases of shoot and root biomass between times t and t+1 are described by equations

(3) and (4) respectively:

∆B s (t) = a (t) • P net (t) (3) ∆B r (t) = ( 1 -a (t) ) • P net (t) (4)
where a (t) is the portion of net primary production allocated to shoot. The computations of P net (t) and a (t) are detailed below.

Photosynthesis

Grasses perform photosynthesis both in their leaves and stems [START_REF] Aschan | Non-foliar photosynthesis -a strategy of additional carbon acquisition[END_REF].

Although stem photosynthetic rate may differ from leaf photosynthetic rate, as well as masssurface ratios, it will be assumed here for simplicity that these quantities are equal among stems and leaves and thus that net primary production P net (t) can be modelled by:

P net (t) = C x A N (t) • SLA • B s (t) -R r • B r (5)
where [START_REF] Kikuzawa | Toward synthesis of relationships among leaf longevity, instantaneous photosynthetic rate, lifetime leaf carbon gain, and the gross primary production of forests[END_REF]. Equation ( 5) has been abundantly used in plant growth modelling (see for instance [START_REF] Foley | Net primary productivity in the terrestrial biosphere: the application of a global model[END_REF]). Note that this approach is still valid if less strong assumptions are used, namely that the stem-leaf ratio is constant during plant ontogeny and that the stem photosynthetic rate responds similarly to light conditions and plant nitrogen status as leaf photosynthetic rate. In this case, A N (t) should be understood as an effective shoot net photosynthetic rate.

Net photosynthetic rate A N (t) has been shown to be linearly related to the nitrogen content of shoot [START_REF] Lambers | Plant physiological ecology[END_REF]. Following [START_REF] Konings | Growth characteristics, nutrient allocation and photosynthesis of carex species from floating fens[END_REF], the following relationship is used:

A N (t) = ( A max (t) / LNC max ) • ( N s (t) / B s (t) ) (6)
where LNC max is the leaf nitrogen that maximises photosynthesis, N s (t) is the nitrogen content of shoot and A max (t) is the net maximal leaf photosynthetic rate in given light conditions:

A max (t) = A sat • f (I r (t) ) -R s (7)
where A sat is the light saturated gross photosynthetic rate, R s is the shoot respiration rate, I r (t) is the incoming irradiance at time t and f (I r ) is a function varying between zero when I r is null and one when it is optimal. This function accounts for the impact of ambient light on photosynthesis. In the following the irradiance is assumed to be constant during the photoperiod and therefore A max (t) is also constant. Overnight the irradiance I r is null, so f is null and A max (t) equals -R d . As we are interested in the first stages of plant growth the decrease of net photosynthetic rate caused by self-shading and the resulting variations in rootshoot scaling is likely to be negligible and so self-shading will be ignored in the following.

Nitrogen uptake

When nutrient supply is non-limiting, nitrogen uptake is only limited by plant physiology and root biomass. So at time t, a plant is able to absorb at most:

N p = U max • B r (t) (8)
where U max is the mass-based root effective uptake capacity.

Roots are assumed to be able to adjust nitrogen uptake so as to match the nitrogen demand of the plant N d [START_REF] Schippers | Competition for light and nitrogen among grassland species: a simulation analysis[END_REF]. The latter corresponds to the amount of nitrogen required for the leaf content of new leaf biomass to be equal to LNC max . N d is thus given by:

N d = ( LNC max • ∆B s (t) ) / a N (t) (9)
where ∆B s (t) is the shoot biomass produced between t and t+1 and a N (t) is the fraction of nitrogen captured between t and t+1 which is allocated to shoot. Thereafter, assimilated nitrogen N u equals the minimum of N p and N d :

N u (t) = min ( N p (t) , N d (t) ) (10)

Nitrogen allocation

Following [START_REF] Dybzinski | Evolutionarily stable strategy carbon allocation to foliage, wood, and fine roots in trees competing for light and nitrogen: an analytically tractable, individual-based model and quantitative comparisons to data[END_REF], it is assumed that a fixed fraction of assimilated nitrogen is allocated to the shoot: a N (t) = a N . An alternative way for modelling nitrogen partitioning would be to use optimal partitioning theory [START_REF] Mäkelä | Optimal co-allocation of carbon and nitrogen in a forest stand at steady state[END_REF]. But to apply this theory, it would be necessary to know the relationship between root uptake efficiency and root nitrogen content or to make some assumptions on the relationship between shoot and root nitrogen content [START_REF] Mäkelä | Optimal co-allocation of carbon and nitrogen in a forest stand at steady state[END_REF][START_REF] Valentine | Modeling forest stand dynamics from optimal balances of carbon and nitrogen[END_REF].

Carbon allocation

An optimal allocation model is used for carbon allocation [START_REF] Dewar | Optimal function explains forest responses to global change[END_REF], in which plants are assumed to allocate assimilates so as to maximize their relative growth rate.

Assuming that biomass and leaf nitrogen content at time t are known, we look for an allocation to shoot a (t) such that RGR (t+1) is maximal. RGR (t+1) is given by:

RGR (t+1) = ∆B (t+1) / B (t+1) = P net (t+1) / B (t+1) (11)
From equations ( 5) and ( 6), RGR (t+1) maximization is equivalent to maximize:

C x (A max (t+1) / LNC max ) • ( N s (t+1) / B s (t+1) ) • SLA • B s (t+1) -R r • B r (12)
Given that:

N s (t+1) = N s (t) + N u (t) (13) RGR (t+1) is maximal when N u (t) is maximal.
In other words RGR (t+1) is maximal when the nitrogen demand at t+1 is equal to the potential uptake N p (t):

LNC max • a (t) • P net (t) = N p (t) (14)
Hence the shoot allocation factor a(t) is given by:

a (t) = N p (t) / (LNC max • P net (t) ) (15)
In summary, the model takes as input seven parameters (A max , R s , R r , LNC max , SLA, U max and a N ). Most of the parameters are commonly measured plant functional traits [START_REF] Kattge | Try -a global database of plant traits[END_REF]. A max , R s , R r , LNC max and SLA were measured in the experiments of [START_REF] Reich | Variation in growth rate and ecophysiology among 34 grassland and savannah species under contrasting n supply: a test of functional group differences[END_REF],

and LNC max and SLA in the experiments of [START_REF] Meziane | Interacting components of interspecific relative growth rate: constancy and change under differing conditions of light and nutrient supply[END_REF]. Note that plant senescence was neglected, since we are interested here in the first stages of plant development.

PARAMETER ESTIMATION AND MODEL SELECTION

In this section, the data set used to test the model ability to reproduce real plant growth dynamics is first presented. The fitting procedure of the remaining unmeasured model parameters is then detailed. The results of this model-data comparison procedure are presented in the last subsection.

Plant growth data

An experiment performed by [START_REF] Reich | Variation in growth rate and ecophysiology among 34 grassland and savannah species under contrasting n supply: a test of functional group differences[END_REF] is used, in which 34 herbaceous and woody plant species were grown in monoculture under controlled environmental glasshouse conditions, for nine weeks after germination. Herein we use 25 of these species (just the herbaceous nonleguminous ones), growing under fertilized conditions (ignoring ambient grown plant) to minimize N limitations. Beginning two weeks after sowing, pots received 30 ml of half strength Hoagland's solution three times per week. Pots were watered as needed between treatment applications to maintain soils near field capacity. Supplemental lighting provided an additional 130 -170 µmol.m 2 .s -1 above ambient light levels during a 14 h.

photoperiod. Each three weeks, plants were harvested and the biomass of the different plant components (leaves, stem and roots) were measured. From these measurements, several quantities were computed: the root weight ratio RWR equal to the root biomass divided by the total plant biomass; and the relative growth rate RGR computed as in [START_REF] Evans | The quantitative analysis of plant growth[END_REF]:

RGR = [ ln(B (t 2 )) -ln(B (t 1 )) ] / ( t 2 -t 1 ) (16)
where t 1 and t 2 are harvesting dates.

Several plant functional traits were also measured: the light saturated photosynthetic rate A sat , the shoot and root respiration rate R s and R r , the leaf nitrogen content LNC, the root nitrogen content RNC and the specific leaf area SLA. RNC is not a model parameter, but it enables to compute the nitrogen allocation coefficient a N with the following equation:

a N = ( N s (t end ) / B s (t end ) ) / [ ( N s (t end ) / B s (t end )) + RNC • (B r (t end ) / B s (t end )) ] (17)
This data set enables us to assess the model ability to explain observed plant growth patterns, when it is strongly constrained by empirically measured plant traits.

Fitting the model on experimental data

Following [START_REF] Goudriaan | Modelling potential crop growth processes[END_REF], the initial shoot-root ratio was set to one. Initial root and shoot biomass values do not affect the computed growth indicator, so they are arbitrarily set to 0.5mg. The maximal nutrient uptake efficiency U max was not measured in the experiments of [START_REF] Reich | Variation in growth rate and ecophysiology among 34 grassland and savannah species under contrasting n supply: a test of functional group differences[END_REF] and thus had to be estimated. Besides, the maximal photosynthetic efficiency A max was measured for an irradiance of 1000 µmol.m 2 .s -1 , while plants were not grown under constant light conditions. The effective A max during the experiments, resulting from the variable light conditions, was therefore estimated. These two parameters were estimated by fitting the plant growth model to the growth data of the experiments of [START_REF] Reich | Variation in growth rate and ecophysiology among 34 grassland and savannah species under contrasting n supply: a test of functional group differences[END_REF]. A distance ε between model predictions and experimental data was defined. It was based on two experimentally measured growth indicators: the plant relative growth rate RGR, and the final root weight ratio RWR. ε was defined as:

2 2 1 2 ) ( ) ( 2 1 obs i i sim i obs i X X X ∑ = - × = ε (18)
where X sim = ( RGR sim , RWR sim ) and X obs = ( RGR obs , RWR obs ) are the growth indicators of the simulations and of the experimental data respectively. This distance ε was computed on a 50x190 grid of parameter values described in Table 1 and the parameter set that minimized this distance was retained. The parameter space was chosen so as to include the reported parameter values found in the literature (Table 1).

Note that the interval chosen for A max does not include the largest values measured by [START_REF] Reich | Variation in growth rate and ecophysiology among 34 grassland and savannah species under contrasting n supply: a test of functional group differences[END_REF]. In this experiment, the photosynthetic efficiency was measured in optimal light conditions (f (I r ) = 1). But plants were not grown in optimal light conditions, so that the effective photosynthetic efficiency is necessarily smaller than light saturated photosynthetic efficiency.

Empirical growth indicators have been measured with observation errors of 5 to 10% [START_REF] Reich | Variation in growth rate and ecophysiology among 34 grassland and savannah species under contrasting n supply: a test of functional group differences[END_REF]. Therefore all sets of parameters leading to ε values smaller than 0.05 are retained.

They form an interval of likely values for A max and U max , which were relatively narrow (Table 2).

Results

The coefficients of variation of the two model parameters A max and U max were smaller than 5 and 10% respectively (Table 2), which means that the growth data used were sufficiently informative to obtain accurate parameter estimates of the minimal version of the model.

Model parameter estimates are realistic compared to the range of values reported in the literature (Table 1). Importantly, the remaining lack-of-fit of the model to data, leading to a residual error ε min (Table 2), was negligible compared to the observed interspecific variations of growth indicators (Figure 1, average NRMSE=3.47%). This means that this simple traitbased model was sufficient to capture interspecific differences in growth rates as well as in root-shoot carbon and nitrogen partitioning.

The growth dynamics predicted by the fitted minimal version of the model consists of a short transient phase during which carbon allocation a (t) varies, followed by a steady-state regime of exponential growth during which carbon allocation is constant (Figure 2A) and root and shoot growth rates scale proportionally (Figure 2B), leading to an allometric coefficient β equal to one (Figure 2C). The predicted duration of the transient phase depends on the species identity, but never exceeds 15 days (Figure 2C). In the 35-day experiments of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF], the first measurements occurred 15 days after germination (shown by the dashed lines in Figure 2B-C). So the duration of transient dynamics in the minimal version of the model is insufficient to explain the observed deviations from isometry observed in the 35day experiments of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF]. These authors suggested that the observed deviation from isometry could come from a progressive appearance of nitrogen limitation in the experimental setting in which nitrogen was added in fixed amounts. They also discussed that a decrease of intrinsic root uptake capacity with root age could contribute to nonisometric growth trajectories. These two additional processes were included in the model to assess their ability to explain observed patterns. It has been further assessed whether a modification of the initial (ontogenetically constrained) shoot-root ratio could significantly contribute to observed non-isometric trajectories.

ALLOMETRIC PREDICTIONS

In this section, the data set used to test the model predictions on root-shoot allometry is first presented. Two additional model ingredients are then introduced: (i) the consideration of nitrogen consumption during the experiment potentially causing some nitrogen limitation for plants, especially at the end of the experiment; and (ii) the inclusion of a decrease in root uptake capacity as they age. The impact of initial shoot-root ratio on allometric patterns was also investigated. Third, the model-data fitting procedure was detailed, as well as the associated test of whether the different models studied are able to reproduce the empirical root-shoot allometries. Fourth, the results of this model-data comparison procedure are presented.

Plant growth data

Given that [START_REF] Reich | Variation in growth rate and ecophysiology among 34 grassland and savannah species under contrasting n supply: a test of functional group differences[END_REF] did not perform detailed measures of allometric relationships (they only performed three sequential harvests), a second data set collected by [START_REF] Meziane | Interacting components of interspecific relative growth rate: constancy and change under differing conditions of light and nutrient supply[END_REF] and [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF] is used. In this experiment, a total of 1150 plants from 22 different herbaceous plant species were grown in hydroponic sand monoculture in factorial combinations of high (1100 µmol.m 2 .s -1 PAR) and low (200 µmol.m 2 .s -1 PAR) irradiance crossed with a full strength and a 1/6 dilution of Hoagland's hydroponic solution. Each plant grew in a separate 1.3 dm 3 container in a growth chamber with 15/9 h light:dark cycles. Each plant grew in one of four resource environments: high (L, 1100 µmol.m 2 .s -1 ) and low (l, 200 µmol.m 2 .s -1 ) irradiance combined with high (N, fullstrength Hoagland's nutrient solution) and low (n, 1/6 dilution) external nutrient concentrations. These four experimental treatments are termed LN, Ln, lN and ln treatments respectively. Each container was filled to field capacity with the nutrient solution three times a day. Plants were harvested and biomass of leaves, stems and roots measured at 15, 20, 25, 30 and 35 days post-germination. Two plant functional traits were measured: LNC and SLA.

RWR was computed, as well as average RGR, computed as the slope of a regression of the natural logarithm of plant dry mass on harvest date. Allometric relationships between shoot and root biomass along plant ontogeny was further measured, using an equation of the form:

B s = α • B r β (19)
which can be re-written as:

ln( B s ) = ln( α ) + β • ln( B r ) (20)
The allometric coefficient β was thus computed as the slope of a regression of the natural logarithm of shoot dry mass on the natural logarithm of root dry mass. This second data set

was not used for model checking, since fewer plant functional traits were measured empirically, and thus it would have been a less conservative test of the model with a larger number of unmeasured model parameters and a smaller number of growth indicators to match. Rather, it was used to confront the model predictions with empirically-measured allometric data.

Adding nitrogen consumption by plants in the containers

Following [START_REF] Engels | Assessing the ability of roots for nutrients acquisition[END_REF], Hane's relationship was used to model the dependency of uptake rate U(t) to substrate concentration:

U(t) = U max • [N] soil (t) / ( K m + [N] soil (t) ) (21)
where K m is the substrate affinity and [N] soil (t) is the nitrogen concentration in soil at time t.

Initially soil nitrogen content equals:

N soil ( t 0 ) = [N] Hoagland • V • C soil (22)
where [N] Hoagland is the nitrogen concentration in the hydroponic solution used in the experiment (0.210 g.L -1 for full strength solution, 1/6 of this value in the low nitrogen treatment), V is the container volume (1.3 L) and C soil is the volumetric soil moisture content remaining at field capacity (about 5% according to Tucker 1999). Then soil nitrogen content is computed at each time step as:

N soil (t) = N soil (t -1) -N u (t -1) (23)
where N u (t) is the amount of nitrogen absorbed by the plant at time t (eq. 12). Every 8 hours each container is filled to field capacity with the nutrient solution, so this dynamics of nitrogen concentration decrease in the container is restarted.

Adding a decrease of root uptake capacity with root age

The model of decrease in root uptake efficiency as they age is based on the observations of [START_REF] Volder | Rapid decline in nitrate uptake and respiration with age in fine lateral roots of grape: implications for root efficiency and competitive effectiveness[END_REF]. Root biomass is divided in several layers. Each layer has its own biomass, age and nitrogen uptake capacity. At the beginning of each time step t, a root layer is added with a biomass corresponding to the newly produced root biomass ∆B r (t -1). A root layer i will have a varying with time nitrogen uptake capacity U i (t) given by:

U max,i (t) = U max • (1 + 2 • e -ρ (t -ti) ) / 3 ( 24
)
where ρ is the decay rate of root nitrogen uptake efficiency, and t i is the time of appearance of the root layer i. Following [START_REF] Volder | Rapid decline in nitrate uptake and respiration with age in fine lateral roots of grape: implications for root efficiency and competitive effectiveness[END_REF], it is assumed that after some days root uptake efficiency stabilizes around one third of maximal efficiency.

Adding variation in initial root-shoot ratio

Since the first measurements in the experiment of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF] occurred 15 days after germination, no information is available on growth trajectories during the very first days of the experiment. Two hypotheses were compared regarding allocation patterns during these first 15 days. (H 0 ) Biomass is optimally allocated, so that the shoot root ratio at first measurement, called initial shoot-root ratio R 0 in the following, is the ratio required for optimal partitioning, R opt . (H 1 ) Because of ontogenetic developmental constraints, the initial shoot-root ratio R ont differs from the optimal one:

(H 0 ) : B s (t 0 ) = R opt • B r (t 0 ) (H 1 ) : B s (t 0 ) = R ont • B r (t 0 ) (25)
These additional model ingredients are used to build three models. The first one does not include nitrogen limitation and assumes that initial shoot-root ratio corresponds to the ratio predicted by optimal partitioning (M 0 ). The second one additionally includes nitrogen consumption and decrease in root uptake capacity with age (M 1 ). The third one additionally authorizes initial shoot-root ratios differing from the optimal one (M 2 ). Each model was fitted to the same data of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF] which have been obtained in four resource environments: high and low irradiance combined with high and low external nitrogen concentration. Importantly, all modelled dynamics are based on optimal partitioning theory: in model M 2 , developmental constraints cease after the first 15 days and only affect the initial conditions.

Fitting models to experimental data

Five parameters, the photosynthetic efficiency A max , the shoot and root respiration rates R s and R r , the nitrogen allocation a N and the maximal nutrient uptake efficiency U max were not measured in the experiments of [START_REF] Meziane | Interacting components of interspecific relative growth rate: constancy and change under differing conditions of light and nutrient supply[END_REF]. Average shoot and root respirations R s mean and R r mean were computed from the data of [START_REF] Reich | Variation in growth rate and ecophysiology among 34 grassland and savannah species under contrasting n supply: a test of functional group differences[END_REF], and the average values R s mean = R r mean = 40 nmol.g -1 .s -1 was used for all the species studied in [START_REF] Meziane | Interacting components of interspecific relative growth rate: constancy and change under differing conditions of light and nutrient supply[END_REF]. It thus remains three parameters to estimate in the minimal version of the model: A max , a N and U max . Up to three additional parameters, ρ, K m and R 0 , need to be estimated in the additional versions of the model described above.

A 25x26x5x7x9x19 grid of parameter values was used (Table 3). The parameter space has been chosen so as to include the reported parameter values found in the literature (Table 3). In M 0 , parameters ρ and K m were set to zero. In models M 0 and M 1 , initial shoot-root ratio R 0 was equal to the ratio predicted by optimal partitioning.

The same distance ε between model predictions and experimental data is used. This time, between four and six model parameters have to be estimated with only two growth indicators.

Consequently, a large array of parameter sets can lead to model predictions matching the two growth indicators. Therefore, rather than trying to estimate the model parameters, the model simulations are filtered so that they fit the available growth indicators (Jabot and Bascompte 2012). As previously, model parameter sets leading to ε values smaller than 0.05 are retained.

These retained realistic simulations are then used to explore the range of allometric relationships between root and shoot biomass that the three models are able to predict. To this end, for each model simulation, the same procedure as in [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF] was used to compute the allometric coefficient β: plant biomass was simulated at 15, 20, 25, 30 and 35 days post-germination and β was computed as the slope of a regression of the natural logarithm of shoot dry mass on root dry mass. To quantify the predictive ability of the models, the relative distance between simulated and empirical values of allometric coefficients β is computed for every retained parameters sets and the sets which lead to the smallest distance, called d min are kept:

                - = 2 2 min ) ( ) ( min obs sim obs d β β β (26)
where β obs and β sim are empirical and simulated values of the allometric coefficients.

Results

Model M 0 was found (as in the first section) to produce β values equal to one (Figure 3A).

Adding a decrease in root uptake efficiency with age and nitrogen consumption (model M 1 ) caused plants to become progressively limited in nitrogen and to allocate an increasing amount of biomass to roots as they grow. This led to β values smaller than one (Figure 3B).

When initial shoot-root ratio was large, plants also allocated more biomass to roots until an optimal ratio was reached (model M 2 ). This also led to β values smaller than one (Figure 3C), although it affected the beginning of the growth dynamics rather than the end as observed with nutrient limitation (Figure 3B). These two processes were found to reinforce each other, since they acted at different growth stages (Figure 3D).

The three models were fitted to the data of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF]. Some models failed at reproducing some of the growth indicators with an average relative error ε smaller than 0.05. Model M 0 did not succeed in reproducing RGR and RWR of 8, 14, 18 and 21 out of the 22 species in the LN, lN, Ln and ln treatments respectively (Figure 4). When the effects of nitrogen limitation were included in the model (model M 1 ), the number of accurately reproduced growth patterns strongly increased: the growth indicators of 18, 20, 20 and 22 out of the 22 species could be reproduced in the lN, LN, Ln and ln treatments respectively. The full model (M 2 ) was able to accurately reproduce the growth indicators of all species, except in the LN treatment (20 of the 22 species).

When models succeeded in reproducing the plant growth indictors, the allometric coefficients β that they predicted were compared to observed ones. Model M 0 failed at reproducing such allometric coefficients, irrespective of the environmental treatments. Adding nitrogen limitation weakly increased the predictive ability of the model: the relative distance d min between empirical and simulated allometric coefficients was equal to zero for only seven species across all treatments, i.e in 8% of the cases. In contrast, model M 2 succeeded in reproducing observed allometric coefficients with d min values almost always equal to zero, except for Rumex acetosa (β = 0.748) in the LN treatment, Deschampsia cespitosa (β = 1.261) in the Ln treatment and Silene cucubalus (β = 0.804), Polygonum lapathifolium (β = 0.906), Plantago major (β = 0.948) and Panicum capillare (β = 1.000) in the lN treatment (Figure 4).

Even if a model fails at accurately predicting the empirical allometric coefficients, it may still make close predictions. β values simulated with model M 1 were closer to empirical values of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF] than those simulated with model M 0 for all species with β values smaller than one (Figure 4). The relative distance d min between empirical and simulated β values was larger than 0.10 in 83% of the cases for model M 0 , while for model M 1 , d min was smaller than 0.025 in 26% of the cases and smaller than 0.10 in 50% of the cases. The full model (M 2 ) had relative errors d min on β values smaller than 0.10 for all species in all treatments, except for Rumex acetosa in the LN treatment.

Model validation

The modelling results make an additional prediction: in cases in which allometric coefficients β are significantly different from one, this should be due to the difference between the initially constrained shoot-root ratio and the ratio required for optimal partitioning. If β is smaller than one, the initial shoot-root ratio should be larger than the optimal ratio and vice versa ([Supplementary information], Figure S2) This final prediction was tested by re-analyzing the data of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF]. The average initial shoot-root ratio and the standard deviation were computed from available biomass data, and a weighted least square regression was used to assess the accuracy of model predictions (Figure 5).

For the LN and ln treatments, the initial shoot-root ratio predicted by the full model (M2) fits the biomass raw data with a good accuracy (r = 0.87 for the LN treatment and r = 0.53 for the ln treatment). For the Ln and lN treatments, model predictions are not so good (r = 0.12 for the lN treatment and r = 0.33 for the Ln treatment). In these treatments only 5 (respectively 6) species out of 21 had a β value significantly different from 1 [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF].

When the regression was performed solely on species with an allometric coefficient significantly smaller than one, the regression coefficient sharply increased to 0.80 for the lN treatment and to 0.63 for the Ln treatment. This means that model predictions regarding initial root-shoot ratios were close to observations in cases in which this ontogenetic constraint was necessary to explain root-shoot trajectories. In cases in which a simpler model was sufficient to account for observations, the uncertainty of the initial root-shoot ratio estimation obscured the predictions (Figure 5).

DISCUSSION

This study aimed at testing whether observed root-shoot allometries during plant development could be explained by the hypothesis of optimal partitioning. The approach was thus twofold.

First, a simple model of plant growth was built based on commonly measured plant functional traits, and it has been tested whether this model was a sufficiently detailed account of plant growth to reproduce various growth indicators of experimental studies. The model succeeded in reproducing these indicators with a very good accuracy (Figure 1). Two plant traits were not empirically measured during the experiments and were estimated to reproduce the two plant growth indicators. Fitted trait values were within the ranges reported in the literature, and were thus realistic. This first study part was essential to discard the possibility that the discrepancy between theoretical predictions and empirical allometries would be due to a poor modelled representation of plant growth.

The second part of the study aimed at using this simple model to predict root-shoot allometry during plant development. The results of previous simpler models were recovered: the optimal partitioning hypothesis led to an isometric growth of roots and shoot [START_REF] Charles-Edwards | Models for mesophyll cell arrangement in leaves of ryegrass (lolium perenne l.)[END_REF], in contrast with the empirical findings of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF]. This steady-state isometry was further found to be preceded by a short transient period of non-isometric growth, during which plants were dynamically adjusting their allocation coefficient if initial root-shoot ratio were not optimal (Figure 2). It was then tested whether adding complementary model ingredients could lengthen the duration of this transient phase and change the shape of the root-shoot allometry. Root senescence was added in the model, as well as the nitrogen consumption by plants in the experimental containers. These first two ingredients improved the ability of the models to reproduce growth indicators when light and nitrogen were limiting. However they were insufficient to explain the results of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF] on allometric trajectories. The initial shoot-root ratio was then varied to represent initial developmental constraints. With this third ingredient, most empirical findings of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF] could be reproduced (Figure 4).

These simulation results show that to be a reasonable approximation of plant allocation scheme, the optimal partitioning framework needs to be complemented by a number of complementary processes which lead to transient phases of allocation adjustment by the plant.

A combination of these processes was found to be necessary to recover the empirical findings of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF]. More precisely, initial shoot-root ratio had to be different from the ratio predicted by optimal partitioning to recover most allometric coefficients which were significantly different from one. [START_REF] Gedroc | Plasticity in root/shoot partitioning: optimal, ontogenetic, or both ?[END_REF] provided evidence that ontogenetic constraints were likely to play a role in allometric trajectories by statistically analyzing plant growth trajectories. The modelling approach proposed here enabled us to mechanistically incorporate various processes that have been suggested in the literature to cause shifts in rootshoot scaling during plant ontogeny, and to test their respective influences quantitatively.

Furthermore, this approach demonstrated that although adding some initial developmental constraints may be needed to recover experimental findings, such developmental constraints were no longer needed during the subsequent plant growth dynamics phase, which was controlled by optimal partitioning mechanisms.

The hypothesis that initial shoot-root ratios may differ from the ratio required for optimal partitioning is supported by the re-analysis of the experiments of Shipley and Meziane (2002) (Figure6). Moreover the finding that initial shoot-root ratio should be larger than 1 during the very first days of plant growth after germination is consistent with observations (Jurado and Westoby 1992;[START_REF] Leishman | The role of large seed size in shaded conditions: experimental evidence[END_REF]. The hypothesis that this initial shoot-root ratio should be at least partially developmentally constrained is also consistent with observations [START_REF] Evans | Comparative root morphology of some pasture grasses and clovers[END_REF][START_REF] Kitajima | Do shade-tolerant tropical tree seedlings depend longer on seed reserves? Functional growth analysis of three Bignoniaceae species[END_REF] and the biological fact that initial plant growth is ensured by the consumption of seed reserves. The explanation modelled here for non-isometric root-shoot trajectories is hence uniquely based on a transient dynamics controlled by optimal partitioning equations. The model proved powerful in the non-constant nutrient conditions of [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF]'s experiments (when looked at an hourly resolution). This model could hence similarly be used to quantitatively predict root-shoot dynamics in response to disturbances [START_REF] Mäkelä | Acclimatation in dynamic models based on structural relationships[END_REF].

The model proposed here does not take into account self-shading, shoot senescence nor water use by the plant. It is thus not able in its present form to assess potential explanations for observed ontogenetic shifts in root-shoot scaling evidenced in longer term experiments (e.g., Mueller et al. 2000) for which these three effects could play an additional role. Similarly, it cannot be used in its present form to assess the potential role of temperature increase or water shortage in biomass allocation. Since these last two environmental pressures are likely to be important according to current climate change scenarios, further model refinements to add temperature and water effects constitute very interesting perspectives. Such future model developments will be eased by the general approach of progressive model building through quantitative assessment that has been developed in this study.

. In this study, optimal allocation equations were based on the idea that relative growth rate should be a good proxy for plant fitness and thus should be optimized by evolution. However, plants have evolved in competitive environments, so that for a plant, maximizing its growth in isolation is not necessarily the best strategy to maximize its growth in competitive conditions [START_REF] Mcnickle | Game theory and plant ecology[END_REF]. Future work addressing this issue would be valuable. They would require the consideration of below-ground competition for soil resources and of aboveground competition for light [START_REF] Tilman | Plant strategies and the dynamics and structure of plant communities[END_REF][START_REF] Schieving | Carbon gain in a multispecies canopy: the role of specific leaf area and photosynthetic nitrogen-use efficiency in the tragedy of the commons[END_REF][START_REF] Gersani | Tragedy of the commons as a result of root competition[END_REF][START_REF] O'brien | Games roots play: effects of soil volume and nutrients[END_REF]. Such work would enable one to understand whether the use of new plant fitness proxies in optimal partitioning modelling could also produce the longlasting transient dynamics that have been here evidenced. S1.

AKNOWLEDGEMENTS

Théophile Lohier was funded by the Institut national de Recherche en Sciences et Technologies pour l'Environnement et l'Agriculture (IRSTEA) and the region Auvergne. [START_REF] Reich | Variation in growth rate and ecophysiology among 34 grassland and savannah species under contrasting n supply: a test of functional group differences[END_REF].

TABLES TABLE 1. Parameter ranges reported in the literature and those used for model calibration with the experiments of

A max nmol.g -1 .s -1 U max mg.g -1 .h -1
Range reported in the literature 150 -900 1 0.22 -0.92 2

Range used for model calibration 100 -600 0.10 -2.00

Grid step used for model calibration 10 0.01

1 [START_REF] Reich | Variation in growth rate and ecophysiology among 34 grassland and savannah species under contrasting n supply: a test of functional group differences[END_REF] ; 2 [START_REF] Maire | Trade-off between root nitrogen acquisition and shoot nitrogen utilization across 13 co-occurring pasture grass species[END_REF] [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF]. and from data with β significantly different from one (r s ).

  FIG. S2. Shoot-root ratio required for optimal partitioning (+) and shoot-root ratio minimizing the distance between empirical and simulated allometric coefficients (○) in (A) the LN treatment,. (B) the Ln treatment., (C) the lN treatment and. (D) the ln treatment. Species are sorted by increasing β values. The abbreviations for species names used in the x-axis are detailed in TableS1.

  FIGURES

  A N (t) is the leaf net photosynthetic rate expressed per unit leaf area, SLA is the specific leaf area, B s (t) is the shoot biomass, R r is the root respiration rate, B r is the root biomass and

	C is a constant accounting for the conversion of assimilated CO 2 into dry matter content. C is
	calculated from the stochiometry of photosynthetic reactions: to synthesize one mole of
	glucose (C 6 H 12 O 6 ) weighting 180g, six moles of carbon dioxide (CO 2 ) are needed, hence C =
	180/6 = 30

TABLE 2 .

 2 Average estimates of model parameters and associated minimal distance between model predictions and empirical data. Standard deviations of parameter estimates which leadto an ε smaller than 5% are indicated inside the parentheses.

	Species	A max	U max	ε min
		nmol.g -1 .s -1	mg.g -1 .h -1	%
	Achillea millefolium (AcM)	262 (9.85)	0.87 (0.064)	0.6

TABLE 3 .

 3 Parameter ranges reported in the literature and those used for model calibration with the experiments of

(Reich et al.,

2003) ; 2[START_REF] Maire | Trade-off between root nitrogen acquisition and shoot nitrogen utilization across 13 co-occurring pasture grass species[END_REF] ; 3[START_REF] Volder | Rapid decline in nitrate uptake and respiration with age in fine lateral roots of grape: implications for root efficiency and competitive effectiveness[END_REF] ; 4 (Morris 1980) ; 5[START_REF] Goudriaan | Modelling potential crop growth processes[END_REF].

SUPPLEMENTARY INFORMATION

TABLE S1. List of abbreviations used for the species studied by [START_REF] Shipley | The balanced-growth hypothesis and the allometry of leaf and root biomass allocation[END_REF].

Abbreviation Species

Acc


OPERATIONAL FRACTIONATION OF ARSENIC IN A MINING SOIL ABSTRACT

Keywords: Arsenic, Mining Soil, Bearing Phases, Fractionation

The environmental assessment of historical mining sites contaminated with various inorganic species, as also the determination of remediation targets require a better knowledge of pollutant-bearing phases. The objective of this paper was to investigate the operational fractionation of arsenic (As) in an As contaminated (ca. 3 wt%) soil collected from a gold mining site in France, where mining activities and smelting processes of gold ores took place until 2004.

Single and sequential chemical extraction procedures were performed to evaluate the potential mobility of As according to the fractions to which it was associated in the contaminated soil, and also the part of As sorbed onto soil particles.

These procedures used provided information on the nature of arsenic fixation processes in the mining soil studied. Results showed that As appeared mainly sorbed (ca. 72%) with iron on the soil particles, in particular amorphous iron oxyhydroxides. These results confirmed a potential risk of arsenic mobilization under acidic and/or reducing conditions, which frequently occur in mining environment.

This research evidenced the importance of the knowledge of As-bearing phases to better predict the mobility of this pollutant and in turn improve the management of arsenic contaminated sites.

INTRODUCTION

Arsenic (As) is a toxic metalloid naturally present in the environment from various sources such asvolcanism and weathering of the bed rock. Arsenic content in natural soils typically ranges from 0.1 to 50 mg.kg -1 [START_REF] Backer | Chemical monitoring of soils for environment quality and animal and human health[END_REF]. Anthropogenic activities also contribute to the geochemical cycling of As in a variety of ways [START_REF] Tamaki | Environmental Biochemistry of Arsenic[END_REF][START_REF] Yan-Chu | Arsenic distribution in soils[END_REF][START_REF] Smedley | A review of the source, behavior and distribution of arsenic in natural waters[END_REF]. As is used in agriculture as a component of pesticides and also used in wood preservation and glassmaking. It is also dispersed into the environment through the burning of fossil fuels. Additionally, due to its geochemistry, As is used as an indicator element in geochemical prospecting for various types of mineral ores (gold, silver, copper, uranium, etc.) [START_REF] Boyle | The geochemistry of arsenic and its use as an indicator element in geochemical prospecting[END_REF]. Soils near mining sites receive significant inputs of As and heavy metals during mining operations and minerals processing [START_REF] Azcue | Impact of abandoned mine tailings on the arsenic concentrations in Moira lake, Ontario[END_REF][START_REF] Hudson-Edwards | Mineralogy and geochemistry of alluvium contaminated by metal mining in the Rio Tinto area, southwest Spain[END_REF]Savage et al., 200). As concentrations in tailing piles and tailing-contaminated soils can reach up to several thousand mg.kg -1 [START_REF] Smedley | A review of the source, behavior and distribution of arsenic in natural waters[END_REF]. However, unlike organic contaminants, As in soils cannot be decomposed chemically or biologically.

Environmental assessment of contaminated sites requires the determination of pollutants leaching potentials. Although a wide variety of leaching tests are available in the literature to determine the potential mobility of pollutants [START_REF] Van Der Sloot | Harmonization of leaching/extraction tests[END_REF][START_REF] Kosson | An Integrated Framework for Evaluating Leaching in Waste Management and Utilization of Secondary Materials[END_REF], very few have been designed to provide information on the nature of their fixation processes in studied matrices.

In this context, an approach was developed based on the use of leaching tests (single and sequential chemical extraction procedures) in conjunction with a mineralogical study. Indeed, knowledge of the pollutants operational fractionation in the studied matrix is required to better understand mechanisms regulating the leaching behavior of inorganic contaminants of interest.

The objective of this paper was to investigate the nature of As distribution and association to the soil constituents. The correlation between mineralogy and leachability of As was examined in a contaminated soil collected from a French gold mining site. Mineralogical analyses were performed following leaching tests to characterize the studied matrix.

MATERIAL AND METHODS

Sampling and preparation of soil

A contaminated soil collected from a gold mining site located in France was used for the study. The site was characterized by its geological heterogeneity and the utilization of a large variety of mining processes for ores extraction and gold purification until 2004. It was already known to present high concentrations of As [START_REF] Clozel | Physical, chemical, and biological treatability of contaminated soils (Final report -BRGM/RP-52065-FR)[END_REF][START_REF] Chatain | Arsenic behavior in mining soil[END_REF][START_REF] Chatain | Characterization of arsenic and other inorganic constituents potential mobilization from soils collected from a gold mining site[END_REF]Chatain et al., 2005a and b;[START_REF] Bayard | Mobilisation of Arsenic from a mining soil in batch slurry experiments under bio-oxidative conditions[END_REF].

Prior to characterizations and liquid-solid partitioning experiments, the soil samples was air-dried at room temperature (20 ± 3ºC) for one day, sieved through a stainless steel 2 mm mesh sieve to remove coarse debris and gravel, homogenized, and finally stored at 4°C in darkness.

Soil characterization

Chemical characterization

A contaminated soil collected from a gold mining site located in France was used for the study. The site was characterized by its geological heterogeneity and the utilization of a large variety of mining processes for ores extraction and gold purification until 2004. It was already known to present high concentrations of As [START_REF] Clozel | Physical, chemical, and biological treatability of contaminated soils (Final report -BRGM/RP-52065-FR)[END_REF][START_REF] Chatain | Arsenic behavior in mining soil[END_REF][START_REF] Chatain | Characterization of arsenic and other inorganic constituents potential mobilization from soils collected from a gold mining site[END_REF]Chatain et al., 2005a and b;[START_REF] Bayard | Mobilisation of Arsenic from a mining soil in batch slurry experiments under bio-oxidative conditions[END_REF].

Prior to characterizations and liquid-solid partitioning experiments, the soil samples was air-dried at room temperature (20 ± 3ºC) for one day, sieved through a stainless steel 2 mm mesh sieve to remove coarse debris and gravel, homogenized, and finally stored at 4°C in darkness.

Mineralogical characterization

X-ray diffraction

Mineralogical characterization was carried on three replicates of the less than 2 mm soil fraction grounded to powder (<50µm)out using a SIEMENS® D500 advance x-ray diffraction (XRD) instrument equipped with a copper anticathode. The samples were scanned from 3° to 70° 2 at a scan rate of 0.02° 2/second. The results obtained were processed using the DIFFRAC Plus EVA © software (BRUKER AXS ® ) and the ICDD ® database (International Centre for Diffraction Data).

Scanning electron microscopy equipped for energy dispersive spectroscopy

Mineralogical investigation was also performed on two replicates of the less than 2 mm soil fraction by scanning electron microscope (SEM) (JEOL ® 840A LGS) using the backscattered electrons (BSE) mode, a voltage of 20 kV, coupled with an X-ray energy dispersive spectrometer (EDS). Samples for SEM-EDS analyses were prepared by attaching the soil particles on an adhesive tape placed on a stub (sample holder of the SEM). No coating of the samples with a conductive material was required prior to the analysis.

Leaching tests

Sequential chemical extractions

A sequential chemical extraction procedure adapted from [START_REF] Tessier | Sequential extraction procedure for the speciation of particulated metals[END_REF] was carried out in triplicates to determine the major arsenic binding phases.

The chemical fractionation used includes the following five fractions:  F1: exchangeable fraction (magnesium chloride solution 1 mol.L-1 adjusted to pH 7.0);  F2: carbonate fraction (sodium acetate solution 1 mol.L -1 adjusted to pH 5.0 with acetic acid);

 F3: Fe-Mn oxide fraction (hydroxylammonium chloride 0.04 mol.L -1 in 25 % (v\v) acetic acid at 96 °C);

 F4: sulfide and organic matter fraction (nitric acid 0.02 mol.L -1 and hydrogen peroxide 8.8 mol.L -1 at 85 °C followed by ammonium acetate solution 3.2 mol.L -1 in 20 % (v\v) nitric acid); and,  F5: residual fraction. This fraction was estimated from the total content obtained using acid digestion.

Single extraction techniques

Single extraction techniques using 1 mol L -1 sodium hydroxide (NaOH) solution, and 0.1 mol L - 1 dipotassium hydrogen orthophosphate (K 2 HPO 4 ) solution were also performed to verify that an As adsorption on Fe oxyhydroxides was the major process of arsenic trapping in the soil [START_REF] Clozel | Physical, chemical, and biological treatability of contaminated soils (Final report -BRGM/RP-52065-FR)[END_REF]. These single extractions were carried out in triplicate using an LS ratio of 10 mL/g and a contact time of 48 hours.

Leachates analyses

All leachates were filtered through 0.45 µm pore size acetate-cellulose Whatman ® membranes. They were subsequently analyzed for concentrations of constituents of interest using ICP-AES.

Sulphates and chlorides were analyzed using IC (ion chromatography, DionexDX320 HPIC ® ).

Total organic carbon (TOC) in solution was measured using a TOC analyzer (Total Organic Carbon analyzer, Shimadzu TOC-5000A ® ) [START_REF] Clozel | Physical, chemical, and biological treatability of contaminated soils (Final report -BRGM/RP-52065-FR)[END_REF][START_REF] Chatain | Arsenic behavior in mining soil[END_REF][START_REF] Chatain | Characterization of arsenic and other inorganic constituents potential mobilization from soils collected from a gold mining site[END_REF]Chatain et al., 2005a and b;[START_REF] Bayard | Mobilisation of Arsenic from a mining soil in batch slurry experiments under bio-oxidative conditions[END_REF].

RESULTS AND DISCUSSION

Chemical characterization

The major physicochemical characteristics and the elemental total contents as determined by acidic digestion of the studied soil are summarized in Table 1.

The soil was mainly characterized by a moisture content of ca. 10 %, ca. 3 wt% arsenic, ca. 9 wt% iron and ca. 1.9 wt% total organic compounds. Its natural pH and ORP measured after 48 hours of contact with deionized water at a L/S ratio of 10 mL/g, were ca. 6.5 and ca. +400 mV vs. NHE, respectively. Moreover, a very low soluble (ca. 0.016 %) fraction of arsenic were obtained from the single extraction procedure using deionized water as the extractant. As (mg kg -1 ) 27,700

Cu (mg kg -1 ) 1,700

Pb (mg kg -1 ) 800

Zn (mg kg -1 ) 400

Mineralogical analysis

The methodology used scanning electron microscopy coupled with energy dispersive spectroscopy analyzer and X-ray diffraction analysis.Results from XRD analysis (Table 2) indicated that the soil was mainly constituted of quartz, gypsum, feldspars, muscovite, calcite, and dolomite. XRD analyses indicated the presence of traces of hematite and jarosite, which are the common pyrite weathering products with which As can be typically bound [START_REF] Clozel | Physical, chemical, and biological treatability of contaminated soils (Final report -BRGM/RP-52065-FR)[END_REF]. The absence of residual mixed sulfide phases in the soil can be due to the weathering and oxidation of pyritic minerals during the post-mining period.

Microscopy analyses performed using the BSE images by SEM coupled with EDS spectroscopy microanalysis can be used to complement the XRD analysis. The major minerals that were identified are quartz, gypsum, and particles containing Fe and As, with a composition in accordance with XRD results. A close similarity between the arsenic and iron cartography was observed in SEM results (Figure 1). Fe and As, highly present in the soil (ca. 3% and ca. 9%, respectively), were not located in any of the main crystallized phases. As seemed to have been sorbed with Fe on the soil particles. The presence of amorphous As-bearing phases, such as amorphous Fe oxyhydroxides, was therefore assumed [START_REF] Chatain | Arsenic behavior in mining soil[END_REF][START_REF] Chatain | Characterization of arsenic and other inorganic constituents potential mobilization from soils collected from a gold mining site[END_REF]Chatain et al., 2005a and b;[START_REF] Bayard | Mobilisation of Arsenic from a mining soil in batch slurry experiments under bio-oxidative conditions[END_REF] and thought to be the major process of arsenic trapping in the soil.

Leaching behavior

Sequential chemical extraction

Results of sequential extractions are summarized in Figure 2. As and Fe appeared to be strongly bound to the soil matrix (ca. 80 wt% were concentrated in the residual fraction). A very low fraction of As and Fe (i.e., ca. < 1 wt%) was extracted as exchangeable and bound-to-carbonates fractions. These results confirmed that As and Fe were strongly bound on soil amorphous oxyhydroxide particles as shown already with the results from the XRD and SEM study.

Moreover, Al appears as the least mobile metal because of its high presence in the residual fraction (F5ca. 97 wt%). No correlations with As or Fe can be established.

Single extractions

Results obtained from single extractions, using 1 mol L -1 NaOH solution, and 0.1 mol L -1 K 2 HPO 4 solution,are summarized in Table 3. They are compared with those acquired during single extraction using demineralized water (natural). Concerning NaOH extraction, about 74 % of the total content in As was extracted within these extreme conditions of pH (ca. 14) achieved. A mechanism of As desorption from Fe oxyhydroxide surface by competition with OH-ions might explain the abruptly increase in As extractability [START_REF] Carbonell-Barrachina | The influence of redox chemistry and pH on chemically active forms of arsenic in sewage sludge-amended soil[END_REF]. On the opposite, Fe solubility remained very low (ca. 0.2% of the total Fe content). Fe solubility may have been controlled by Fe oxyhydroxide precipitates.

These results suggest that most of As present in the soil is sorbed onto Fe-oxyhydroxides, which is consistent with mineralogical analysis. The remaining fraction (ca. 26% of the total As content) would represent the part co-precipitated with Fe and/or bound to resistant compounds (i.e., silicates or sulfides).

Results of extraction of As from the contaminated soil by K 2 HPO 4 solution are also shown in Table 3. The addition of this solution was found to have no significant effect on pH and ORP conditions and Fe release, as compared to those obtained under "natural conditions" with deionized water. However, a drastic increase in As solubilization (as much as 400-fold at pH ca. 7), and in phosphates adsorption onto Fe oxyhydroxides (ca. 76%) was evidenced. These results are consistent with studies performed by [START_REF] Alam | Extraction of arsenic in a synthetic arsenic-contaminated soil using phosphate[END_REF][START_REF] Clozel | Physical, chemical, and biological treatability of contaminated soils (Final report -BRGM/RP-52065-FR)[END_REF], that confirmed that As is for the most part reversibly sorbed onto Fe oxyhydroxides.

CONCLUSIONS AND PERSPECTIVES

In this research, an As contaminated soil collected from a gold mining site has been considered. The assessment of the operational fractionation of As present in this soil has been performed by using leaching tests in conjunction with mineralogical characterization.

The first results shown in this study indicate that despite a high total concentration of As in the studied soil, this release with deionized water is very limited. This was principally explained by its low solubility and by the stability of its bearing solid phases under natural slightly acidic conditions.

The multiple approach (leaching and mineralogy) used in this work provided information on the nature of As fixation processes in the mining soil studied, and a understanding of geochemical processes that control As reactivity in the subsurface environment.

Indeed, results showed that most of As appeared reversibly sorbed (ca. 72%) with Fe on the soil particles, in particular Fe oxyhydroxides. These results confirmed a potential risk of As mobilization over the long term under specific leaching conditions (i.e. pH or ORP gradient, chelation effect, etc.) which frequently occur in mining environment.

This initial study will lead to further research aimed at answering or confirming the processes considered. Metal speciation must be investigated further to provide the data needed for a better knowledge of As-bearing phases to better predict the mobility of this pollutant and in turn improve the management of arsenic contaminated sites.

Figure 1 -

 1 Figure 1 -SEM-EDX mapping of the 200 µm-2 mm fraction of the arsenic contaminated soil (12X)

Figure 2 -

 2 Figure 2 -Results fromsequential chemical extractions performed on the fines < 2 mm of the soil (conc. in wt%)

Table 1

 1

	-Physicochemical characteristics (LS ratio: 10 mL/g; contact time:
	48 hours) and total content of trace elements as determined by acid digestion
	of the studied soil
	Moisture content (%)	10
	pH (H 2 O)*	6.5
	ORP (mV vs. NHE)*	400
	Si (%)	20.3
	Fe (%)	8.9
	Al (%)	4.3
	Ca (%)	3.7

Table 2 -

 2 Mineral composition by XRD analysis (Abundance: ++++ very abundant; +++ abundant; ++ present ; + traces)

	Mineral	Abundance
	Quartz	++++
	Micas	++
	Feldspar	++
	Chlorite	+
	Kaolinite	+
	Calcite	++
	Dolomite	++
	Calcite	++
	Dolomite	++
	Gypsum	+++
	Jarosite	++

Table 3 -

 3 Results from single extractions (LS ratio: 10 mL/g; contact time: 48 hours)

	Extracting solution	As (conc. in wt%)	Fe (conc. in wt%)	pH	ORP (mV vs. NHE)
	Deionized water	0.0016	0.0001	6.5	+400
	1 mol L -1 NaOH solution	73.51	0.17	13.8	+60
	0.1 mol L -1 K 2 HPO 4 solution	4.89	0.0034	7.5	+430

ACKNOWLEDGMENTS

The research presented in this paper was supported by the Association RE.CO.R.D. (Waste Research Cooperative Network, France), and the Région Rhône-Alpes.


Anastasia Volkova

Thibault Hilaire

Christoph Lauter

Determining Fixed-Point Formats for a Digital Filter Implementation using the Worst-Case Peak Gain measure

HAL is

I. Introduction

The implementation of a digital filter can be considered as a path from filter specifications to a physical implementation (DSP device). This process includes generation of a transfer function, which defines the relation between the inputs and outputs of a filter. Further one needs to choose a filter structure to be eventually implemented, i.e. a computational scheme, a specific order of filter evaluation. Finally, there is software and hardware code generation.

However, on each step various issues arise. First of all, every transfer function must be discretized, which implies the modification of filter. Secondly, while in infinite precision filter structures (e.g. Direct Forms [START_REF] Oppenheim | Digital Signal Processing[END_REF], State-Space [START_REF] Friedland | Control Systems Design: An Introduction to State-Space Methods[END_REF], Wave [START_REF] Fettweiss | Wave digital filters: Theory and practice[END_REF] etc.) are equivalent, they are no longer the same objects in finite precision. Since the computational schemes are different, the round-off errors and consequently quality of structures vary. Finally, software and hardware implementations are performed under constraints, e.g. power consumption, area, output error, etc. And on each step of the filter implementation the degradation due to finite precision must be taken into account in order to produce a reliable implementation.

To unify the above-described process of filter implementation for any Linear Time-Invariant (LTI) filter and provide reliable implementation, we develop an automatized filter generator. Its work-flow is described in the Fig. 1.

In this paper we focus on one of the various steps: the Fixed-Point (FxP) Algorithm generation step and wordlength optimization. On this step we already know target specifications This work has been sponsored by the Agence nationale de la recherche grant ANR-13-INSE-0007-02 MetaLibm. and software implementation constraints. Therefore, during this optimization, we naturally look for a trade-off between area (for hardware implementation) and computational errors. If we take more bits, the implementation cost arises. If we take less bits, there exist a risk of an overflow. Depending on the application, one or even both criteria can be neglected. However, in the general case on each step of the optimization procedure we need to determine rigorously the Fixed-Point Formats (FxPF), i.e. the most and the least significant bits positions, for all variables in the filter implementation. These formats must allow no overflow and take into account computational errors.

There exist various approaches on determining FxPF, such as Affine Arithmetic [START_REF] Lopez | Improved intervalbased characterization of fixed-point LTI systems with feedback loops[END_REF], or the common approach of numerous simulations [START_REF] Kim | Fixed-point optimization utility for C and C++ based digital signal processing programs[END_REF]. However, these methods are not suitable for our needs since they do not give any mathematical guarantee on the absence of overflow for any possible input and tend to overestimate. We, on the other hand, propose a rigorous algorithm completely based on mathematical proofs.

In Section II we describe the fixed-point arithmetic in question and theoretical background for our approach. In Section III we define the problem for a filter in a statespace representation and introduce a two-step algorithm. Error analysis of the algorithm is presented in Section IV. Finally, numerical results are provided before conclusion.

Notation: Throughout the article matrices are in uppercase boldface, vectors are in lowercase boldface, scalars are in lowercase. All matrix inequalities and absolute values are considered element-by-element. The matrix 1 denotes a matrix of ones with the appropriate size.

m + 1 - s w -2 m 2 0 2 -1 2 m-1 2
Fig. 2. Fixed-point representation (here, m = 5 and = -4).

II. Basic bricks

A. The Worst-Case Peak Gain theorem

Let H := (A, B, C, D) be a Bounded-Input Bounded-Output stable MIMO LTI filter in state-space representation:

H x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) + Du(k) (1)
where u(k) ∈ R q×1 is the input vector, y(k) ∈ R p×1 the output vector, x(k) ∈ R n×1 the state vector and A ∈ R n×n , B ∈ R n×q , C ∈ R p×n and D ∈ R p×q are the state-space matrices. We will use the following approach on deducing the output interval for a filter: suppose all the inputs are guaranteed to be in a known interval:

∀k, |u i (k)| ūi , i = 1, . . . , q.
(2)

Then, there exists a N 0 such that ∀k ≥ N 0 the output y(k) lies in the interval:

|y(k)| H ū, (3)
where H denotes the so-called Worst-Case Peak Gain (WCPG) matrix of the filter [START_REF] Hilaire | Reliable implementation of linear filters with fixed-point arithmetic[END_REF]. It can be computed as the 1 -norm of the impulse response of the filter. The N 0 can be determined using reasoning similar to the WCPG error analysis in [START_REF] Volkova | Reliable evaluation of the worst-case peak gain matrix in multiple precision[END_REF]. There always exists an input sequence u(k) such that the equality in (3) is reached.

There exist numerous other approaches to determine the output interval, such as involving Affine Arithmetic [START_REF] Lopez | Improved intervalbased characterization of fixed-point LTI systems with feedback loops[END_REF], which may lead to overestimations, or simulations [START_REF] Kim | Fixed-point optimization utility for C and C++ based digital signal processing programs[END_REF], which in turn are inefficient and not rigorous.

It can be shown that the WCPG approach gives the smallest interval containing all the possible values of y(k) and that the bound (3) can be attained. Moreover, in [START_REF] Volkova | Reliable evaluation of the worst-case peak gain matrix in multiple precision[END_REF] the authors have published an algorithm to evaluate the WCPG at arbitrary precision.

B. Fixed-Point arithmetic for filter implementation

The considered filters are implemented using signed FxP Arithmetic, with two's complement representation [START_REF] Newmann | First draft of a report on the edvac[END_REF], [START_REF] Finley | Two's complement[END_REF]. Let t be such a FxP number. It can be written as t = -2 m t m + m-1 i= 2 i t i , where its FxPF (m, ) gives the position of the most (MSB) significant and least (LSB) significant bits respectively, and t i ∈ B = {0, 1} is the i th bit of t as shown in Fig. 2. The wordlength w is given by w = m -+ 1. The quantization step of the representation is 2 . Further in the article the LSB position sometimes will be substituted directly with mw + 1.

The real number t is represented in machine by the integer T , such that

T = t • 2 -, where T ∈ [-2 w-1 ; 2 w-1 -1] ∩ Z. Note that this interval is not symmetric.
Let y(k) ∈ R be an output of a digital filter. Given the wordlength w y , determining the Fixed-Point Formats for y(k) means to find a MSB vector m y such that for each k

y i (k) ∈ [-2 m y i ; 2 m y i -2 m y i -w y i +1 ]. (4)
Obviously, we are interested in the least possible MSB, which guarantees (4): any greater MSB would yield to a larger quantization step, hence round-off noise, at constant wordlength.

III. Determining the Fixed-Point Formats

A. Exact filter

The problem of determining the FxPF for a filter H implementation can be formulated as follows. Let H be a filter in a state-space representation [START_REF] Oppenheim | Digital Signal Processing[END_REF]. Suppose all the inputs to be in an interval bounded by ū.

Given the wordlength constraints vector w x for the state and w y for the output variables we look for a FxPF for x(k) and y(k) such that for any possible input no overflow occurs. We wish to determine the least MSB vectors m y and m x such that

∀k, y(k) ∈ [-2 -m y ; 2 m y -2 m y -w y +1 ], (5)
∀k, x(k) ∈ [-2 -m x ; 2 m x -2 m x -w x +1 ]. ( 6
)
Remark 1. Since the filter H is linear and input interval is centered at zero, the output interval is also centered in zero. Therefore, further we will seek to determine the least m y and

m x such that ∀k, |y(k)| 2 m y -2 m y -w y +1 , (7)
∀k, |x(k)| 2 m x -2 m x -w x +1 . (8)

B. Applying the WCPG to compute MSB positions

Applying the WCPG theorem on the filter H yields a bound on the output interval:

|y i (k)| ( H ū) i , i = 1, . . . , p. (9)
Denote the bound vector ȳ := H ū. We can determine the FxPF for the output of a LTI filter H utilizing the following lemma.

Lemma 1. Let H = (A, B, C, D) be a BIBO-stable MIMO LTI filter and ū be a bound on the input interval. Suppose the wordlengths w y are known and w y i > 1, i = 1, . . . , p.

If for i = 1, . . . , p the MSBs are computed with

m y i = log 2 (ȳ i ) -log 2 1 -2 1-w y i (10)
and the LSBs are computed with y i = m y i + 1 -w y i , then for all k |y i (k)| 2 m y i -2 m y i -w y i +1 and m y i is the least.

Proof. We look for the least m y such that (7) holds. Since the bound ȳ can be reached, it is sufficient to require:

ȳi 2 m y i -2 m y i -w y i +1 (11)
Solving this inequality for m y i we obtain that the smallest integer, which satisfies the above inequality is

m y i = log 2 (ȳ i ) -log 2 1 -2 1-w y i . (12)
C. Modification of filter H to determine bounds on the state variable x Using Lemma 1 we can determine the FxPF for the output of a filter. In order to determine the FxPF for the state variable we modify the filter H by vertically concatenating the state vector and the output and include necessary changes into the state matrices.

Denote vector ζ(k) := x(k) y(k) to be the new output vector.

Then the state-space relationship (1) takes the form:

H ζ          x(k + 1) = Ax(k) + Bu(k) ζ(k) = I C x(k) + 0 D u(k) . ( 13
)
Hence the problem is to find the least MSB vector m ζ such that (element-by-element)

∀k, |ζ(k)| 2 m ζ -2 m ζ -w ζ +1 . ( 14
)
Now, applying the WCPG theorem on the filter H ζ and using Lemma 1, we can deduce the MSB positions of the state and output vectors for an implementation of the filter H.

D. Taking rounding errors into account

However, due to the finite-precision degradation what we actually compute is not the exact filter H ζ but an implemented filter H ♦ ζ :

H ♦ ζ            x ♦ (k + 1) = ♦ x Ax ♦ (k) + Bu(k) ζ ♦ (k) = ♦ ζ I C x ♦ (k) + 0 D u(k) (15)
where the Sums-of-Products (accumulation of scalar products on the right side) are computed with some rounding operator ♦ . Suppose, this operator ensures faithful rounding [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]:

|♦ (x) -x| < 2 (16)
In [START_REF] Lopez | Formatting bits to better implement signal processing algorithms[END_REF], [START_REF] De Dinechin | Sum-of-product architectures computing just right[END_REF] it was shown that such an operator can be implemented using some extra guard bits for the accumulation.

Denote the errors due to operator ♦ as ε x (k) and ε y (k) for the state and output vectors, respectively. Essentially, the vectors ε x (k) and ε y (k) may be associated with the noise which is induced by the filter implementation. Then the implemented filter can be rewritten as

H ♦ ζ          x ♦ (k + 1) = Ax ♦ (k) + Bu(k) + ε x (k) ζ ♦ (k) = I C x ♦ (k) + 0 D u(k) + 0 I ε y (k) , (17)
where

|ε x (k)| < 2 x , ε y (k) < 2 y .
It should be remarked that since the operator ♦ l is applied ε x (k) x(k) -x ♦ (k) and ε y (k) y(k) -y ♦ (k). As the rounding also affects the filter state, the x ♦ (k) drifts away from x(k) over time, whereas with ε x (k) we consider the error due to one step only.

It can be observed that at each instance of time the state and output vectors are computed out of u(k) and errorvectors, which can be considered as inputs as well. Thanks

H ⇣ H ⇣ ⌃ (k) u(k) ⇣(k) ⇣ (k) ✓ "x(k) "y(k) ◆ m ⇣ Fig. 3. Implemented filter decomposition.
to the linearity of the filters, we can decompose the actually implemented filter into a sum of the exact filter and an "errorfilter" H ∆ as shown in Fig. 3. Note that this "error-filter" is an artificial one; it is not required to be "implemented" by itself and serves exclusively for error-analysis purposes.

The filter H ∆ is obtained by computing the difference between H ♦ ζ and H ζ . This filter takes the rounding errors

ε(k) := ε x (k) ε y (k)
as input and returns the result of their propagation through the filter:

H ∆                ∆ x (k + 1) = A∆ x (k) + I 0 ε(k) ∆ ζ (k) = I C ∆ x (k) + 0 0 0 I ε(k) , (18)
where, the vector ε(k) is guaranteed to be in the interval bounded by ε := 2 ζ . Once the decomposition is done, we can apply the WCPG theorem on the "error-filter" H ∆ and deduce the output interval of the computational errors propagated through filter:

∀k, ∆ ζ (k) H ∆ • ε. (19)
Hence, the output of the implemented filter is bounded with

ζ ♦ (k) |ζ(k)| + ∆ ζ (k) . (20)
Applying Lemma 1 on the implemented filter and using (20) we obtain that the MSB vector m ♦ ζ can be upper bounded by

m ♦ ζ i = log 2 H ζ • ū i + ( H ∆ • ε) i (21) -log 2 1 -2 1-w ζ i .
Therefore, the FxPF (m ♦ ζ , ♦ ζ ) guarantee that no overflows occur for the implemented filter.

Since the input of the error filter H ∆ depends on the FxPF chosen for implementation, we cannot directly use (21). The idea is to first compute the FxPF of the variables in the exact filter H, where computational errors are not taken into account, and use it as an initial guess for implemented filter H ♦ ζ . Hence, we obtain the following two-step algorithm:

m := log 2 H ζ • ū i + ( H ∆ • ε) i -log 2 1 -2 1-w ζ i .
Handling floating-point analysis of multiplications and additions in ( 21) is trivial using approach by Higham [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF]. The difficulty comes from the WCPG matrices, which cannot be computed exactly. However both approximations H ζ and H ∆ , even if computed with arbitrary precision, bear some errors ε WCPG ζ and ε WCPG ∆ that satisfy

0 H ∆ -H ∆ ε WCPG ζ • 1 (22) 0 H ζ -H ζ ε WCPG ∆ • 1 (23)
Introducing the errors on the WCPG computations into the formula (21) we obtain that what we actually compute is

m ζ i ♦                    m + log 2                 1 + ε WCPG ζ q j=1 ū j + ε WCPG ∆ n+p j=1 ε j H ζ ū i + ( H ∆ ε) i                                    . (24)
The error term in (24) cannot be zero (apart from trivial case with zero ū). However, since we can control the accuracy of the WCPG matrices, we can deduce conditions for the approximation m ζ i ♦ to be off by at most one.

ε WCPG ∆ < 1 2
H ∆ • ε i p+n j=1 εi (25) ε WCPG ζ < 1 2 H ζ • ū i q j=1 ūi , ( 26
)
where

H := | D| + |CB| + |C AB|, then 0 m ζ i ♦ -m ζ i ♦ 1. ( 27
)
Proof. Proof by construction, we reason as follows: since the error-term caused by the WCPG floating-point evaluation is positive and the ceil function is increasing, then

m ζ i ♦ -m ζ i ♦ 0, (28)
i.e. the floating-point approximation m ζ i ♦ is guaranteed to never be underestimated. However, it can overestimate the MSB position by

m ζ i ♦ -m ζ i ♦                 m -m -1<• 0 + log 2                 1 + ε WCPG ζ q j=1 ū j + ε WCPG ∆ n+p j=1 ε j H ζ ū i + ( H ∆ ε) i                                    . (29)
The approximation m ζ i ♦ overestimates at most by one bit if and only if the error term is contained in the interval [0, 1), i.e. if

0 log 2                 1 + ε WCPG ζ q j=1 ū j + ε WCPG ∆ n+p j=1 ε j H ζ ū i + ( H ∆ ε) i                 < 1. (30)
Hence, using the above condition we can deduce the upper bounds on the ε WCPG ζ and ε WCPG ∆ :

0 ε WCPG ζ q j=1 ū j + ε WCPG ∆ n+p j=1 ε j H ζ ū i + ( H ∆ ε) i < 1. (31)
Since all the terms are positive, the left inequality is always true. The right inequality in (31) is satisfied for instance if

ε WCPG ζ q j=1 ū j H ζ • ū i < 1 2 ε WCPG ∆ n+p j=1 ε j ( H ∆ • ε) i < 1 2 . ( 32
)
Rearranging terms we obtain following inequalities on the WCPG computation with error:

ε WCPG ζ < 1 2 • H ζ ū i q j=1 ū j ε WCPG ∆ < 1 2 • ( H ∆ ε) i n+p j=1 ε j . (33)
Unfortunately, the above results cannot be used in practice, since they depend themselves on the exact WCPG matrices.

It can be shown that H is a lower bound of the WCPG matrix. We can compute this matrix exactly. Obviously,

H ∆ • ε i p+n j=1 εi ( H ∆ • ε) i p+n j=1 εi (34) H ζ • ū i q j=1 ūi H ζ • ū i q j=1 ūi . ( 35
)
Hence, if the WCPG matrices in the right sides of (33) are substituted with their lower bounds, the condition (31) stays satisfied and we obtain bounds (25) and (26).

Analogously, Lemma 2 can be applied to the computation of m ζ i with the terms concerning filter H ∆ set to zero.

In order to guarantee the output MSB position to be optimal, an instance of the Table Maker's Dilemma [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF] must be solved. This is due to the simplification by using the triangle inequality in (20), which basically requires both filters H ζ and H ∆ reach the worst-case bound with the same input sequence.

V. Complete algorithm

The two-step algorithm, presented in subsection III-D takes into account accumulation of computational errors in a filter over time and Lemma 2 presents error-analysis of the MSB position computation procedure. However, one additional fact has not been taken into account.

In most cases the MSB vectors m ζ (computed on Step 1) and m ♦ ζ (computed on Step 2) are the same. However, in some cases they are not, which can happen due to one of the following reasons:

VI. Numerical results

The above described algorithm was implemented as a C library, using GNU MPFR1 [14] version 3.1.12, GNU MPFI2 version 1.5.1 and the WCPG library [START_REF] Volkova | Reliable evaluation of the worst-case peak gain matrix in multiple precision[END_REF].

Consider following example: let H be a random stable filter with 3 states, 1 input and 1 output. Suppose the inputs are in an interval absolutely bounded by ū = 5.125. In this example we set all the wordlength constraints to 7 bits, however our library supports the multiple wordlength paradigm.

The results of the work of our algorithm can be observed in Table I. On the Step 1 we deduce the MSB positions for the filter H, which does not take computational errors into account. These MSBs serve as an initial guess for the Step 2. We compute the vector εζ and construct the "error-filter" H ∆ . Taking into account the error propagation yields changes in the MSB positions: the rounding errors force the third state x 3 (k) to pass over to the next binade. However, moving the MSB yields larger quantization step and an addition step is required. Step 3 verifies that the FxPF deduced on the Step 2 guarantees no overflow. We can verify the result by tracing the state and output vectors in two cases: the exact vectors y(k) and x(k); and the quantized vectors y ♦ (k) and x ♦ (k). The quantized vectors are the result of an implementation of the filter H with FxPF deduced on the Step 1. For each state and output we take an input sequence which makes the state respectively output attain the theoretical bound. However, this sequence is not guaranteed to maximize the quantized error.

In Fig. 4 we can observe that the exact vector y(k) attains the bound computed with the WCPG theorem. The quantized output stays within the bound ȳ♦ for the implemented filter but passes over the bound of the exact filter. However, as the bound ȳ is far from the next binade, the MSB position is not changed and the computational errors do not result in overflow.

For the third state x 3 (k), its bound x3 is very close to the 2 m x 3 -2 x 3 and taking into account the filter computational errors yields passing to the next binade. It can be clearly observed in Fig. 5 that the quantized state passes over the value 2 m x 3 -2 x 3 . If the computational errors were not taken into account, the initial format (deduced on the Step 1) would result in overflow.

Therefore, we observed that deducing the FxPF without taking into account the computational errors can result in overflow but our approach guarantees a reliable implementation without overestimation.

VII. Conclusions

We give an algorithm to determine the FxPF for all variables of a LTI filter, where the computation errors are taken into account. We ensure, by construction, that no overflow occurs. However, the computed MSB positions can overestimate at most by one bit. In order to guarantee no overestimation at all an instance of the TMD must be solved. Our algorithm can be extended to any LTI filter realization (Direct Forms, Lattice filters, etc.) using the Specialized Implicit Framework (SIF) tool [6]. SIF enables unification of various LTI filter realizations for analysis, comparison and FxP code generation. Therefore, this work is an essential basic brick for the automatic filter-to-code conversion.

Fig. 1 .

 1 Fig. 1. Automatic filter generator flow.

Step 1 :

 1 Determine the FxPF (m ζ , ζ ) for the exact filter H ζ Step 2: Construct the "error-filter" H ∆ , which gives the propagation of the computational errors induced by format (m ζ , ζ ); then, compute the FxPF (m ♦ ζ , ♦ ζ ) of the actually implemented filter H ♦ ζ using (21). The above algorithm takes into account the filter implementation errors. However, the algorithm itself is implemented in finite-precision and can suffer from rounding errors, which influence the output result. All operations in the MSB computation will induce errors, so what we actually compute are only floating-point approximations m ζ and m ζ ♦ . In what follows, we propose error-analysis of the floating-point evaluation of the MSB positions via (10) and (21). IV. Error analysis of the MSB computation formula Let us consider case of m ζ i ♦ and show afterwards that m ζ i is its special case. To reduce the size of expressions, denote

Lemma 2 .

 2 If the WCPG matrices H ζ and H ∆ are computed such that (22) and (23) hold with

•

  the accumulated rounding error due to the FxPF ( m ζ , ζ ) makes output of the actually implemented filter pass over to the next binade; or • the floating-point approximation m ♦ ζ is off by one. Moreover, if the MSB position after Step 2 of the algorithm is increased, the LSB position moves along and increases the error. Therefore, the modified format must be re-checked. Hence, the FxPF determination algorithm gets transformed into the following three-step procedure: Step 1: Determine the FxPF ( m ζ , ζ ) for the exact filter H ζ ; Step 2: Construct the "error-filter" H ∆ , which shows the propagation of the computational errors induced by format ( m ζ , ζ ); then, compute the FxPF ( m ♦ ζ , ♦ ζ ) of the actually implemented filter H ♦ ζ ; Step 3: If m ♦ ζ i == m ζ i , then return ( m ♦ ζ , ♦ ζ ); otherwise m ζ i ←-m ζ i + 1 and go to Step 2.

Fig. 4 .

 4 Fig. 4. The exact and quantized outputs of the example. Quantized does not pass over to the next binade.

3 Fig. 5 .

 35 Fig. 5. The exact and quantized third state of the example. Quantized output passes over to the next binade.

http://www.mpfr.org/

https://gforge.inria.fr/projects/mpfi/