Jean-Daniel Bontemps

Jean-Christophe Hervé

Jean-François Dhôte

Long-term changes in forest productivity: a consistent assessment in even-aged stands

Keywords: growth trends, dominant height, statistical modeling, mixed-effects models, Fagus sylvatica

The objective of the study was to provide a stand-scale assessment of long-term productivity changes in even-aged stands and to depict their chronology over the twentieth century.

We focused on dominant height growth as a proxy for productivity, reconstructed from stem analyses in temporary plots. Height increments from two generations of stands were compared.

Stands were associated in pairs to ensure accurate control of intrinsic site fertility conditions. The historical evolution of the growth rate was estimated using a statistical modeling approach based on a mixed-effects model, with a control of both site and developmental stage effects. We placed emphasis on a model formulation that leads to a meaningful interpretation of growth rate evolution.

We applied the methodology to a sample of 14 stand pairs and 84 stem analyses of Common Beech in northeastern France. An accelerated increase in growth rate was identified, reaching +50% over the twentieth century. The trend also displayed short-term growth anomalies.

INTRODUCTION

Increases in the growth of forest tree species have been acknowledged in the Northern Hemisphere and in other parts of the world (Spiecker et al. 1996[START_REF] Jacoby | Tree rings, carbon dioxide, and climatic change[END_REF][START_REF] Boisvenue | Impacts of climate change on natural forest productivityevidence since the middle of the 20 th century[END_REF]. With increasing consensus on their existence, emphasis has then been placed on the environmental causality of such changes (Nelleman and Thomsen 2001[START_REF] Nabuurs | Stemwood volume increment changes in European forests due to climate changea simulation study with the EFISCEN model[END_REF][START_REF] Milne | A comparison of two modelling studies of environmental effects on forest carbon stocks across Europe[END_REF][START_REF] Kahle | Causes and consequences of forest growth trends in Europe[END_REF]. However, the magnitude and historical dynamics of these increases remain poorly understood to date [START_REF] Spiecker | Growth trends in European forests -Do we have sufficient knowledge? P. 157-169 in Causes and consequences of accelerating tree growth in Europe[END_REF]. The main reasons for this lie in recurrent deficiencies of growth-trend analyses: a choice of inconsistent productivity indicators, a limited or partial control in sampling designs of the factors which usually affect productivity, and a lack of a combined modeling strategy for the estimation of growth change.

Productivity indicators. Forest productivity is defined as a wood biomass or volume increment per unit of time and area (Assmann 1970, P. 160-161). Because most forests undergo regular but undocumented thinning operations, there is inherent difficulty in determining forest productivity from repeated forest inventory measurements. Noticeable exceptions arise from a limited number of permanent plot networks, covering a time period ranging from several decades to a century [START_REF] Eriksson | Yields of Norway spruce (Picea abies (L.) Karst.) in two consecutive rotations in southwestern Sweden[END_REF][START_REF] Pretzsch | Growth trends of forests in Southern Germany, P. 107-131[END_REF][START_REF] Elfving | Trends of tree growth in Swedish forests 1953-1992: an analysis based on sample trees from the National Forest Inventory[END_REF]. Hence, the only way to address growth-trend issues is to consider past growth of trees using retrospective techniques, such as ring width measurements from increment cores (Cook and Kairiukstis 1990, P. 26-28), or past height measurements from stem analysis [START_REF] Curtis | A stem analysis approach to site index curves[END_REF]. Because increment cores are far easier and cheaper to sample, radial growth is the first ranking data source on the topic (Spiecker et al. 1996, Jacoby and[START_REF] Jacoby | Tree rings, carbon dioxide, and climatic change[END_REF]. Nevertheless, there are weaknesses inherent in radial growth indicators: (i) they depend on between-tree competition, including the effects of silvicultural schedules or local disturbances, (ii) their measurement at breast height may underestimate tree volume increment [START_REF] Bouriaud | Is ring width a reliable proxy for stem-biomass increment? A case study in European beech[END_REF]. In contrast, dominant height is not affected by competition [START_REF] Frothingham | Height growth as a key to site[END_REF][START_REF] Lanner | On the insensitivity of height growth to spacing[END_REF] and has been used as a traditional proxy for productivity of even-aged stands [START_REF] Eichhorn | Relationships between dominant height and stand volume (in German)[END_REF], Assmann 1970, P. 161-163). The dominant height of stands can be reconstructed from stem analysis if applied to dominant trees. However, stem analysis is destructive and time consuming, and very few growth-trend studies based upon this proxy have been published [START_REF] Untheim | Has site productivity changed? A case study in the Eastern Swabian Alb, Germany[END_REF][START_REF] Kahle | Causes and consequences of forest growth trends in Europe[END_REF].

Sampling designs. The context of retrospective analyses renders possible the conception of sampling designs well-grounded for the investigation of growth trends. They must permit comparison of growth intensities during distinct periods in time, with control of intrinsic site fertility conditions. This requires selecting stands of different generations, with a sufficiently large time range to identify long-term trends [START_REF] Spiecker | Growth trends in European forests. Has site productivity changed?[END_REF]. A classical approach suggests that the balance of age and fertility factors is ensured through regional samples of large size (Badeau et al. 1996, Schneider and[START_REF] Schneider | Growth trends of trees. Regional study on Norway spruce (Picea abies[END_REF]. However, rotation age usually increases with decreasing site fertility. Consequently, the oldest stands are mostly located on the poorest sites whereas a wider range of sites is covered by the younger ones. This may result in a positive bias in the detected longterm trend [START_REF] Becker | Long-term changes in forest productivity in northeastern France: the dendrocological approach, P. 143-153 in Forest decline and atmospheric deposition effects in the French mountains[END_REF]. The latter flaw can be avoided by sampling clusters of stands of different ages, close to each other in space. This procedure grants very close proximityif not identityof topographic, soil, and mean climate conditions, but has been rarely used [START_REF] Untheim | Has site productivity changed? A case study in the Eastern Swabian Alb, Germany[END_REF][START_REF] Lebourgeois | Height and radial growth trends of Corsican pine in western France[END_REF][START_REF] Vejpustkova | Growth trends of spruce in the Orlické hory Mts[END_REF].

Statistical modeling approach. The assessment of growth trend magnitude and dynamics requires isolation of the historical signal from growth data, where it is confused with site and age, and with competition effects if the growth indicator selected depends on competition. As a quantitative and multi-dimensional problem, it should be addressed through a statistical modeling approach.

However, the aim of most earlier studies was to detect rather than to quantify growth trends, and very few of them were based on a modeling approach (Elfving and Nyström, Tome et al., Mielikäinen and Timonen, among contributions gathered in Spiecker et al. 1996[START_REF] Elfving | Trends of tree growth in Swedish forests 1953-1992: an analysis based on sample trees from the National Forest Inventory[END_REF][START_REF] Goelz | Long-term growth trends of red spruce and fraser fir at Mt. Rogers, Virginia and Mt. Mitchell, North Carolina[END_REF][START_REF] Gschwantner | Growth changes according to the data of the Austrian national forest inventory and their climatic causes[END_REF]). In addition, the estimation of the current date effect leads to considering short-term growth increments (annual to pluri-annual) instead of totals. Otherwise, the chronology estimate at a given date would integrate former changes recorded in growth and would not be interpretable (e.g. Elfving andTegnhammar 1996, Goelz et al. 1999).

The need for modeling approaches has been pointed out in [START_REF] Spiecker | Growth trends in European forests -Do we have sufficient knowledge? P. 157-169 in Causes and consequences of accelerating tree growth in Europe[END_REF], but no decisive step ahead has been accomplished in more recent studies [START_REF] Vejpustkova | Growth trends of spruce in the Orlické hory Mts[END_REF][START_REF] Perez | Growth trends and relationships with environmental factors for Scots pine [Pinus sylvestris (L.)] in Brandembourg[END_REF][START_REF] Lopatin | Long-term trends in height growth of Picea obovata and Pinus sylvestris during the past 100 years in Komi Republic (north-western Russia)[END_REF]). The quantitative description of long-term trends therefore remains a largely unexplored issue to date.

In France, the available contributions on changes in forest productivity mostly result from dendrochronological approaches to radial growth, and report positive trends [START_REF] Badeau | Long-term growth trends of trees: Ten years of dendrochronological studies in France[END_REF]. In this paper, we therefore revisit the subject of changes in productivity through a case study on even-aged stands of Common Beech (Fagus sylvatica L.) in the northeastern area. Our objective is to provide a consistent assessment of the chronology and intensity of such past changes. Dominant height growth was selected for its insensitivity to stand density and its relationship to stand productivity, and was reconstructed from stem analyses. The sampling design was defined according to the paired-plots method, based on neighboring stands of different ages located in the same site fertility conditions. The historical signal in dominant height growth was extracted using an original statistical modeling analysis of height increments, allowing to factor out the effects of site and developmental stage. Specific attention was paid to the interpretability of the model and of the historical signal extracted, which to our knowledge has never been proposed earlier.

MATERIALS

Sampling strategy: the paired-plots method

The identification of the environment-related historical changes in growth implies control of the developmental stage and intrinsic site fertility, which affect dominant height. The paired-plots sampling method [START_REF] Untheim | Has site productivity changed? A case study in the Eastern Swabian Alb, Germany[END_REF][START_REF] Lebourgeois | Height and radial growth trends of Corsican pine in western France[END_REF] consists in associating neighboring stands, which should differ only in age. It thus provides an efficient way to sample stands whose ages are balanced with permanent environmental conditions. Control of the site fertility can be further strengthened by additional measurements of local environmental indicators. Also, greater proximity in genetics (in natural stands) and management practices between the associated plots is an implicit benefit of the method.

We applied the method to pure and regular stands of Common Beech within its semi-continental production range in northeastern France. Beech was selected because management-related impacts on growth evolution are likely very weak: there is no existing breeding program, Beech is attested always to have been regenerated naturally [START_REF] Hüffel | Les méthodes de l'aménagement forestier en France[END_REF], and the intensity of management has been reported as traditionally moderate [START_REF] Polge | The influence of thinnings on the growth constraints of Beech (in French)[END_REF]. Moreover, the sampling was restricted to State forests to ensure structure and management continuity through time.

To vary site conditions within the sample, 14 stand pairs were sampled in 1998 within the considered geographic domain. Old stands were selected that were close to the standard rotation age (150-180 years). Young stands were chosen half as old (70-80 years) in order to ensure a correct statistical separation of ontogenic/environmental effects on growth and a comparison of growth within a sufficient age range. The prior control of within-pair environmental conditions included a maximum distance of 500 m, and the identity of topography, parent rock, and soil. Stand pairs were selected after an in situ check of stand characteristics, topography and soil (depth, texture, and stoniness), completed by humus forms [START_REF] Jabiol | Humus under its all forms[END_REF]) and understorey vegetation. For an accurate comparison of soil indicators, soil analyses were also carried out in 2007 (see below), and were completed by vegetation surveys conducted in the summer season in 400 m 2 area plots, to predict water and nutrient status from Ellenberg indicators [START_REF] Ellenberg | Zeigerwerten von pflanzen in Mitteleuropa[END_REF][START_REF] Hawkes | Using Ellenberg indicator values to assess soil quality in British forests from ground vegetation: a pilot study[END_REF]. Stand characteristics are reported in Table 1 andmapped in Figure 1

Dominant height sampling and reconstitution

The height growth trajectory of a tree can be reconstructed from its origin using stem analysis [START_REF] Curtis | A stem analysis approach to site index curves[END_REF], based upon ring counts of disks sawn at regular intervals along a tree stem, up to the terminal bud. When stem analysis is applied to a sample of dominant trees, the calculation of a dominant height growth trajectory is possible. Dominant height was defined as the mean height of the 100 thickest trees/ha at 1.30 m height. In each stand, dominant trees were sampled according to [START_REF] Duplat | Modelling dominant height growth of Sessile oak (Quercus petraea Liebl. in France. Regional variability and effect of the recent period (1959-1993)[END_REF]: to avoid local sampling bias [START_REF] Pierrat | Estimation de la moyenne des valeurs les plus élevées d'une population finie: Application aux inventaires forestiers[END_REF][START_REF] Garcia | Estimating top height with variable plot sizes[END_REF], the five thickest trees of 0.06-ha circular plots were selected, of which the first, third, and fifth ones were cut. Disks were first sawn every 4 m along logs, starting at a constant 0.30 m, then every 2 m on reaching the upper parts of stems and crowns. Age was counted from the stump height. In total, 84 stem analyses were performed.

These stem analyses provided successive height measurements separated by several years. To compute the mean trajectory of dominant trees, additional height-age pairs were calculated by linear interpolation of individual curves. As tree ages were seldom identical within a stand, curves were averaged by date to preserve growth synchronization over time. Calculation dates were chosen to include as many initial (instead of interpolated) measurements as possible, with an average resolution criterion of 7-10 years. Stand age was defined as the median tree age. Growth trajectories are plotted in Figure 2. Seven singularly small increments were detected (a fraction of the old-stand curve of pair 8 was affected, see Figure 2) and were removed. The number and duration of increments considered for modeling are reported in Table 2.

Post-sampling assessment of site control

The close control of mean climate in the paired-plots method is implicit. To further check the quality of pairing for factors varying rapidly in space, environmental indicators were computed from in situ measurements, soil analyses, and vegetation surveys. They were also used as an indication of site fertility variations between pairs. They comprised:

-soil indicators of nutrition and water availability, based on soil analyses carried out in each stand following recommendations by [START_REF] Gégout | Soil analysis in forest context: Phytoecologist choices within the frame of station typologies or scientific studies (in French)[END_REF]: (a) soil layer textures were analyzed to provide 1 m depth soil water capacity (SWC100) calculated from texture capacity estimates by [START_REF] Bruand | Use of class pedotransfer functions based on texture and bulk density of clods to generate water retention curves[END_REF] and stone layer contents, and (b) pH, cations (Ca, K, Mg, Na, Al, H, Fe, Mn, cmol+.kg -1 ), total nitrogen (N), organic carbon (C), and phosphorus (P) concentrations (g.kg -1 )

were measured in the first soil layer. Additional synthetic nutrition indicators were then computed:

C/N ratio for nitrogen availability, cation exchange capacity (CEC), concentration of bases (S = Ca + K + Mg + Na), and base saturation rate (S/CEC).

-Vegetation-based Ellenberg indicators of local water and nutrition status [START_REF] Ellenberg | Zeigerwerten von pflanzen in Mitteleuropa[END_REF] including: nitrogen availability (Nel), basicity (Rel), and humidity (Fel). These indicators are defined for botanic species of Europe and are expressed on a relative scale ranging from 1 to 9, reflecting the increase in factors. For each indicator, averages were computed over each vegetation survey.

Means and standard-deviations for the main environmental indicators are given in Table 3.

MODEL FORMULATION

Model structure

The model was initially formulated in the form of a differential equation (continuous time). It was then integrated with respect to time to be fitted to growth increments.

Under a stable environment, dominant height growth rate is usually expressed as a function of stand developmental stage (age or size, denoted DS) and of permanent site conditions (S) (Equation 1):

dt dH 0 = f (DS, S) (1)
Figure 2 depicts an overall sigmoid pattern for height-age trajectories. Because they offer a good compromise between flexibility, robustness, and structural interpretation [START_REF] Zeide | Analysis of growth equations[END_REF], threeparameter sigmoid-generating differential equations are usually used for expressing Equation (1).

Growth intensity and trajectory asymptotic limit are parameterized by two scale parameters and a shape parameter controls the trajectory curvature. Rewriting Equation (1) within this frame leads to:

dt dH 0 = r f (DS, K, m) (2)
where r denotes the growth rate (meters/year) and K the height asymptotic upper limit (meters), and m is the shape parameter (dimensionless). More precisely, r is the growth rate occurring at a conventional developmental stage DS r such that f (DS r ) = 1, which can be modified through a reparameterization procedure. The dependence of the increment level on r is immediate from Equation (2). We thus considered r as the primary site parameter. For further coherence, Equation ( 2) was parameterized for r to stand as the maximal growth rate (denoted R), reached at the inflexion point of growth paths. Because K may also vary according to site, the significance of a site variation of that parameter was also investigated. A correlation between r and K was tested.

If historical environmental changes have modified site conditions over time, their effects can be summarized as functions of calendar date, applied as modifiers of any site parameter in the model (potentially r and K). Because the present sampling is retrospective, late growth is not documented for the younger stand generation (Figure 2). Consequently, the asymptote parameter K cannot be reliably estimated for that generation, as well as any dependence on calendar date (the estimate of K is assumed common to each stand pair, and is estimated mainly from the older stand data). On the contrary, R is sensitive to the early stages of height growth and can be soundly estimated in both generations of stands. The effect of calendar date was thus investigated on parameter R.

The productivity change was expressed as a function of time g(t). To keep consistent with the proportional formulation of the increment to parameter R [START_REF] Garcia | Forest modellingnew developments[END_REF], g(t) was defined as a multiplier of the increment:

dt dH 0 = g(t) R f (DS, K, m) (3)
Equation ( 3) suggests that g applies in a homogeneous way, inducing the same relative change in the R value in time whatever the initial level of R (i.e. there is no site/growth-trend interaction in a multiplicative formulation).

There are two reasons why the developmental stage might be modeled as an effect of size (H 0 ) rather than age. First, if a historical enhancement of growth exists, it implies that mensurations of successive generations of stands taken at the same age differ. Because growth is not independent of size [START_REF] Shinozaki | A quantitative analysis of plant form -The pipe model theory. I. Basic analyses[END_REF][START_REF] Valentine | Tree-growth models: Derivations employing the pipe-model theory[END_REF][START_REF] West | A general model for ontogenetic growth[END_REF], the use of age as a proxy for DS would bias the comparison of increments. Second, the replacement of DS by H 0 in Equation ( 3) suggests a dynamic interpretation of the growth trend: if we consider a fixed increment starting from a fixed initial height, a growth enhancement (g > 1) implies that less time is needed to achieve the same final height for any R value. Consequently, g acts as a time flattening factor. For instance a value of 1.5 for g is equivalent to a reduction by one third of the time needed to reach a given size with an equal initial value. Such an interpretation is of particular relevance to the present study, but is not permitted when age stands for DS. The following formulation was therefore assumed:

dt dH 0 = g(t) R f (H 0 , K, m) (4)
Time t was further taken as t = (datedate b ), where date denotes the calendar date and date b is a reference base date for trend estimation, taken as 1900. R was consequently denoted R b , representing the site conditions prevailing at the base date:

dt dH 0 = g(t) R b f (H 0 , K, m) (5)
Equation ( 5) is an ordinary differential equation with separable variables. Its integration over any increment time period [t -1 , t] (assuming a closed-form solution for H 0 ) leads to:

                  t u t u du u g t H T T t H 1 ) ( ) ( ) ( 1 0 1 0 with     ) m , K , ( R ) ( 0 b 0 0 H f dH H T (6)

Functional representation

As the mathematical representation of developmental stage influences the estimation of g (Equation ( 5)), special attention was paid to the asymptotic behavior of growth equations chosen, to represent the slow convergence of tree height growth processes (Figure 2, [START_REF] Zeide | Analysis of growth equations[END_REF]. We first considered Richards equation [START_REF] Richards | A flexible growth function for empirical use[END_REF], which has exponential convergence of height in time and is widely used in forestry. We also selected two equations with less asymptotic slowing (power convergence of time), including the Hossfeld equation [START_REF] Woollons | The Hossfeld function: An alternative model for depicting stand growth and yield[END_REF][START_REF] Zeide | Analysis of growth equations[END_REF]) and the Korf equation [START_REF] Zarnovican | Korf's growth equation (in French)[END_REF]). These were further denoted as R-, H-, and K-equations. The equations were parameterized for R b to represent the maximal growth rate. The differential forms of growth equations and expressions for T and H 0 (t) are given in the Appendix (Equations ( 7) to ( 15)). Due to the lack of theory on the growth trend pattern, g was expressed in three different forms:

(i) a linear function of t (mean trend): 1 + d 1 t (16)

(ii) a quadratic function of t, to test for a putative accelerated pattern: 1 + d 1 t + d 2 t 2 (17)

(iii) a higher-order polynomial-type function of t, to detect historical variations of decennial range. The accuracy of common high-degree polynomials is restricted by their oscillatory nature and the impact of their local behavior on the entire target range [START_REF] Burden | Numerical analysis[END_REF]. We thus considered cubic splines, which are piecewise polynomial functions of degree 3 and have the property of optimum curvature [START_REF] Lange | Numerical analysis for statisticians[END_REF]).

The spline function was written as a function of t centered on t = 0, with nodes equally-spaced at nyear intervals (between 10 and 20 years for the present purpose), and was parameterized with a base time interval [0, n]:

                2 1 0 3 k 3 1 k k 3 3 2 2 1 ) 0 , k min( pm ) 0 , k max( p d d 1 ) ( k k k n t n t t d t t t g (18)
where d 1 , d 2 , d 3 , p k , and pm k are spline parameters to be estimated, collected in a parameter vector  g , and n, k 1 , and k 2 externally specify the width and number of spline intervals necessary to describe the whole period range covered. The spline function extremities over the time period covered by available increments are not constrained by any end condition.

As g varies slowly over a usual increment time step, its integration in Equation ( 6) was approximated (except in the case of Equation ( 16), where it holds strictly) as:

  1 g 1 1 θ , 2 ) (                t t t t g du u g t u t u (19)
This made the statistical adjustment of parameters associated to cubic terms easier.

STATISTICAL METHODOLOGY

Mixed-effect models as a frame

The sampling design corresponds to longitudinal data structured according to two hierarchical levels: a stand-pair level and a stand-within-pair level, next designated levels 1 and 2. The fertility parameter R b should meet with site variations (level 1) that can be enriched by an additional withinpair variation (level 2) for a control of residual site differences among stands. A simultaneous site variation in parameter K was also tested. As mentioned above, the incompleteness of the sampling design allowed for a variation of K at level 1 alone. The shape parameter m was assumed constant.

Mixed-effects models [START_REF] Davidian | Nonlinear models for repeated measurement data[END_REF] provide an efficient way of fitting longitudinal data with structured variation in some parameters. Rather than using dummy-variable parameterization, site parameters are considered as random variables for which expectancy and variance are estimated under a given distribution assumption (e.g. [START_REF] Lappi | A height prediction model with random stand and tree parameters: an alternative approach to traditional site index methods[END_REF].

Individual estimates are then obtained by a specific posterior procedure [START_REF] Davidian | Nonlinear models for repeated measurement data[END_REF]. From Equation ( 6) the model is non-linear in the parameters and was fitted using multi-level non-linear mixed-effects models, adjusted by maximum likelihood. A normal distribution for both errors and random parameters and their independence were assumed [START_REF] Lindström | Nonlinear mixed effect models for repeated measures data[END_REF].

The model was fitted on successive non-overlapping growth intervals [START_REF] Borders | Forest growth models: parameter estimation using real growth series. P 660-667 in Forest growth modeling and prediction[END_REF], and thus predicted a height at year t conditionally to height at the former year t -1 available in data. This data structure has been recurrently recommended in literature [START_REF] Borders | Forest growth models: parameter estimation using real growth series. P 660-667 in Forest growth modeling and prediction[END_REF][START_REF] Wang | Parameter estimation of base-age invariant site index models: Which data structure to use ?[END_REF]). Due to the choice of mixed-effects models, the error-in-variable problem (Wang et al. 2007, Cieszewski and[START_REF] Cieszewski | Parameter estimation of base-age invariant site index models: which data structure to use? -A discussion[END_REF] was not addressed. The present statistical modeling method can be viewed as a "varying parameter" method [START_REF] Cieszewski | Parameter estimation of base-age invariant site index models: which data structure to use? -A discussion[END_REF] applied to height increments: all successive observations of a given forest plot are used for estimating site-dependent parameters.

From a statistical point of view, the model defined by Equation ( 6) applied to successive growth intervals [t i-1 , t i ] is written as:

        i t u t u i i du , u g t H T T t H i i                  1 g b 1 0 1 0 θ m K, , R , (20)
with:

R b = R b0 + R b1 + R b2 , K = K 0 + K 1 (21) R b1 ~ N (0, ² Rb,1 ) , K 1 ~ N(0,    R b2 ~ N (0, ² Rb,2 ) ,  i ~ N (0, ²  cor (R b1 , K 1 ) =   (23)
where  g , R b0 , K 0 , m and ² are fixed parameters, R b1 and K 1 are level-1 random components of R b and K and R b2 that of R b at level 2. Parameters  Rb,1 ,  Rb,2 and  1 are random parameter standarddeviations and   is the eventual level-1 correlation between R b1 and K 1 .

The homoskedasticity assumption was furthermore relaxed. A variance function was introduced, assuming that the residual standard deviation around the increment was proportional to a power of the expected increment:

        λ 2 1 0 0 2 σ V    i i i t H t H ε (24)
where  and  are variance function parameters, and

) ( ˆ0 t H
stands for the prediction from model (20).

As compared with a fit on age-height series where structural autocorrelation in successive observations is generated (Duplat andTran-Ha 1997, Nord-Larsen 2006), the present data structure made an autocorrelation function in the error model unnecessary [START_REF] Borders | Forest growth models: parameter estimation using real growth series. P 660-667 in Forest growth modeling and prediction[END_REF]. Models were fitted using the nlme procedure of S-PLUS [START_REF] Pinheiro | Mixed-effects models in S and S-PLUS[END_REF]. Model comparisons were drawn according to the Akaike Information Criterion (AIC). Nested models were further compared using the Likelihood Ratio Khi-2 Test (LRT). Plots of residuals versus fits/variables were used to check for adequate fitting or abnormal curvatures while testing growth equations and date effect representations.

Estimation of mean growth changemodeling steps

Step 1. The accuracy of Equations ( 13) to (15) (Appendix) was compared by a fit without any date effect. A random variation of parameter R b was allowed at levels 1 and 2. Three fits were performed: (i) no random variation for parameter K and (ii) random variation for parameter K at level 1 with correlation or (iii) without correlation to R b .

Step 2. The effect of date was introduced in the successive functional forms of Equations ( 16) to ( 18), the latter with 20-year node intervals.

Step 3. A variance function was introduced into the model to take the error structure into account.

Through examination of residual plots, two forms were tested (Equation ( 24)), as either a proportional ( = 1) or power function of height increments. In the following, fits were restricted to the most accurate growth equation and the random variation structure resulting from step 2. The cubic spline function was then tested with 15-and 10-year node intervals.

Step 4. The model with the most convenient node interval length from step 2 was fitted again with each growth equation, to further check for their relative consistency and also to consider the sensitivity of date effect estimation to equation choice. A bilateral 95% confidence interval for cubic-spline estimate of date effect was computed according to the derivation presented in the Appendix.

Interaction between growth change and site

As deviations from the multiplicative formulation suggested in Equation (3) may be possible, the date/site interaction was tested by deriving an indirect individual stand-pair index for productivity changedenoted Ibased on a jackknife-like approach. It was computed as follows: for each sample, Equation ( 20) was adjusted as many times as available stand-pairs, with successive one-toone pair exclusions. The productivity level for the year 2000 without pair k was designated as P -k , while the mean productivity level over the sample was denoted P. Assuming that pair k has experienced a sharper change in time than the mean sample, P -k should be less than P. Therefore P -P -k is positive and provides a consistent definition of I, denoted I k for stand k:

I k = P -P -k (25)
The date/site interaction was investigated by testing the relationship between index I k and level-1 estimates of R b .

RESULTS

Main characteristics of sample and quality of pairing

Stand age criteria were adequately met, with a mean difference between generations of around 75 years. The ages of older stands ranged within sound limits (Table 1). Stands were mainly located on brown soils (European Commission 2005). The main environmental characteristics (Table 3) indicated a fair average nutritional status (pH of 5.0, base saturation rate of 47%, and C/N ratio of 16.2). Soil water capacity at 1m depth (SWC) reached 140mm. In addition, the sample variations in environmental conditions were not negligible, as demonstrated by coefficients of variation (15% for C/N, 20% for pH, and 30% for SWC).

Geographical distances between paired stands are given in Table 1. Paired-stands were located very close to one another (mean distance: 0.16 km). The within-pair site differences were investigated using environmental indicators. No singular pair was detected when individual anomalies were sought. Systematic between-generations differences were then assessed using two-sided paired ttests. The mean values of main indicators and p-values for t-tests are presented in Table 3. One significant difference was found (p = 0.04) for the Ellenberg indicator of basicity (Rel), higher in the younger generations (5.55/5.18). It was not confirmed on soil pH. A trend towards higher nitrogen availability in younger stands was also detected, but was not significant (Table 3).

Estimation of mean growth change

A first assessment of growth change was provided by computing a paired t-test of the within-pair difference between the total height of young stands and that of older stands taken at the same age. A systematic and very significant difference of 5.7m on average was found (p < 10 -4 , mean comparison age: 75.4 years).

The main modeling steps and parameter estimates with the final growth equation are summarized in Table 4, also including a comparison of different growth equation fits from step 2.

Step 1. The Requation initially appeared to be the least adapted, although with a slight difference from the K-and H-equations, which had similar accuracy (all equations were within a 5 AIC unit interval). The closeness of the equations' behaviors and the high estimates for K (50 to 80 m) seemed in good accordance with the low-curvature pattern observed for old-stand growth curves (Figure 2). The introduction of random variation for K at level 1 provided significantly better fits for the R-and Hequations (p = 10 -4 and 6.0 10 -4 respectively, around -10 AIC units), and to a lesser extent for the K-equation (p = 0.03), consistently with higher estimates for K. The order of magnitude for variation in K was around 10% and remained stable in the next steps. The correlation between random variations of R b and K at level 1 was found to be positive (0.5 to 0.7 depending on the growth equation fitted), but did not lead to significantly worse fits when set to zero (p = 0.3 to 0.5).

The best fit was obtained with the H-equation, closely followed by the R-and K-equations (+4 and +9 AIC units respectively, see Table 4).

Step 2. As detected by the paired t-test, residuals plotted against calendar date displayed a positive trend. An acceleration after the 1950s and growth anomalies were also noticeable. Consistently, the introduction of a linear trend in the model (Equation ( 16)) induced significantly better estimates (p = 10 -4 ). The trend amounted to a +35% in growth rate over the century (Table 4). Quadratic and 20-year node interval cubic spline effects of date greatly improved the goodness of fits (p < 10 -4 in each case, -15 and -30 AIC units respectively). The H-equation remained the best, preceding the R-and then the K-equation (+6 and +14 AIC units respectively). The level-2 random variation of R b logically flattened. The random structure was checked again and confirmed.

Step 3. A graphical analysis of residuals revealed a heteroskedasticity pattern. A variance function (Equation ( 24)) was introduced into the model, now restricted to the H-equation, and led to a significant improvement (p < 10 -4 , ). Smaller cubic-spline node intervals were finally tested. The 15-year interval turned out to be the most appropriate (p = 0.002, no additional gain with the 10-year node interval) and was retained.

Height increments were filtered out from the estimates site and developmental stage effects and were superimposed onto the estimate of g (and its 95% bilateral confidence interval) in Figure 3.

This was achieved using the transformation based on a combination of Equations ( 6) and ( 19):

                     2 1 1 1 0 0 i i i i i i t t g t t t H T t H T ˆ (26)
-5.9 6.1 5.0 6.9 5.0 4.8 5.4 5.0 4.9 < 10 -4 <10 -4 < 10 -4 5 10 -3

(0.22) 1(Rb,K) a - 0.64 - - - - - - - -  (RSE b ) (m/
Table 3 (continued).

Figure 3 revealed a positive but irregular evolution of dominant height growth. A sharp increase occurred in the second part of the century, whereas no clear pattern appeared in the first decades.

Growth anomalies were detected in spline variations and standardized increments, the latest of which centered on the 1940s and 1990s. These observations were consistent with the significant goodness of fit gained from the cubic-spline representation of the calendar date effect and the low accuracy provided by a linear effect. The increase in growth rate reached +60% around 1980, and then oscillated between +50 and +60%.

Step 4. The final model was fitted again with each growth equation. The sensitivity of the estimate of g to both the growth equation and spline node interval selection are illustrated in Figure 4. Growth equation selection quite unaffected the estimation of g, with a maximum deviation between estimates within 5% on the g scale. It appeared that g was more sensitive to the spline node interval: at low resolution, the identification and magnitude of growth anomalies were impacted. The deviation between estimates rose to 10% on the g scale. As an example, the 20-year node interval spline failed to detect the anomaly in the 1940s. Accordingly, it was significantly less accurate than the 15-year one. Overall, the order of magnitude of g and the location of accidents in time were quite robust to the node interval resolution.

Interaction between growth change and site

Any deviation from the multiplicative formulation of date/fertility effects should be absorbed by within-pair random variations of parameter R b . A preliminary test (not presented) was therefore to plot the young-stand estimates of R b against those of old stands. No deviation from the first bisecting line was detected. Index I k (Equation ( 25)) was then calculated based on a quadratic effect of date to simplify the fitting process. With respect to the magnitude of variation of g and R b , the variability of I k was found to be very small (standard deviation of 3.9% on the g scale). I k showed no clear relationship with level-1 estimates of R b (Figure 5). Their correlation was found negative and not significant (-0.26, p = 0.35 using a Pearson t-test). trees have remained dominant during their whole lifespan. Nevertheless, trees of lower social ranks may become dominant, following local thinning or mortality events [START_REF] Delvaux | Contribution à l'étude de l'éducation des peuplements. Partie 1: Acquisition de la position dominante dans les jeunes peuplements équiennes d'Epicéa[END_REF][START_REF] Harcombe | Tree life tables. Simple birth, growth and death data encapsulate life histories and ecological roles[END_REF]. When moving backward from the sampling date, retrospective height trajectories would increasingly diverge from the true one, underestimating height at earlier stages, and concurrently overestimating the growth rate [START_REF] Magnussen | Recovering time trends in dominant height from stem analysis[END_REF]Penner 1996, Raulier et al. 2003). Whereas the issue is poorly covered in the literature, we believe that this impact should be limited. First, dominant trees were defined as the thickest and not the highest trees, which prevented sampling of opportunist trees that may favor allocation to height growth (contrary to [START_REF] Magnussen | Recovering time trends in dominant height from stem analysis[END_REF]. Second, rank progressions rarefy when thinning intensity is low [START_REF] Pardé | From 1882 to 1976/80: The silviculture experiment plots of Beech in the Haye State forest (in French)[END_REF], which is reported to be true for Beech forests in northeastern France [START_REF] Polge | The influence of thinnings on the growth constraints of Beech (in French)[END_REF]). Because we compared pairs of stands sampled at different ages, further concern may be raised about whether the intensity of such bias is age-dependent and may impact the comparison of growth. Basically, the bias originates from social rank shift events, and by nature it may increase with age, whereas the opposite seems impossible. Thus, the older stands may be more affectedbut not lessthan the younger ones, which would at worst result in underestimation of the true magnitude of the growth trend. Therefore, the +50% increase reported cannot be artifactual, and rather stands as a minimum if the bias exists.

Growth modeling: functional representations and model structure

Growth equation

Special attention was paid to selection of the growth equation and its impact on the final estimate of growth evolution. As depicted in Figure 4, the latter was weak. Actually, the assessment of growth change results from the comparison of the early-growth trajectories that are common to both generations. Hence, the asymptotic behavior of growth equations mostly ensures an adequate fitting of late-growth trajectories of mature stands. Another feature of equation fits was the significant between-stand variation of the asymptote parameter K, which was further found to be uncorrelated to R b . Such variation was in accordance with trajectory patterns (Figure 2), which do not especially suggest any convergence towards an upper limit, nor a strong tightening of the set of growth curves.

Usual datasets seldom uncover the asymptotic phase of growth patterns and may practically call into question the use of asymptotic models as the most accurate type for fitting current growth curves. Alternative growth equations based on assumptions of oblique-asymptotes or power trajectories have been suggested (Duplat andTran-Ha 1997, Cieszewski 2003) and may deserve more attention.

Hypothesis of a multiplicative effect of calendar date

The hypothesis was based on: (i) the need to keep consistent with the usual representation of the effect of site as a multiplicative parameter in growth equations, (ii) the dynamical interpretation offered by such formalism. In addition, no clear deviation from that hypothesis was highlighted by the analysis of index I k and within-pair R b variations. From a causal viewpoint, this suggested that no site factor would have constrained/saturated the response of growth to environmental changes.

However, such patterns may be hardly detectable within the usual range of site conditions of the present sample.

Cubic spline representation of historical change in dominant height growth rate

The spline function provided a detailed chronology of dominant height growth rate, in good agreement with the fluctuations depicted in the observations. It allowed us to detect a medium-term acceleration pattern starting in the 1950s and repeated growth anomalies over the century. However, its endmost increase in the 1990s might be partly artifactual, due to the absence of an end condition.

The use of spline functions may also raise concerns regarding the positioning of anomalies. Figure 3 depicted rather symmetrical variations of the spline around its local minima. If we accept that growth recovery following perturbations (say for instance climatic) may be progressive, the variation should be asymmetric by nature [START_REF] Bréda | Forest tree responses to extreme drought and some biotic events: Towards a selection according to hazard tolerance?[END_REF], and possibly result in a lagged positioning of the spline function extrema. Some confirmation may be found in the spline behavior in the 1900s and 1990s, but it remains tenuous.

Intensity of changes in productivity

In Figure 3, an index of around 1.6 in 1980 means that current height growth rate level is 60% higher at that date, relative to 1900. It is also equivalent to say that only 62.5% of the time (1/1.6) would be required to reach the same dominant height in the virtual permanent conditions of year 1980 as in those of 1900, corresponding to a time contraction factor of 37.5%. However, this does not mean that young stands growing at this virtual date would be 60% higher than in 1900 (see Figure 2), which would hold only if height trajectories were linear. Actually, the shift in growth has been progressive. The average level of g over the twentieth century amounts to 1.27, which is equivalent to a contraction factor of 21%.

The univalent relationship that links dominant height to total stand production (in volume or basal area per hectare), whatever site conditions, empirically holds at a regional scale [START_REF] Assmann | The principles of forest yield study[END_REF]).

Thus, it is plausible to assume its stability under environmental changes, and it means that the contraction factor would apply straightforwardly to production levels. Yet, because the order of magnitude reported for changes is so high and has strong potential consequences for forest management, further work is required to test the stationariness of the dominant height growthproductivity relationship over time. An implication of the present work is also that classical site index curves can no longer be used for defining site index as an updated measure of site fertility.

Based on the height growth chronology, date-dependent adaptive site index-curves were thus derived for forest management applications in Beech stands [START_REF] Bontemps | Dominant height growth of Beech in Northern France: deriving site index curves with integration of long-term trends[END_REF]).

Regarding growth trends on Common Beech, the following results were found in the literature: in the Swabian Alb in Germany, [START_REF] Untheim | Has site productivity changed? A case study in the Eastern Swabian Alb, Germany[END_REF] identified an increase in dominant height of +2.9 to +3.7m at a constant age of 80 years over the period 1930-2000. A comparison of site index of successive stand generations in Denmark led to a similar estimate of around +3.7m over a one century interval [START_REF] Skovsgaard | Increasing site productivity during consecutive generations of naturally regenerated and planted Beech (Fagus sylvatica L.) in Denmark, P. 89-97[END_REF]. Both estimates were within the range of those reported (Table 1). From a dendrochronological approach in northeastern France, [START_REF] Badeau | Long-term growth trends of Fagus sylvatica L. in northeastern France. A comparison between high and low density stands[END_REF] found a growth rate increase by around +100% over the century, twice as high as our estimate, and seemingly unrealistic. Because the latter study differed from ours in both the sampling design, growth indicator and signal extraction method, no immediate explanation for this difference can be found. However, a very similar acceleration pattern was detected in the mid century following the decennial anomaly of the 1940s [START_REF] Badeau | Long-term growth trends of Fagus sylvatica L. in northeastern France. A comparison between high and low density stands[END_REF].

CONCLUSIONS

We reported an original methodology for the estimation of long-term productivity changes in evenaged stands and its application to common Beech in northeastern France. It was intended to address earlier inconsistencies [START_REF] Spiecker | Growth trends in European forests -Do we have sufficient knowledge? P. 157-169 in Causes and consequences of accelerating tree growth in Europe[END_REF] and was based on combining the paired-plots method as an accurate way to control site fertility, the use of dominant height as a proxy for productivity, and a statistical modeling approach to estimate the chronology of dominant height growth. An accelerated increase of up to 60% at the end of the twentieth century was found. This is half the estimate reported in the same area using a dendrochronology approach [START_REF] Badeau | Long-term growth trends of Fagus sylvatica L. in northeastern France. A comparison between high and low density stands[END_REF]. The average increase over the century was +27%.

Appendix

Confidence interval for cubic-spline estimate of g(t)

We rewrite g (Equation ( 18)) as: where p is the vector of spline parameters p k and f t is the vector of f k (t) which stands for t, t 2 , t 3 , and then the [max (tn k, 0)] 3 and [min (t + n k, 0)] 3 terms (Equation ( 18)). Note that g is linear in p.

The estimator for p is a maximum likelihood estimator (ML) and is thus asymptotically unbiased and normally distributed [START_REF] Lindström | Nonlinear mixed effect models for repeated measures data[END_REF], which applies to ML , ˆt  , the derived estimator of t  , and: With the variance-covariance matrix of vector p denoted p Σ ˆ, an estimate of 2 t  is given by:

t p t t t t t f Σ f f ] p [ V ˆ2   
The variance-covariance matrix estimate for fixed effects is obtained from the statistical fit of the model, from which the sub-matrix of dimension q for the p k can be extracted. Such an estimate of (2000). Hence, with n obs and r denoting the total number of observations and model parameters respectively, n p the number of stand pairs at level 1 (between-stands variation), and assuming random variation of 2 parameters (R b and K) within that level, it becomes: df = n obs -(n p + r -2) Finally a bilateral confidence interval conditional to t at level 1- is given by:

Figure 1 .

 1 Figure 1. Location of sampled stands. Each dot corresponds to a stand-pair.

Figure 2 .

 2 Figure 2. Dominant height growth curves reconstituted from stem analyses. Dashed lines: old stands; full lines: young stands. Arrow: disturbance detected on the old stand of pair 8.

  There are recurrent contributions asking which parameter(s) should stand for the effect of site. From the traditional Algebraic Difference Approach (ADA,Bailey and Clutter 1974), anamorphic (proportional) or polymorphic (single asymptote) site curves are derived, depending on which parameter, K or r respectively in the present formulation, is allowed to vary. In the more recent Generalized Algebraic Difference Approach (GADA, Cieszewski and Bailey 2000) advanced polymorphic curves are derived (with both different shapes and multiple asymptotes), based on simultaneous variations in several site-specific parameters (usually two) constrained by specific relationships. Most of the time, the selection of the site parameter structure is based on the fitting accuracy of respective formulations (e.g. Nord-Larsen 2006).

Figure 3 .

 3 Figure 3. Historical change in dominant height growth rate. Increments are centered on their median date and have been filtered out from site and developmental stage influences (see text).Intervals indicate increment duration. Filled and blank dots: young and old stand increments respectively; thick-full line: 15-yr node interval cubic-spline estimate of g(t) conventionally set to unity at year 1900; dashed line: 95% confidence interval for g estimate; thin-full line: quadratic estimate of g(t).

Table 4 .

 4 Characteristics of main models fitted with the Hossfeld growth equation. Characteristics of models fitted with the Richards and Korf equations at the end of modeling step 1 are indicated for comparison (see text). For the last model, standard-errors of parameters are provided under parentheses. a index figures denote random variation level for standard deviation/correlation estimates, b in the absence of the variance function, c number of model parameters, d log-likelihood, e Akaike Information Criterion, f p-value associated with likelihood ratio tests (LRT) between successive nested models based on the Hossfeld equation, * not provided.

Figure 4 .

 4 Figure 4. Sensitivity of date effect estimation to growth equation and cubic spline node interval. H: Hossfeld equation; K: Korf equation; R: Richards equation. The final model selected combines the H-equation and 15-yr cubic spline node interval.

Figure 5 .

 5 Figure 5. Relationship between stand-pair productivity change index (I k ) and level-1 estimates of site parameter (R b1 ).

  Integration of each equation in any time interval [t -1 , t] leads to a closed-form solution for H0:

  for the Student distribution and df stands for the number of degrees of freedom for the estimators of p.Degrees of freedom calculation for a hierarchical mixed-effects model coincides with that of classical hierarchical analysis of variance and was as indicated in Pineiro and Bates

. Table 1. Location and age of stands. a

  Mean geographic coordinates of stand pairs (ED 50 system), b mean elevation of stand pairs in meters asl, c stand age in 1998.

	Stand pair			Elevation b	Within-pair distance	Stands age c
		Forest	Location a	(m)	(km)	(yr)		Age difference
						Young	Old
	1	Haye	6° 05' E, 48° 39' N	400	0,11	72	136	64
	2	"	6° 05' E, 48° 40' N	370	0,20	66	137	71
	3	"	6° 07' E, 48° 39' N	395	0,23	58	143	85
	4	Sarrebourg	7° 00' E, 48° 44' N	325	0,10	53	109	56
	5	Hesse	7° 04' E, 48° 40' N	325	0,07	63	157	94
	6	Lemberg	7° 17' E, 49° 00' N	295	0,11	84	142	58
	7	Mouterhouse	7° 24' E, 49° 01' N	370	0,30	53	132	79
	8	Goendersberg	7° 26' E, 49° 07' N	360	0,33	47	184	137
	9	Morimond	5° 42' E, 48° 03' N	440	0,18	56	124	68
	10	La Petite Pierre	7° 18' E, 48° 51' N	330	0,08	39	122	83
	11	FC Fislis	7° 24' E, 47° 31' N	480	0,16	90	169	79
	12	Ban d'Uxegney	6° 25' E, 48° 10' N	405	0,15	75	122	47
	13	Sainte Hélène	6° 39' E, 48° 19' N	340	0,07	65	131	66
	14	Fraize	6° 24' E, 48° 21' N	365	0,12	84	153	69
	Mean			370	0.16	64.6	140.1	75.4
	Standard-			50	0.08	15	20.2	21.7
	deviation

Table 2 . Growth increment characteristics. Standard deviations are given in parentheses. Sample size Mean duration (yr) Maximum duration

 2

	Old stands	185	9.3 (4.0)	21.5
	Young stands	171	5.0 (2.1)	13.0
	Total	356	7.2 (3.9)

Table 3 . Main environmental characteristics of sample and between-generation comparisons.

 3

	pH	pH c	5.09 (0.95)	5.11	5.07	0.67
	Base saturation rate	S/CEC (%) c	53 (36)	51	55	0.41
	Ellenberg basicity index	Rel	5.32 (1.2)	5.55	5.18	0.04
	Carbon to nitrogen ratio	C/N (kg.kg -1 ) c	16.0 (3.0)	15.8	16.3	0.36
	Ellenberg nitrogen index	Nel	5.30 (0.82)	5.46	5.15	0.29
	Phosphorus concentration P (g.kg -1 ) c	0.044 (0.041)	0.039	0.050	0.79
	Soil water capacity (1 m)	SWC100 (mm)	141 (48)	148	135	0.56
	Ellenberg humidity index Fel	5.38 (0.25)	5.37	5.39	0.78

a Standard-deviations are given in parentheses, b p-value associated with the paired t-tests, c mean of indicator for the first soil layer. Environmental indicator Sample mean a Young stands mean Old stands mean p b (paired t-test)

DISCUSSION

Sampling strategy and properties of the paired-plots method

Control and variation of fertility

The control of fertility in stand pairs relied upon the in situ assessment of environmental indicators.

The quality of pairing was further assessed from soil analyses and vegetation indicators. No singular pair was detected. From vegetation surveys, significantly more acidic conditions (p = 0.04) were found in the older stands (Rel, Table 3), but this was not reflected in the soil analyses. Stand ageing in Common Beech is associated with an acidification of the upper soil layers [START_REF] Bureau | Effect of the mixture on the humic epipedon and vegetation diversity in managed beech forests (Normandy, France)[END_REF][START_REF] Godefroid | Variation in the herb species response and the humus quality across a 200-year chronosequence of beech and oak plantations in Belgium[END_REF]. This difference in the acidity status may thus result from the difference in the maturation stage between paired stands. Conversely, no other difference was identified among all indicators measured from soil analyses and vegetation surveys, which reinforced the confidence in the effective control of site conditions.

The variation of site conditions across stand-pairs was also part of the sampling strategy. The pairing constraint may inherently limit the range of fertility that can be covered. However, significant variations in nutritional and water indicators were found (Table 3). Also, the 95% level confidence interval for parameter R b (mean and random standard-deviation estimates, Table 4) ranged from 0.32 to 0.58 m.yr -1 , corresponding to a ratio of almost 2 between extremes and suggesting that a substantial range of site conditions was encompassed.

Stands history

The comparison of distinct stand generations could be affected by uncontrolled variations in time other than environmental conditions, such as those associated with stand management. In France, an overall conversion of coppice with standards (CWS) started in the 1830s and intensified in the mid nineteenth century [START_REF] Hüffel | Les méthodes de l'aménagement forestier en France[END_REF]). Yet, Beech trees do not grow from stumps [START_REF] Oswald | Natural regeneration[END_REF]) and present stands thus result from natural regeneration, irrespective of their age. In addition, any impact of possible historical changes in silviculture is limited here, as dominant height is quite insensitive to density. Beech stands have moreover been managed at traditionally high densities in France, and thinning practices have intensified only recently [START_REF] Polge | The influence of thinnings on the growth constraints of Beech (in French)[END_REF].

Bias in dominant height reconstitution

The reconstruction of dominant height trajectories from stem analyses assumes that the sampled difficult to detect from the low number of replicates of the method and the site range encompassed.

The method may be used to assess regional or species-related variations in past trends and to gain insight into the environmental factors involved through comparisons of their chronologies. Also, it could be similarly applied to other productivity-oriented indicators such as radial growth.

ACKOWLEDGEMENTS

We gratefully thank the French Ministry for Agriculture and Fisheries (MAP) as well as the French Forest Service (ONF) for providing funding and support to the present study. The first author was also funded by a PhD grant from the MAP. We also wish to thank : C. Richter, J.-L. Dupouey, B.