Martin Maxence

forest engineer, Ph. D Raymond Patricia

Ph. D Epron Daniel

Martin Barrette

Stéphane Tremblay

Daniel Dumais

Éric Saulnier

l'Environnement Specialty Forests and their Environment Identifying and preserving old-growth attributes in mixed forest stands dominated by yellow birch (Betula alleghaniensis) and balsam fir (Abies balsamea) in the context of ecosystem management in Québec

Table of contents

Table of figures

Dbh : Diameter at breast height.

DRF : Direction de la Recherche Forestière.

EFE : Écosystème Forestier Exceptionnel.

MFFP : Ministère de la Forêt, de la Faune et des Parcs.

SSAM-II : Système Sylvicole Adaptés à la forêt Mélangée, protocole de la phase II.

Introduction 1.1 The Direction of Forest Research

This study was made at the Direction de la Recherche Forestière (DRF) of the Ministère de la Forêt, de la Faune et des Parcs du Québec (MFFP). This branch aims to develop new ways of management for the Québec forest in a context of sustainable management. For this reason, the DRF drives or take part in numerous research projects The work presented in this report was made at the Service de la sylviculture et du rendement des forêts of the DRF. Most precisely in the mixed forest research team of the pole "Sylviculture et rendement des forêts naturelles".

The context of the ecosystem management

The government of Québec has recently taken the decision to apply the ecosystem management to the whole public forest (90% of Québec forest). It can be explained by three reasons : the appearance of environmental issues related with forest exploitation, an improvement of the ecological knowledge and an increase in public sensibility for the protection of our natural heritage [START_REF] Grenon | Manuel de référence pour l'aménagement écosystémique des forêts au Québec -Module 1 -Fondements et démarche de la mise en oeuvre, Québec[END_REF].

The ecosystem management can be defined as the contribution of an ecologic vision to forest management. Its goal can be resumed as the decrease of the gap between the managed forest and pre-industrial forest and the conservation of natural variability of key elements in natural forests. Its application is in general a management close to the natural dynamic (ibid.).

The yellow birch -conifer domain

Our study focuses on two ecosystems of the mixed forest named "Balsam fir -Yellow birch stand" and "Yellow birch -Balsam fir stand". The distinction is made according to the most dominating species between yellow birch and balsam fir, mostly because of the ecological conditions in the stand. The expression "Yellow birch -conifer stands" will be used in this study to represent both of them. They are the main ecosystems of the mixed forest, an ecotone of 98 600 km², 12% of the total forest superficies of Québec. It represents the transition between the broadleaved temperate forest of the south and the boreal coniferous forest of the north. Yellow birch -conifer stands can also be found in the broadleaved temperate forest but in minor proportions. At the most advanced and most stable stage of natural dynamic, yellow birch -conifer stands follow an irregular or complex structure, mainly dominated by yellow birch (Betula alleghaniensis), balsam fir (Abies balsamea) and red spruce (Picea rubens).

The dynamic is mainly created by the regular creation of small gaps. Following the site ecological characteristics, the leading species can be either yellow birch (MJ type according to the ecological classification of Québec) or balsam fir (MS type) (Ministère des Ressources Naturelles, 2013a).

Before this state, yellow birch -conifer stands following several steps. The first strictly forest state is composed of intolerant species : paper birch (Betula papyrifera), trembling aspen (Populus tremuloides) or red maple (Acer rubrum). Balsam fir comes next, progressively, and create a mixture with paper birch and red maple to a step we can consider as intermediate. Red spruce and yellow birch appears after this stage, creating the stability state (ibid.).

According to this transitory aspect, knowledge is still missing on yellow birch -conifer stands. It mostly affects the old-growth stands, for which few information are available [START_REF] Doyon | Écologie et sylviculture de la forêt mixte : Qu'avons-nous appris au cours de la dernière décennie? Colloque de transfert de connaissances[END_REF]). Yet, this kind of forest represents the main part of preindustrial forest in the mixed forest domain (ibid. ; [START_REF] Barrette | Reconstitution historique du paysage préindustriel de la région écologique des hautes collines du Bas-Saint-Maurice[END_REF][START_REF] Vaillancourt | Avis scientifique portant sur les forêts mûres et surannées[END_REF]. It is therefore not a secondary stand but a main component of the natural forest.

1.4 Old-growth and naturalness

Old-growth stands, an important ecosystem

Conservation of old-growth stands characteristics is a main issue of the ecosystem management [START_REF] Grenon | Manuel de référence pour l'aménagement écosystémique des forêts au Québec -Module 1 -Fondements et démarche de la mise en oeuvre, Québec[END_REF]. Forest exploitation in its main aspects tends to rejuvenate the stands. Focusing on the optimal harvest period, forestry stops the maturation and old-growth stages [START_REF] Rossi | Evaluer la naturalité[END_REF].

It is now well-admitted that old-growth stands have a strong importance for biodiversity at the landscape scale as well as at the stand scale [START_REF] Vaillancourt | Avis scientifique portant sur les forêts mûres et surannées[END_REF]. Thus, it is primordial to balance society needs with environmental preservation. The development of process of forest management able to protect old-growth attributes is an urgency. However, management tools allowing the realization of this objective are still rare. In order to achieve the goal previously presented, it appears necessary to develop efficient and ergonomic protocols which can meet the managers expectations.

The evaluation of naturalness, an inspiration

Naturalness is a complex notion for which it is still difficult to found a well-admitted definition [START_REF] Winter | Forest naturalness assessment as a component of biodiversity monitoring and conservation management[END_REF]. However, the one presented by [START_REF] Vallauri | Biodiversité, naturalité, humanité. Application à l'évaluation des forêts et de la qualité de la gestion[END_REF] seems to be the most relevant : "Naturalness is the expression in one place of Nature, its biodiversity, its organization, its complexity and its spontaneous and self-governing dynamic". On the other side, in Québec, naturalness is considered as the gap between managed and natural forest (Groupe d'experts sur la sylviculture intensive des plantations, 2013).

The second definition use the pre-industrial forest as a reference, without consideration for the development stage. Thus, a present forest regeneration can be compared with the pre-industrial regeneration in a naturalness study. On the contrary, the first definition emphasis the dynamic and the absence of human impact : naturalness is stronger when the development were spontaneous and autonomous on a long period. A forest free from major or anthropogenic disturbances since centuries will have an higher naturalness than a spontaneous but young forest.

Thus, a proximity exists between old-growth forest and naturalness. Many parameters used for the evaluation of this notion are also consistent for the study of ancient forest [START_REF] Rossi | Evaluer la naturalité[END_REF]. Even if it is important to remind that differences exist between these notions, methods used for naturalness estimation seems to be source of instruction for the creation of an evaluation protocol of the conservation of old-growth characteristics.

Aim of the study

Our work will mostly consist in creating a method allowing the characterization of oldgrowth elements in a stand and their evolution after forest harvests.

The studied sites following different sylvicultural treatments. We will also try to define those which respect the most the ecosystem management exigences with a satisfying conservation of old-growth elements. This project concerns specifically the yellow birch -conifer domain. Another aim will be to improve the still fragmentary knowledge on the old-growth stands of this ecologic group.

Material and methods

CPI-M1 et SSAM-II sites

General characteristics

These two sites are two DRF study zones, located between the cities of Saint Raymond and Rivière-à-Pierre (annex I). They are both natural but managed mixed stands. But they differ in their potential vegetation, due to different ecological characteristics. The stand named "SSAM-II*1 " is in a context where yellow birch -balsam fir are favoured (ecological class MJ) when the stand named "CPI-M1*" favoured balsam fir -yellow birch stands (ecological class MS). A complete presentation of these experimental designs will be done in the part 2.1.2 and 2.1.3. The past treatments were low intensity selection cuts. The last harvest has been made 20 years ago. In SSAM-II, some plots even show no cutting clues, probably because of their low accessibility.

On each site different treatments have been made on a equal number of plots, grouped into subgroups. The experimental designs CPI-M1 and SSAM-II does not share any common treatment. A simplified presentation of the scenarios applied on the two sites site can be seen on annex II. One of the main objective of those experimental designs is the conservation of oldgrowth characteristics. The hammer finish made for all the cuts took into consideration the limits that such a decision requires. Its concrete effect is a minimal preservation of dead-wood (10-15 snags per hectare in CPI-M1 protocol ; at least 10 snags with a dbh* superior to 19 cm per hectare for SSAM-II) and the conservation of a minimum rate of trees with a wildlife value like cavity trees (5 to 10 trees per hectare in CPI-M1 ; at least 6 per hectare in SSAM-II). Moreover, specific recommendations have been made in SSAM-II for the red spruce : it was unauthorized to harvest any healthy red spruce with a dbh under 35 cm.

The experimental design CPI-M1

The principal objective of this experimental design is to estimate the efficiency of the irregular shelterwood method (CPI*) for the conservation of mixed stands characteristics and composition while being economically interesting. The conservation of old-growth characteristics is an other aim of this protocol.

The site is divided between 20 plots of 70 x 70 m equally distributed between five treatments : irregular shelterwood system with continuous cover (CPI-CP) ; irregular shelterwood system extended (CPI-RL) ; shelterwood system (CPR*) ; clearcut (CPRS*) and control plots were no treatments have been made. Some of these forestry processes occur in several steps spread over twenty years (CPI-RL and CPR). At the moment this study has been realized, we only have observed the effects of the establishment cuts for two of these treatments and not those of the final cut.

The experimental design SSAM-II

This protocol studies the effect of selection cuts made on small tree groups on the dynamic and conservation of old-growth attributes in yellow birch -conifer stands.

Another aim is the study of red spruce regeneration, a declining species for which the conservation in managed forest is complex [START_REF] Dumais | Management for red spruce conservation in Québec: The importance of some physiological and ecological characteristics -A review[END_REF].

The site is divided between 20 plots of 80 x 80 m where 4 treatments (5 plots per treatment) are equally applied. Three of them consist in a selection cut by groups. The size of these gaps varies and their distribution is not homogeneous. They differ by the harvest intensity, defined by the residual basal area (strong : 13 m²/ha ; medium : 16 m²/ha ; low : 19 m²/ha). The last scenario is used as reference.

Collected data

Our study is mainly based on data already collected on our experimental designs, before the harvest and one year after. On each plot a complete inventory of living and dead trees with a dbh superior to 9 cm has been realized. A wide variety of parameters has been collected such as species, vitality (dead, declining, alive…), vertical situation of the crown, degradation step for the snags… An evaluation of coarse woody debris has also been realized for each situation except for CPI-M1 before harvest. This protocol used two perpendicular and linear 40m long transects . Each debris crossed by the transect with a diameter on this section superior to 9 cm has been noted. Then, the volume was converted with the Von Wagner method (Von Wagner, 1968) : V = π² * Σd² / 8L. d being the diameter of each debris and L the length of the transect. All additional data required by our methodology has been recorded during summer 2014. For these reasons, they will only affects the context after harvest.

EFE data

The EFE* labels forests with exceptional characteristics. The reason can be their rarity, their ancientness or the presence of rare or threatened plant species. Their status is validated after a field inventory realized by the Direction de l'Aménagement et de l'Environnement Forestier (DAEF) agents of the MFFP, previously Ministère des Ressources Naturelles (MRN). These inventories follow a strict methodology, close to what is used in our experimental designs (diameter, species, snags…). Thus, a comparison can be easily done. To allow the validation of the stand as an EFE, the inventory has to be made where the exceptional aspect is the stronger. For this reason, they do not represent the whole stand in its complete complexity. Data from 15 of these inventories have been transmitted by the DAEF. They all concern yellow birch -conifer stand of the ecological class MJ (MS class is very rare as an EFE in our ecological context) distributed on the whole domain of the western mixed forest.

Each inventory has been made on a 400 m² area.

Old-growth clues selection

Literature review on naturalness evaluation

In order to realize a methodology for the evaluation of old-growth characteristics, we studied 6 different protocols of naturalness evaluation : [START_REF] Schnitzler | Lessons from natural forests as keys for sustainable management and improvement of naturalness in managed broadleaved forests[END_REF][START_REF] De | Quantifier la valeur écologique des milieux pour intégrer la conservation de la nature dans l 'aménagement des forêts : une démarche multicritères[END_REF][START_REF] Haye | Mise en place d'un protocole d'évaluation de la naturalité des forêts gérées en irrégulier[END_REF][START_REF] Winter | Relative Quantitative Reference Approach for Naturalness Assessments of forests[END_REF][START_REF] Greg | de Greg St-Hilaire, L'enrichissement en épinette blanche à la forêt Montmorency, un traitement sylvicole intensif de restauration écosystémique dans la sapinière à bouleau blanc de l'est, présenté en dépôt initial, en août[END_REF][START_REF] Rossi | Evaluer la naturalité[END_REF]. These studies show a wide variety of clues (annex III). A selection had to be done to retain only the most relevant for our study. This choice has been made following three main axes directions :

-Simplicity : we try to obtain an ergonomic methodology, which requires a small amount of data and is easily reusable. Moreover, we are also constrained by our study limits. -Old-growth characteristics : As point 1.4.2. shows, the concept of naturalness goes further than old-growth aspects. If some clues can be pertinent for the evaluation of the first, they may not apply to the former. -Québec context : the main part of naturalness literature comes from Europe,in a context where forests are mostly recent, fragmented and where exploitation has been strong and regular. But the context is different in Québec and some parameters hold no relevance here.

Final parameters selection

Because of the conditions previously shown, we will mainly use dendrometric parameters. The advantages of these elements are the simplicity of their measurement and their common use for forest inventories. The utilization of our methodology for other study will be not constrained by a requirement in very specific skills.

Our protocol can be divided into 4 main parts (table 1). The first focuses on the diversity in tree species (diversity) and their proportion in basal area at the standlevel (succession). It allows us to locate stands in their theoretical dynamic. An oldgrowth forest should be dominated by species of the final succession stage [START_REF] Rossi | Evaluer la naturalité[END_REF]. However, mature stands of intermediate stages can also be interesting for biodiversity (Vaillancourt & al., 2009). These data also permit their identification.

Group Criterion Description

Diversity Specific richness in tree species (nb.)

Parameter obtained with a comparison of the basal area of different groups of tree species according to the yellow birch -conifer stand dynamic :

- The second part concerns the elements allowing the development of fauna, flora and fungus on trees : the microhabitats. We consider an important part of forest vertebrates and invertebrates to be depending on these microhabitats [START_REF] Blanc | Identifying suitable woodpecker nest trees using decay selection profiles in trembling aspen (Populus tremuloides)[END_REF]. Forest exploitation has often been the cause of their decrease. It could have been a

direct action -some microhabitats are considered as a defect that motivates the harvest of the tree -or an indirect one -with the stand rejuvenation or the suppression of the senescence stage [START_REF] Rossi | Evaluer la naturalité[END_REF]. Their abundance in natural forest, and especially old-growth forests is acknowledged [START_REF] Möller | Microhabitats in lowland beech forests as monitoring tool for nature conservation[END_REF][START_REF] Rossi | Evaluer la naturalité[END_REF]. We defined three microhabitat classes : cavity originally excavated by woodpeckers (Woodpecker lodges), cavity created by other means (Natural cavities) and the last class which regroups other kind of microhabitats in living trees only (Other microhabitats). A same tree can be contained in these tree class. A more detailed presentation is shown in annexe IV.

The third focus on the structure created by the living trees in the stand. It contains some general elements (basal area, density, mean and maximum diameter) but also more specific ones. Especially the stand largest diameters class density. With these parameters, we can study these declining elements in managed forest but also the structural complexity of the whole stand. An old-growth forest should contain a significant part of large diameter elements [START_REF] Rossi | Evaluer la naturalité[END_REF] but should also follow a complex structure (Vaillancourt & al., 2009).

The last part relates to the part of dead wood in the stand, focusing both on snags and on coarse woody debris. This element is very important for biodiversity [START_REF] Tremblay | Analyse de la régénération après la CPRS dans les peuplements de la sapinière à bouleau jaune de Charlevoix et du Bas-Saguenay. Bibliothèque des PDFD[END_REF][START_REF] Angers | L'enjeu écologique du bois mort -Complément au Guide pour la description des principaux enjeux écologiques dans les plans régionaux de développement intégré des ressources et du territoire, Québec, pour le ministère des Ressources naturelles et de la Faune[END_REF] and is threatened by forest exploitation [START_REF] Rossi | Evaluer la naturalité[END_REF]. Same pattern applies for both classes. First, an indicator of global composition (density for snags, volume for coarse woody debris). Then the part of two large diameter class.

Statistical analysis

Each test has been conducted with the R software, version 3.03. All results, were considered significant when p-value ≤ 0.05 and considered illustrating a trend when pvalue ≤ 0.1.

Qualitative variables have been tested with a Khi² test. For quantitative variables, the test was different according to the hypothesis we were working on. The choice made for each assumption will be presented in the following points. Some of these tests required variances equality and data normality. For each of them, a Bartlett and a Shapiro-Wilks tests have been realized to lead us to the test choice.

Comparison of CPI-M1 and SSAM-II before the harvest

We tried to establish if some significant differences existed at the initial state between the experimental designs CPI-M1 and SSAM-II and on which parameters. The objective was to estimate the relevance of a comparison between these two stands.

In order to do so, we used for each quantitative variable a non-paired Student test to verify if the data respected normality and equality of the variance, or, if not, a Wilcoxon test.

Effect of the different treatments

The experimental designs CPI-M1 and SSAM-II are respectively divided between five and four forest treatments. For each parameter, a first test has been conducted on each site before the harvest. The aim was to discover existence of some pre-existing differences due to the natural variability of forest stands. In this case, these results should have been taken into consideration during the interpretation. Then, a second test has been made to observe the effects of the different treatments on each experimental site.

The parameter with a variances equality and a data normality has been studied with an analysis of variance test (Anova). For the other variables, and according to the low number of data (20 for each site), we choose to realize the non-parametric Kruskall-Wallis test but also an Anova. If the two tests showed same results (significant or not), we considered it valid. In case of results divergence, we considered them as invalid.

For each parameter giving a significant result, we then realized a multiple comparison analysis. Its aim was to define the treatments where differences occurred.

Results

Initial state

Comparison of CPI-M1 and SSAM-II

The stands studied in the experimental designs CPI-M1 and SSAM-II are very different. Every parameter, except the diversity, show significant differences or tendency between the two sites (table 2). CPI-M1 experimental design presents a stand with an high basal area and high density (30.49 m²/ha and 1115.6 stems/ha) but with a lack of mature elements. However, 73% of the stems recorded (dbh ≥ 9 cm) have a dbh under 30 cm and these of 50 cm or more are quite inexistent. Same pattern applies for the snags. The density in CPI-M1 is superior compared to SSAM-II (256 snags/ha in the first and 107.2 in the second) but again, they are mainly composed of snags with a dbh under 30 cm. SSAM-II plots show significantly higher diameter stems for living trees as well as for snags. Only the living stems with a diameter superior or equal to 60 cm show a tendency (p-value = 0.05513). The reason is probably the natural low density of these elements which limits the test efficiency. The results of the mean and maximum diameters confirm this image of SSAM-II : a richer in mature elements experimental design.

Even if the specific richness is the same, the repartition of the different species is different. The stable state (dominated by the yellow birch -red spruce -balsam fir group) is, on average, the most commonly found on the two sites. Nevertheless, the strong proportion of red maple in CPI-M1 causes a domination of the intermediate group characterized by red maple, white birch and balsam fir in 6 plots. This group is absent of SSAM-II. This explains the significant difference at the Khi² test (p-value = 0.02068).

As a summary of these observations, it appears the stands studied on experimental designs CPI-M1 and SSAM-II are strongly different. The old-growth characteristics are significantly more present in the second one. Due to the absence of coarse woody debris data for CPI-M1 before the harvest, we cannot display results on this aspect.

Initial variability between the treatments

For all parameters, there is no difference between plots grouped by their future treatment in CPI-M1 (Annex V). The comparison between treatments after the harvest will not be influenced by a pre-existing variability. The experimental design SSAM-II has different results for two of its parameters: snags with a dbh ≥ 30 cm (p-value = 0.021) and coarse woody debris with a diameter ≥ 30 cm (p-value = 0.0634).

For the snags, this difference occurred between the control and the STR13 and STR19 treatments. The control plots contain twice more snags with a dbh ≥ 30 cm compared to these two treatments. However, there is no significant difference with the snags with a dbh ≥ 50 cm. This disparity seems to occur only for snags with a dbh between 30 and 49 cm. For the coarse woody debris, we see a significant difference between the control and the STR16 treatment. The values are strongly inferior for the second. The absence of any coarse woody debris with a diameter ≥ 30 cm explain this result. These differences always occurred over the control treatment. Thus, it will be easy to consider them in the post-harvest analysis. The lack of information about old-growth forests in the yellow birch -coniferous domain strongly appears here. Except for the data from the EFE, all information come from other ecosystems but that we can consider as close enough for a comparison. For the stand characterized by the ecological type MS (CPI-M1), we referred to data coming from the balsam fir -paper birch domain [START_REF] Desponts | Structure de sapinières aménagées et anciennes du massif des Laurentides (Québec) et diversité des plantes invasculaires[END_REF][START_REF] Desponts | The eastern boreal old-growth balsam fir forest : a distinct ecosystem[END_REF][START_REF] Roberge | Comparaison des caractéristiques des gros chicots entre une sapinière primitive et une sapinière aménagée sur la péninsule gaspésienne , Québec[END_REF][START_REF] Mccarthy | Age and Size Structure of Gap-Dynamic , Old-Growth Boreal Forest Stands in Newfoundland[END_REF][START_REF] Aakala | Trees dying standing in the northeastern boreal old-growth forests of Quebec : spatial patterns, rates, and temporal variation[END_REF][START_REF] Tremblay | Analyse de la régénération après la CPRS dans les peuplements de la sapinière à bouleau jaune de Charlevoix et du Bas-Saguenay. Bibliothèque des PDFD[END_REF]. This ecosystem is the one following the yellow birch -conifer ecosystem toward the north. For the stand characterized by the ecological type MJ (SSAM-II), we will use the work of Shwarz & al. (2001). This forest is located in the New Hampshire (USA) and is dominated by American beech (Fagus grandiflolia) but, due to an altitudinal gradient, is close to the yellow birch -balsam fir ecosystem (HBef and Bowl stands). We will also use the work of [START_REF] Fortin | Évolution de la structure diamétrale et de la composition des peuplements mixtes de sapin baumier et d'épinette rouge de la forêt primitive après une coupe à diamètre limite sur l'Aire d'observation de la rivière Ouareau[END_REF] on the pre-industrial yellow birch -balsam fir forest despite some methodological differences which limits the comparison.

Comparison of our results with the values of other old-growth forests

The comparison between data collected in different old-growth balsam fir -paper birch forests is difficult because of their strong variability (table 3). The mean diameter in some stands [START_REF] Tremblay | Analyse de la régénération après la CPRS dans les peuplements de la sapinière à bouleau jaune de Charlevoix et du Bas-Saguenay. Bibliothèque des PDFD[END_REF] is half of the others [START_REF] Desponts | Structure de sapinières aménagées et anciennes du massif des Laurentides (Québec) et diversité des plantes invasculaires[END_REF]. The density in some stands [START_REF] Desponts | The eastern boreal old-growth balsam fir forest : a distinct ecosystem[END_REF]) can be three times superior to another one [START_REF] Desponts | Structure de sapinières aménagées et anciennes du massif des Laurentides (Québec) et diversité des plantes invasculaires[END_REF]. And the ratio of snags with a dhp ≥ 20 cm between to stands can be up to fourteen [START_REF] Desponts | Structure de sapinières aménagées et anciennes du massif des Laurentides (Québec) et diversité des plantes invasculaires[END_REF][START_REF] Tremblay | Analyse de la régénération après la CPRS dans les peuplements de la sapinière à bouleau jaune de Charlevoix et du Bas-Saguenay. Bibliothèque des PDFD[END_REF]. These important differences between stands all considered as old-growth forest is a piece of information: old-growth forests dominated by balsam fir can take a wide variety of shapes. We also underline the lack of information about living stems or snags with a dhp superior to 30 cm. It is a strong limit for the study and the comparison of the most mature elements of the forest. The experimental design CPI-M1 seems to be close to some of these old-growth forests. There is a strong similarity between its basal area,density and maximum diameter and these observed by [START_REF] Mccarthy | Age and Size Structure of Gap-Dynamic , Old-Growth Boreal Forest Stands in Newfoundland[END_REF]. The proportion of stems with a dhp superior to 20 cm and 30 cm is also close to the values of [START_REF] Tremblay | Analyse de la régénération après la CPRS dans les peuplements de la sapinière à bouleau jaune de Charlevoix et du Bas-Saguenay. Bibliothèque des PDFD[END_REF] even if the mean diameter is lower. Snag density is superior in CPI-M1 to the other stands and the wide range of results for snags with a dhp ≥ 20 cm does not allowany comparison. We can only note that our results are situated between the lowest and the highest result. Yet, the density of snags with a dbh ≥ 30 cm in CPI-M1 is lower than in the stands described by [START_REF] Roberge | Comparaison des caractéristiques des gros chicots entre une sapinière primitive et une sapinière aménagée sur la péninsule gaspésienne , Québec[END_REF] and [START_REF] Tremblay | Analyse de la régénération après la CPRS dans les peuplements de la sapinière à bouleau jaune de Charlevoix et du Bas-Saguenay. Bibliothèque des PDFD[END_REF]. If we use the data collected in control plots after the harvest, density in trees containing a woodpecker lodge is close to the results of [START_REF] Roberge | Comparaison des caractéristiques des gros chicots entre une sapinière primitive et une sapinière aménagée sur la péninsule gaspésienne , Québec[END_REF]. But they only represent a 4 plots result. Even if there were no significant differences between the control plots and the others, we cannot ensure that these results represent the whole experimental design before harvest.

HBef and Bowl stands presented by [START_REF] Schwarz | Structure and composition of three northern hardwood-conifer forests with differing disturbance histories[END_REF] shows very similar results with SSAM-II (table 4) for the stems with a dbh from 10 cm to 50 cm. For the higher diameters, the density is lower in SSAM-II. We could then consider that this site has a lack in strong diameter trees. We can also see a close repartition of yellow birch and red spruce between these stands. The main difference occurs for balsam fir, quite absent in the New Hampshire study, and for European beech, almost inexistent in SSAM-II. There is also see a strong difference in composition in the work of [START_REF] Fortin | Évolution de la structure diamétrale et de la composition des peuplements mixtes de sapin baumier et d'épinette rouge de la forêt primitive après une coupe à diamètre limite sur l'Aire d'observation de la rivière Ouareau[END_REF], especially for red spruce. This specie seems to be a very important part of the preindustrial stands.

The data from the EFE are not very useful for a stand-level reflexion. According to their small measurement area, the transcription of the values in a one hectare density gives very strong results that do not correspond to the true variability of yellow birch -balsam fir stands. These forests are a patchwork of different and small units [START_REF] Doyon | Caractérisation de la structure et du dynamisme des peuplements mixtes à bouleau jaune : pour une sylviculture irrégulière proche de la nature[END_REF] but the EFE data only depicts the oldest patch. However, we can extract some information about the potential maximum values in diameter or age of the old-growth elements.

Results after harvest

Results in CPI-M1

We can observe many significant differences between treatments after the harvest (table 5). Differences apply on the following parameters : density in stems with a dbh ≥ 20 cm containing natural (p-value = 0.0154), density in stems with a dbh ≥ 20 cm containing other microhabitats (p-value = 0.0208), basal area (p-value = 2.43e-08), density in living stems (p-value = 9.63e-06), living stems with a dbh ≥ 30 cm (p-value = 0.0195) and ≥ 40 cm (p-value = 0.0883), mean diameter (p-value = 0.0004), maximum diameter (p-value = 0.000121) and snags density (p-value = 0.000435). For all the results concerning structure or dead wood, they always follow the same pattern : they either oppose control or CPRS with all of the remaining treatments.

The divergences with control are about density, basal area and cavity trees. Such results are logic for the first two parameters. They are indeed the most impacted by the harvest and even a soft treatment will have a strong impact on them. As for the last, this result is more surprising (we were expecting a gradual decrease correlated with the intensity of the harvest). This result can be explained by a difference in natural cavity repartition between plots before the harvest.

Those occurring in opposition with the CPRS impact more parameters : basal area, density, density in living stems with a dbh ≥ 30 cm, mean diameter, maximum diameter, density in living stems containing other microhabitats and density in snags. Such results are due to the extreme impact of CPRS. Table 5 Values obtained for each treatment after harvest in CPI-M1 and results of the statistical analysis. ( "." : p-value ≤ 0.1 ; " * " : p-value ≤ 0.05 ; " ** " : p-value ≤ 0.01 ; " *** " : p-value ≤ 0.001). Letters represent significant differences between the treatments. For the parameter "Succession", "St." means "Stable" and "Int." means "Intermediate".

CPRS and control are two extremes and the situation of the other treatments compared to them is an interesting result. If the difference is only due to the CPRS, other treatments then offer a good conservation of the related parameters. When results are situated between the two limits, it then means a significant impact but softer than those of a clearcut. And if no differences between any treatments are witnessed, CPRS and control included, then these elements were too rare to be significantly impacted by the harvest. This last point mainly relates to large diameter trees. It is consistent with the point 3.1.1, where the lack of mature elements was raised. The different treatments have significant impacts on parameters in the experimental design SSAM-II (table 6) such as : basal area (p-value = 0.00227), density of living stems (p-value = 0.0299), density in living stems with a dbh ≥ 30 cm (p-value = 0.0116), density in snags with a dbh ≥ 30 cm (p-value = 0.00496), density in living stems with a dbh ≥ 20 cm containing at least one natural cavity (p-value = 0.0165). Contrary to CPI-M1, only one of these differences relates to the control versus the remaining treatments : density in snags with a dbh ≥ 30 cm. It is important to remind that a significant difference between the plots designed for each treatment before harvest has already been showed in point 3.1.2. It concerned the STR13 and STR19 treatments. This difference is only new for STR16 treatment. If we take a look at the multiple comparison analysis, we can see the gap growing for STR13 : Pr(>|t|) = 0.0352 before the harvest and Pr(>|t|) = 0.00469 after. On the opposite, the results are balanced for STR19 : Pr(>|t|) = 0.0270 before the harvest and Pr(>|t|) = 0.01806 after.

Results in SSAM-II

So the impact appears different between the treatments. It is strongest for STR13 and STR16 than for STR19, even if in the last case the initial low density in this element can explain this result.

There is few significant differences regarding the microhabitats whereas means are very different between treatments. The high variability between plots of the same group explains this result. Natural cavities are the only parameter impacted by the harvest (pvalue = 0.0165). It is stronger for STR13 and STR16 treatments where they are absolutely absent. The extreme of this result is surprising and a pre-existing lack before harvest can have an important influence on this result. Some data collected before the harvest shows that habitats that can be considered as natural cavities was already rare in STR13 and STR16. But this methodology was different of the one used in our study and therefore limits the comparison.

For the remaining parameters, we observe that the strongest differences occurred between the STR13 treatment and the control. For the two other treatments, behaviour is more complex even if STR19 is mostly close to the control. Because of the graduation in intensity of the different treatments, with STR19 being the lowest and STR13 the highest, these results are logic. The absence of significant difference for the most mature elements of the stand can prove, on the other side, a satisfying conservation of them in each treatment. Their low occurrence (especially stems with a dbh ≥ 50 cm) can also have an impact on this result still being better than in CPI-M1. Our statistical tests was not efficient for rare events.

Even in the absence of significant differences, we see a strong and unexpected variability of coarse woody debris after the harvest. The debris with a diameter superior to 30 cm disappear in three treatments and the mean volume in the control is divided by two. We should not see such important variations in a stand in free evolution and on a short period. A change in transect line between the two inventories can explain this result. So we can consider the methodology used for the estimation of coarse woody debris as non-satisfying.

Discussion

An overview of old-growth yellow birch -conifer stands

Stands dominated by balsam fir

The study of the literature on old-growth forest dominated by balsam fir and the results of the CPI-M1 experimental design show the variety of forms that this ecosystem can take. The strong dynamic of the main species of the stand can explain this observation.

The regeneration of balsam fir can persist under the canopy cover for decades and up to one century. When the canopy opens, its development is strong and the competition important. But for trees that have reached the dominant stage, the decline comes quickly [START_REF] Mccarthy | Age and Size Structure of Gap-Dynamic , Old-Growth Boreal Forest Stands in Newfoundland[END_REF]. Balsam fir is mostly a short longevity specie, around 100 years [START_REF] Hébert | Gap dynamics in Balsam Fir, Abies balsamea -Yellow Birch, Betula alleghaniensis, forests of Québec[END_REF]. There are older specimens in the boreal forest (up to 260 years-old) but their longevity is due to exceptional environmental conditions or a long waiting period under the overstory [START_REF] Mccarthy | Age and Size Structure of Gap-Dynamic , Old-Growth Boreal Forest Stands in Newfoundland[END_REF]. Natural disturbances have also a significant impact on the dynamic of balsam fir, especially spruce budworm epidemics which strongly affect the balsam fir (Ministère des Ressources Naturelles, 2013a). For these reasons, stands dominated by balsam fir tend to follow a complex structure due to the regular occurrence of mortality by small groups or by individual (ibid.).

In the studied balsam fir -yellow birch stands, the stand structure is complex but balsam fir rarely grows over 40 cm of dbh and even less above 50 cm [START_REF] Burns | Sylvics of North America : 1[END_REF]. Only red spruce and yellow birch are found in these highest diameter class in our experimental design. The lack of reference for such elements in the stands described in the literature appears as a consequence of balsam fir characteristics, mainly unable to reach this class because of its high sensibility to disturbances (spruce budworm, windfall). Red spruce and yellow birch (species absent in balsam fir -paper birch stands) do not provide a high number of large stems in CPI-M1 stands, especially compared with SSAM-II. Three factors can explain this statement : the ecological conditions in CPI-M1 are not optimal for yellow birch and red spruce (Ministère des Ressources Naturelles, 2013a) ; some natural disturbances rejuvenated the stands decades ago ; species with highest diameter were harvested first by forest exploitation.

The important proportion of typical species of the intermediate development stage seems to support the hypothesis of past disturbances, but combination of these three factors can also be expected.

As a conclusion, it seems difficult to give a single representation of old-growth stands dominated by balsam fir in the yellow birch -conifer domain. The dynamic of balsam fir is the main factor explaining this difficulty. If we can agree with the natural complex structure of these stands [START_REF] Doyon | Caractérisation de la structure et du dynamisme des peuplements mixtes à bouleau jaune : pour une sylviculture irrégulière proche de la nature[END_REF]), yet their apparent rarity in very mature elements differs to our vision of old-growth stands.

Stands dominated by yellow birch

The stands of the experimental design SSAM-II are very different from those of CPI-M1. Even if we found expected complex structure in each site [START_REF] Fortin | Évolution de la structure diamétrale et de la composition des peuplements mixtes de sapin baumier et d'épinette rouge de la forêt primitive après une coupe à diamètre limite sur l'Aire d'observation de la rivière Ouareau[END_REF][START_REF] Doyon | Caractérisation de la structure et du dynamisme des peuplements mixtes à bouleau jaune : pour une sylviculture irrégulière proche de la nature[END_REF], SSAM-II is richer in old-growth and large-diameter elements. This experimental design is closer to the usual representation of old-growth forests than CPI-M1 (Bauhus & al., 2008).

The proximity of our results in SSAM-II with those of the old-growth forests in New Hampshire [START_REF] Schwarz | Structure and composition of three northern hardwood-conifer forests with differing disturbance histories[END_REF] bring the hypothesis of a stems larger than 30cm at the dbh deficit. High diameter trees are not a specificity of southern stand, as proved by the results from the EFE. A past exploitation can be the explanation of this difference.

We also have to keep in mind that the study of Schwarz & al. is not in the same ecological context, even if they are close. Thus, American beech is an important specie in New Hampshire stands and is almost absent in SSAM-II whereas it is the contrary for balsam fir. Differences in climatic conditions, warmer in the american state, can also explain these differences.

The results from the EFE gave us interesting information about the composition or the most ancient and mature patches that composed yellow birch -balsam fir stands [START_REF] Doyon | Caractérisation de la structure et du dynamisme des peuplements mixtes à bouleau jaune : pour une sylviculture irrégulière proche de la nature[END_REF]. Even if the sampling area (400 m²) was twice the mean patch area in these forests [START_REF] Hébert | Gap dynamics in Balsam Fir, Abies balsamea -Yellow Birch, Betula alleghaniensis, forests of Québec[END_REF], it does not enable a general vision of the stand in its whole complexity. Further studies are required to improve the knowledge on this ecosystem.

However, the study of [START_REF] Fortin | Évolution de la structure diamétrale et de la composition des peuplements mixtes de sapin baumier et d'épinette rouge de la forêt primitive après une coupe à diamètre limite sur l'Aire d'observation de la rivière Ouareau[END_REF] shows a different composition of the natural stands. In pre-industrial forests, the species have a distribution strongly different from that observed in SSAM-II. Especially for red spruce, which represents half of the basal area in pre-industrial forest but only 21% in SSAM-II. Even if the decrease of red spruce in Québec is known [START_REF] Dumais | Management for red spruce conservation in Québec: The importance of some physiological and ecological characteristics -A review[END_REF], we observe in the New Hampshire oldgrowth forest of [START_REF] Schwarz | Structure and composition of three northern hardwood-conifer forests with differing disturbance histories[END_REF] a distribution of this species close to SSAM-II (figure 5.2.). Thus, red spruce can potentially be a more important specie in yellow birch -conifer stands but its proportion in SSAM-II stays in the natural variability of oldgrowth forests.

Old-growth forest in the yellow birch-conifer domain

Old-growth forests can appear under different aspects in the yellow birch -conifer domain, according to the dominant species in the stand. If yellow birch is the main specie, we stay close of old-growth characteristics usually expected in temperate broadleaved forests. The structure is complex and high diameter living or dead trees are common. If balsam fir is the dominant specie then we are not in this representation of old-growth forest anymore. Its composition is strongly variable on many parameters (density, snags…) and mature elements lack. The complex structure is the only certainty about these stands. This natural behaviour is a potential explanation for the prevalence of stands of the MJ ecological class in the yellow birch -conifer domain in the studied EFE : the MS stands are not consistent with the usual representation of oldgrowth forest and are therefore rarely selected.

Yet, the lack of knowledge about old-growth stands in the yellow birch conifer domain (Doyon & Varady-Szabo, 2013) have made our work more difficult. The experimental designs SSAM-II and CPI-M1 are insufficient to provide an accurate picture of oldgrowth forests. Furthermore, the ecological differences between these sites and those studied in the literature also limit the comparison. More important studies must be made in the yellow birch -conifer domain to provide more precise clues for the ecosystem management.

Impact of the harvest on old-growth elements

Experimental design CPI-M1

Among the different management processes applied in CPI-M1, one result is obvious : the strong impact of the CPRS. As a clearcut, its intensity compared to the other treatments was predictable. Such harvests can even change, or at least severely slow down, the natural dynamic of the stand (Arcahmbault & al., 1997, Tremblay, 2009). Thus, a decrease of these practices seems necessary. The other treatments are not much different from one another, control excepted. The low density of trees with natural cavities seems to be the only common effect between the four harvest methods. Again, the lack of previous data limits the interpretation for this parameter.

Two points may explain these results. First, the initial density in mature and large elements was low in the majority of plots. Then, these densities were close to the minimal conservation standard for old-growth elements used during the hammer finish (Déry & Leblanc, 2005). The harvest of these elements was then strongly limited for each treatment, except for the CPRS. For the stems with a dbh ≥ 50 cm, the comparison with clearcut does not show any result because of their initial low density.

It is also necessary to take into account the next harvests. For two treatments, CPI-RL and CPR, it was only an opening of the canopy for the development of the regeneration.

In the next decades, the remaining will also be harvested. We can therefore expect an impact close to the CPRS for the old-growth elements. However, for the CPI-CP treatment, the exact same harvest will be repeated around 20 years. The aim of this scenario is to always keep at least 40% of the canopy. This treatment is expected to enable a great conservation of old-growth elements [START_REF] Raymond | The Irregular Shelterwood System : Review, Classification, and Potential Application to Forests Affected by Partial Disturbances[END_REF][START_REF] Doyon | Écologie et sylviculture de la forêt mixte : Qu'avons-nous appris au cours de la dernière décennie? Colloque de transfert de connaissances[END_REF].

Moreover, the repetition of harvest aiming to keep a complex structure can also lead to an standardisation and a loss in natural elements of these stands [START_REF] Angers | Comparing composition and structure in old-growth and harvested (selection and diameter-limit cuts) northern hardwood stands in Quebec[END_REF].

The strong decrease of living trees with other microhabitats than cavities and with a dbh ≥ 20 cm seems to be a clue of this threat. These stems are those with the highest probability to develop some cavities or lodges in the future. Their regular harvest has a high probability to lead to a decrease in the number of the microhabitats in the stand, and so to impact the biodiversity. So long-term survey are necessary to corroborate this hypothesis.

Experimental design SSAM-II

Most of the differences we can observe between treatments are predictable : the impact will be much stronger as the harvest will be intense. Yet, all the treatments tend to preserve the complex structure and the old-growth elements. Stems with a dbh ≥ 40 cm do not show important variations. The snags with a dbh ≥ 30 cm are the only parameter where treatments seems to have a strong impact. But in fact, they only enlarge a preexisting gap between control and the other treatments. However, if we are only looking at the intensity of the harvest, we can notice that its effect on this element is stronger in STR16 whereas we expected this result for the treatment STR13. Maybe the standard for the conservation of dead wood has limited the impact for this last scenario. However, it seems difficult to give an interpretation for the microhabitats. We consider that the lack of any natural cavities in the treatments STR13 and STR16 cannot be explained only by the harvest. An initial lack in these elements is very likely. The same logic is true for woodpecker lodges in STR19. Their low density does not seem consistent with the low impact of this scenario. Nevertheless, the results for living trees with microhabitat follow the intensity of the harvest : their density decrease when this parameter increase.

In general, these treatments confirm the confidence in selection cuts for the conservation of old-growth characteristics (Colloque sur les vieilles forêts boréales, 2009 ; [START_REF] Doyon | Écologie et sylviculture de la forêt mixte : Qu'avons-nous appris au cours de la dernière décennie? Colloque de transfert de connaissances[END_REF]. However, with the intensity of the harvest increases the minimal period between two cuts for a good regeneration of the stand. A selection cut remains very different from the natural mortality of a stand [START_REF] Doyon | Caractérisation de la structure et du dynamisme des peuplements mixtes à bouleau jaune : pour une sylviculture irrégulière proche de la nature[END_REF] and does not have the same impact than individual mortality or regeneration by small gaps [START_REF] Doyon | Enjeux écologiques de la forêt feuillue tempérée québécoise, Québec, pour le ministère des Ressources naturelles et de la Faune[END_REF]. A low-intense and frequent treatment seems to be the best way to maintain a dynamic close to the natural one while keeping some plasticity in the management. A soft harvest allows a more precise choice among the stems so a better conservation of oldgrowth elements. Again, a long-term survey is necessary to be sure that the repetition of the harvest does not lead to a standardisation or a modification of the stand [START_REF] Angers | Comparing composition and structure in old-growth and harvested (selection and diameter-limit cuts) northern hardwood stands in Quebec[END_REF].

Comparison of these results with old-growth standards

The interpretation of the standard submitted by Déry & Leblanc ( 2005) is relatively free. They do not show any precise diameter categories, preferring the use of the word "large". Thus, they advise maintaining (per hectare) 10 to 15 "large" snags, 5 to 10 "large" living stems, with microhabitats or not, and at least 5 m³ of "large" coarse woody debris. The use of this word instead of a specified threshold seems to be an important limit.

If we consider as "large" a dbh of 40 cm for living trees (Ministère des Ressources Naturelles, 2013b) and 30 cm for a snag [START_REF] Tremblay | Analyse de la régénération après la CPRS dans les peuplements de la sapinière à bouleau jaune de Charlevoix et du Bas-Saguenay. Bibliothèque des PDFD[END_REF], we note the compliance with this standard in CPI-M1 as shown in point 3.2.1. However the decrease of these thresholds have a strong impact: a living tree with a dbh of 32 cm can be considered as large [START_REF] Fortin | Évolution de la structure diamétrale et de la composition des peuplements mixtes de sapin baumier et d'épinette rouge de la forêt primitive après une coupe à diamètre limite sur l'Aire d'observation de la rivière Ouareau[END_REF] and a snag can be considered large at 20 cm [START_REF] Roberge | Comparaison des caractéristiques des gros chicots entre une sapinière primitive et une sapinière aménagée sur la péninsule gaspésienne , Québec[END_REF]. It would allow an intensification of the harvest without a decrease of old-growth characteristics. Referring ourselves to point 3.2.1 results we can observe this situation to be non-sensical. Last, this standard is not adapted for SSAM-II.

Even with the largest diameter categories, every threshold is substantially exceeded.

The one from the Ontario Tree Marking Guide [START_REF] Omnr | Ontario Tree Marking Guide, Version 1.1[END_REF]) is strongly different. It advises to keep at least 12 cavity trees per hectare, half with a dbh ≥ 40 cm, the others ≥ 20 cm. Among these stems, 6 should contain woodpecker lodge with at least one with a dbh ≥ 40 cm. Every 4 hectares one supercanopy tree with a dbh ≥ 60 cm should be maintained. The standard also advises not to harvest until the basal area exceeds 28 m²/ha and to keep at least a residual basal area of 20 m²/ha (including 3 m²/ha of trees with a dbh between 50 and 60 cm and 2 m²/ha of trees above 60 cm).

About microhabitats, the results in the control plots of SSAM-II are above the thresholds. CPI-M1 controls plots are also close to them, except a lack in cavity trees with a dbh ≥ 40 cm, certainly because of their low density in these stands. The advices for the harvest are close to our previous conclusions but the basal area thresholds for trees with a dbh above 50 cm seems to be inadequate with our stand, especially for CPI-M1 and to a lower extent for SSAM-II. Its values are half of those recommended by the OMNR. A difference of productivity between our stands and those used for the definition of these threshold can explain this result. Moreover, the advice of one supercanopy tree can be relevant in balsam fir -yellow birch stands but not in stands where yellow birch is the main specie.

As a conclusion, the Déry & Leblanc's standard can be an interesting tools for stands dominated by balsam fir but its lack in clear thresholds is an important limit. The standard developed by the OMNR is very interesting and its advice for microhabitats could be used in mixed forest but the thresholds related with the basal area of the stand do not appear adapted to this ecosystem. It underlines again the necessity to improve our knowledge about old-growth stands in the yellow birch conifer domain for a development of relevant thresholds.

Critic of our methodology

Our methodology is mostly based on data usually collected during forest inventories. For this reason, its use and repetition seems very simple even if some differences between CPI-M1 and SSAM-II results show the limits of some parameters. For example the high diameter trees in CPI-M1 are too rare to have a significant impact on statistical analysis but still are important for a right vision of the stand. The transformation of such data into a presence/absence figure can be a solution.

Applying this methodology in other stands of the yellow birch -conifer domain is necessary to improve our knowledge on this ecosystem. It will lead to the development of relevant old-growth threshold but will however show the limits of our method and then will permit its amelioration (for example a new protocol for coarse woody debris evaluation with the Von Wagner method). The high variability of its results is probably due to the use of different transects in the two collect sessions. Thus, even if it is a timeefficient method, it seems better to choose a more precise protocol.

The absence of some parameters usually linked with the notion of old-growth forests can also be seen as an error. As an example an analysis of the vertical structure [START_REF] Bauhus | Silviculture for old-growth attributes[END_REF] or the availability in dead wood depending of the rotting classes [START_REF] Winter | Relative Quantitative Reference Approach for Naturalness Assessments of forests[END_REF]. But this lack can be explained by different reasons. No data about the vertical structure has been collected in our experimental designs and they require a complex methodology [START_REF] Haye | Mise en place d'un protocole d'évaluation de la naturalité des forêts gérées en irrégulier[END_REF]. This is the opposite of our objective to offer a simple methodology. Moreover, we consider that the presentation of the different diameter classes allows the evaluation of the stand horizontal complexity but also, indirectly, its vertical complexity. The use of rotting class was an interesting parameter but a clear presentation of its results was complex. Furthermore, we consider that a better knowledge about deadwood in the yellow birch -conifer domain is necessary to use as a relevant parameter.

Improving the ecological knowledge on the old-growth mixed forest of yellow birch and balsam fir is another important aim. A large part of forest biodiversity depends on oldgrowth characteristics (Colloques sur les vieilles forêts boréales, 2009). That is why we have to certify the efficiency of these elements protection in managed forests. Moreover, some methodologies for naturalness evaluation use wildlife indicator and are therefore very important to understand in the connexion between old-growth forests and biodiversity [START_REF] Greg | de Greg St-Hilaire, L'enrichissement en épinette blanche à la forêt Montmorency, un traitement sylvicole intensif de restauration écosystémique dans la sapinière à bouleau blanc de l'est, présenté en dépôt initial, en août[END_REF]). Yet, methodologies using dendrometric parameters appear more efficient due to their capacity to be used by most forest workers and not only by specialists.

Conclusion

This study highlights the complexity of defining old-growth yellow birch and balsam fir forests. Their characteristics can be very different according to the dominant species of the stand. Our work showed the lack of knowledge about these ecosystems yet they represent an important part of pre-industrial mixed forests. This can be interpreted as a limit for a correct application of the ecosystem management in mixed forest.

Our methodology seems relevant to enhance this knowledge while focusing on data commonly collected by forest workers. Moreover, a largest utilization is necessary to identify its limits and so being able to improve it. This study also underlines the efficiency of continuous-cover treatments such as selection cuts or irregular shelterwood systems to maintain old-growth attributes. However, long-term analysis are necessary to ensure the real effect of such treatments on stands. Moreover, the current standards for conservation of old-growth elements are not satisfying in a mixed forest context. Further studies are necessary to improve our knowledge about old-growth mixed forests and about the best ways to maintain these elements in managed stands. It will provide essential information for the application of the ecosystem management in this ecosystem.

Annexes

Regular shelterwood method (CPR) :

The CPR is a regeneration process consisting in an harvest in several steps. It starts with a partial opening of the canopy. Then, other partial harvest can be made until the final harvest, where all the remaining stand is cut. The duration of this process does not exceed 1 / 5 of stand revolution.

Irregular shelterwood method (CPI) :

The CPI follows the same steps than the CPR but in a period exceeding the fifth of the revolution. The final harvest is also optional. But three different CPI types exist and two were applied in our study :

-Continuous cover (CPI-CP) : its aim is to maintain a continuous cover in the stand. Thus, a harvest shall keep at least 40% of the commercial stems. -Extended (CPI-RL) : This process is very close to the CPR. The main difference is the longer period between the opening of the canopy and the final harvest (more than 1 / 5 of the revolution)

Clearcut (CPRS) :

In a clearcut, all the commercial stems of the stand are harvested. But forest machines can only use predefined skid trails to limits the damages on soil and regeneration.

Selection cut (CJPG) :

A selection cut regroup three different aims in one cut : harvest, amelioration and regeneration. For this reason, a selection harvest just applies to a part of the stand. In our context, it was made by gap creations. Their repartition and their size where heterogeneous (maximum diameter : 25m

Summary

The conservation of old-growth elements in managed forest is an important aim of the ecosystem management in Québec. This work focus on their representation in yellow birch (Betula alleghaniensis) and balsam fir (Abies balsamea) stands of the mixed forest domain. We studied two experimental sites, one dominated by yellow birch and the other by balsam fir, where was estimated the effect of different forest treatments on oldgrowth characteristics. Our method shows the possibility to represent them with a small number of parameters, mainly dendrometric. Continuous cover treatments seems to be the most interesting systems for the conservation of these elements. However, our study raised the lack of information concerning old-growth stands in mixed forest. Yet, these data are essential for a relevant application of the ecosystem management. For this reason, further researches on this subject are necessary.

La conservation d'éléments propres aux forêts anciennes dans les peuplement aménagés est un enjeux important de l'aménagement écosystémique au Québec. Cette s'étude s'est principalement intéressée à leur expression dans les peuplements de bouleau jaune (Betula alleghaniensis) et de sapin baumier (Abies balsamea) dans le domaine de la forêt mixte. Pour cela, nous avons travaillé sur deux sites d'expérimentation, l'un dominé par le bouleau jaune et le second par le sapin baumier, visant à estimer l'effet de différents traitements forestiers sur les caractéristiques de vieilles forêts. Notre méthode montre qu'il est possible de les représenter en se basant sur un nombre réduits de paramètres, essentiellement dendrométriques. Les itinéraires de gestion sylvicole à couvert continu apparaissent alors comme les systèmes les plus adaptés à leur conservation. Notre étude a cependant mis en évidence le manque d'informations disponible sur les vieux peuplements du domaine de la forêt mixte. Ces données sont pourtant essentielles pour la réalisation concrète de l'aménagement écosystémique. De nouvelles recherches sur ces écosystèmes sont donc nécessaires.

Annex I :

 : Cartography of the experimental designs CPI-M1 et SSAM-II Localisation of the two experimental designs (source : Google.ca) Experimental design SSAM-II (source : C. Bourdon) of the different treatments in CPI-M1 et SSAM-II (source : Ministère des Ressources Naturelles, 2013b)

Table 1

 1 Presentation of the parameters studied in the old-growth evaluation. ................ 6

Table 2

 2 Mean values obtained in the experimental designs SSAM-II and CPI-M1 before the harvest and results of the statistical analysis. .................................................... 9

Table 3

 3 Bibliographical synthesis of the characteristics of old-growth stands dominated by balsam fir, compared with the results of CPI-M1 before harvest. ...................... 10

Table 4

 4 Comparison of the results of SSAM-II before harvest with data from the EFE and with the work of Schwarz & al. (2001). ............................................................ 11Table5Values obtained for each treatment after harvest in CPI-M1 and results of the statistical analysis. . .

.............................................................................................. 13

Table 6 Values obtained for each treatment after the harvest in SSAM-II and results of the statistical analysis. . ......................................................................................... 14 List of appendices Annex I Cartography of the experimental designs CPI-M1 et SSAM-II………..………….I Annex II Presentation of the different treatments in CPI-M1 et SSAM-II……………..…..II Annex III Naturalness clues raised in the literature review……………………………..….III Annex IV Microhabitats categories…………………………………………………………....V Annex V Results of the different studied parameters in the experimental designs CPI-M1 and SSAM-II before the harvest, grouped by treatments……………………………..VI Acronyms Note : The French words or expressions are shown in italics CPI : Coupe Progressive Irrégulière CPR : Coupe Progressive Régulière. CPRS : Coupe avec Protection de la Régénération et des Sols. CPI-M1 : Coupe Progressive Irrégulière en peuplement Mixte n˚1.

Table 1 Presentation of the parameters studied in the old-growth evaluation.

 1

	Dynamic	Succession	Young (Yo.) : paper birch, trembling aspen and red maple -Stable (St.) : balsam fir, red spruce and yellow birch -Intermediate (Int.) : balsam fir, red maple and paper birch
			-Companion (Comp.) : other species (sugar maple, american beech...)
			The group with the highest basal area is considered as the dominant succession stage.
		Woodpecker lodge	Number of trees, alive or dead, containing at least one woodpecker lodge (nb./ha)
	Microhabitats	Natural cavity	Number of trees, alive or dead, containing at least one natural cavity (nb./ha)
		Other	Number of living trees only with other microhabitats than the one previously defined
		microhabitats	(nb./ha)
		Density	Number of living stems with a dbh ≥ 9 cm (nb./ha)
		Basal area	Calculated basal area of the stand (m²/ha)
	Structure
	Dead wood	Volume of coarse	Volume in coarse woody debris with a diameter ≥ 9 cm at the transect intersection
		woody debris	(m³/ha)
		Debris ø ≥ 30 cm	Part of the coarse woody debris volume with a diamter ≥ 30 cm (%)
		Débris ø ≥ 50 cm	Part of the coarse woody debris volume with a diamter ≥ 50 cm (%)

Mean diameter

Mean diameter at breast height of the living trees (cm) Maximum diameter Highest diameter at breast height in the stand (cm) Stems dbh ≥ 30 cm Number of living stems with a dbh ≥ 30 cm (nb./ha) Stems dbh ≥ 40 cm Number of living stems with a dbh ≥ 40 cm (nb./ha) Stems dbh ≥ 50 cm Number of living stems with a dbh ≥ 50 cm (nb./ha) Stems dbh ≥ 60 cm Number of living stems with a dbh ≥ 60 cm (nb./ha) Snags Number of snags with a dbh ≥ 9 cm (nb./ha) Snags dbh ≥ 30 cm Number of snags with a dbh ≥ 30 cm (nb./ha) Snags dbh ≥ 50 cm Number of snags with a dbh ≥ 50 cm (nb./ha)

Table 2 Mean values obtained in the experimental designs SSAM-II and CPI-M1 before the harvest and results of the statistical analysis

 2

		Site	SSAM-II	CPI-M1
		Variable	Mean	Mean
			Stand.-dev.	Stand.-dev.	p.value
	Dynamic	Diversity	6	6	1 n.s.
		(nb. species)	0,9	1,1
		Succession	Stable	Stable	0.02068 *
		Stable (nb. plots)	19	14
		Intermediate (nb. plots)	0	6
		Compaign (nb. plots)	1	0
	Structure	Maximum age	180	113	3.462e-06 ***
		(years)	46,3	38,7
		Basal area	25,25	30,49	2.847e-05 ***
		(m²/ha)	3,5	3,5
		Mean diameter	21	18	0.001231 **
		(cm)	3,5	1,3
		Maximum diameter	56	48	0.007511 **
		(cm)	8,1	10,5
		Density	536,6	1115,6	9.823e-09 ***
		(nb./ha)	172,7	287,0
		Stems dbh ≥ 30 cm	118,1	68	2.218e-06 **
		(nb./ha)	26,8	29,9
		Stems dbh ≥ 40 cm	40,6	13,4	4.236e-06 ***
		(nb./ha)	15,9	12,1
		Stems dbh ≥ 50cm	6,9	1,6	0.0005829 ***
		(nb./ha)	5,7	4,0
		Stems dbh ≥ 60 cm	2,2	0,4	0.07755 .
		(nb./ha)	3,7	1,2
	Deadwood	Snags	107,2	256	2.015e-10 ***
		(nb./ha)	53,4	54,9
		Snags dbh ≥ 30 cm	29,7	17,4	0.01405 *
		(nb./ha)	12,1	16,5
		Snags dbh ≥ 50 cm	5,0	1,4	0.007188 **
		(nb./ha)	3,3	4,8

( "." : p-value ≤ 0.1 ; " * " : p-value ≤ 0.05 ; " ** " : p-value ≤ 0.01 ; " *** " : p-value ≤ 0.001).

Table 3 Bibliographical synthesis of the characteristics of old-growth stands dominated by balsam fir, compared with the results of CPI-M1 before harvest.

 3

Table 4 Comparison of the results of SSAM-II before harvest with data from the EFE and with the work of Schwarz & al. (2001).

 4

Table 6 Values obtained for each treatment after the harvest in SSAM-II and results of the

 6

	statistical analysis. ( "." : p-value ≤ 0.1 ; " * " : p-value ≤ 0.05 ; " ** " : p-value ≤ 0.01 ; " *** " : p-value ≤
	0.001). Letters represent significant differences between the treatments. For the parameter "Succession",
	"St." mean "Stable".

Results of the different studied parameters in the experimental designs CPI-M1 and SSAM-II before the harvest, grouped by treatments.

		Treatment			STR13					STR16			STR19	TEM
		Variables			Mean					Mean			Mean	Mean
				Stand.-dev.	Stand.-dev.	Stand.-dev.	Stand.-dev.	p.value
	Dynamic	Diversity		6,2					6,4				6,0	6,8	0.549	n.s.
		Treatment (nb. species)		CPI_25 0,4	CPI_40 0,9	CPR_50	1,2	CPRS	0,8	Témoin
		Variables Succession		Mean St.		Mean St.	Mean	St.	Mean	St.	Mean	0.3679 n.s.
		Stable (nb. plots)	Stand.-dev. Stand.-dev. Stand.-dev. Stand.-dev. Stand.dev 4 5 5 5	p.value
	Dynamic	Diversity Intermediate (nb. plots)	6	1		6			0		7		0	6,8	0	6	0.612	n.s.
	Structure	(nb. species) Basal area		1,4 25,78		0,8 25,20		1,2		24,83	1,0	25,21	1,4	0.983	n.s.
		Succession (m²/ha)		St. 3,07		St.	3,50		St.		2,70	St.	5,23	St.	0.9169 n.s.
		Stable (nb. plots) Mean diameter		3 21,6		3		20,6		3		20,6	2	20,2	3	0.925	n.s.
		Intermediate (nb. plots) (cm)		1 4,3		1			3,1		1		3,0	2	4,4	1
	Structure	Basal area Maximum diameter	31,42 56,4		29,59 54,7	30,41	53,0	29,19	61,1	31,83	0.83 0.471	n.s. n.s.
		(m²/ha) (cm)		4,85 9,3		2,12 10,9		4,13		5,3	3,54	6,0	3,62
		Mean diameter Density		17,3 540		18,1	559		17,0		483	19,0	565	16,7	0.313 0.89	n.s. n.s.
		(cm) (nb./ha)		1,7 266		1,6	143		0,9		75	2,8	199	0,7
		Maximum diameter Stems dbh ≥ 30 cm	44,7 118,8		54,4 108,8		44,5		133,8	51,6	111,3	44,3	0.537 0.482	n.s. n.s.
		(cm) (nb./ha)		2,8 41,7		16,1 25,2		6,8		13,7	15,9	19,5	3,8
		Density Stems dbh ≥ 40 cm	1210 43,8		1007 38,8		1184		41,3	905	38,8	1272	0.356 0.96	n.s. n.s.
		(nb./ha) (nb./ha)		380 21,2		247 20,4		288		16,3	243	6,8	220
		Stems dbh ≥ 30 cm Stems dbh ≥ 50cm	63,0 8,8		73,0	5,0		58,0		6,3	87,0	7,5	59,0	0.664 0.782	n.s. n.s.
		(nb./ha) (nb./ha)		30,4 7,1		34,6	5,2		19,2		7,7	34,5	2,8	34,6
		Stems dbh ≥ 40 cm Stems dbh ≥ 60 cm	11,0 2,5		13,0	2,5		7,0		1,3	25,0	2,5	11,0	0.29 0.8944 n.s. n.s.
		(nb./ha) (nb./ha)		6,0 3,4		8,9		5,6		11,5		2,8	17,4	3,4	11,5
	Deadwood	Stems dbh ≥ 50cm Snags		0,0 91,3		2,0 137,5		2,0		73,8	4,0	126,3	0,0	0.632 0.202	n.s. n.s.
		(nb./ha) (nb./ha)		0,0 51,5		4,0	65,6		2,3		42,0	8,0	38,6	0,0
		Stems dbh ≥ 60 cm Snags dbh ≥ 30 cm	0,0 21,3 a		1,0	28,8 a	0,0		1,0 22,5 a	0,0 46,3 a	0.573 0.045	*	n.s.
		(nb./ha) (nb./ha)		0,0 7,1		2,0	18,5		0,0		11,4	2,0	16,3	0,0
	Deadwood	Snags Snags dbh ≥ 50 cm	279,0 7,5		216,0 2,5	280,0	5,0	263,0	5,0	242,0	0.479 0.468	n.s. n.s.
		(nb./ha) (nb./ha)		41,4 5,2		62,4	3,4		44,9		5,2	69,8	5,2	59,9
		Snags dbh ≥ 30 cm Coarse woody debris	17 40,67		20	35,38		12		41,58	27	54,97	11	0.352 0.58	n.s. n.s.
		(nb./ha) (m³/ha)		10 29,31		20 19,18		9		18,14	12	22,44	7
		Snags dbh ≥ 50 cm (nb./ha) Debris ø ≥ 30 cm (%)		3,0 6,0 0,16 ab 0,17	0,0 0,0	0,07 a 0,10	0,0 0,0		3,0 3,8 0,29 b 0,09	1,0 0,31 b 2,0 0,07	0.52 0.0158 *	n.s.
	Debris ø ≥ 50 cm Results in CPI-M1 at the initial state 0,00 (%) 0,00				0,00 0,00				0,00 0,00	0,00 0,00	n.a.	n.s
	Results in SSAM-II at the initial state
											VI VII

). III Annex V :

All the terms followed by the symbol * are explained in the presentation of the acronyms.


Luca Sanguinetti

Antonio A D'amico

Michele Morelli

Merouane Debbah

Random Access in Uplink Massive MIMO Systems: How to exploit asynchronicity and excess antennas

Massive MIMO systems, where base stations are equipped with hundreds of antennas, are an attractive way to handle the rapid growth of data traffic. As the number of users increases, the initial access and handover in contemporary networks will be flooded by user collisions. In this work, we propose a random access procedure that resolves collisions and also performs timing, channel, and power estimation by simply utilizing the large number of antennas envisioned in massive MIMO systems and the inherent timing misalignments of uplink signals during network access and handover. Numerical results are used to validate the performance of the proposed solution under different settings. It turns out that the proposed solution can detect the collisions with a probability higher than 90%, while providing reliable timing and channel estimates at the same time. Moreover, numerical results demonstrate that it is robust to overloaded situations.

I. INTRODUCTION

Massive MIMO is considered as one of the most promising solution to handle the dramatic increase of mobile data traffic in the years to come [START_REF] Marzetta | Noncooperative cellular wireless with unlimited numbers of base station antennas[END_REF], [START_REF] Larsson | Massive MIMO for next generation wireless systems[END_REF]. The basic premise behind massive MIMO is to reap all the benefits of conventional MIMO, but on a much greater scale: a few hundred antennas are used at the base station (BS) to simultaneously serve many tens of user equipment terminals (UEs) in the same frequencytime resource using a time division duplexing protocol. It is worth observing that, contrary to what the name "massive" suggests, massive MIMO arrays are rather compact; 160 dualpolarized antennas at 3.7 GHz fit into the form factor of a flat-screen television [START_REF]A flexible 100-antenna testbed for massive mimo[END_REF]. The benefits of massive MIMO in terms of area throughput, power consumption and energy efficiency have been extensively studied in recent years and are nowadays well understood [START_REF] Larsson | Massive MIMO for next generation wireless systems[END_REF], [START_REF] Marzetta | Massive MIMO: An introduction[END_REF], [START_REF] Bjornson | Massive MIMO: ten myths and one critical question[END_REF]. On the other hand, the potential benefits of massive MIMO in the network access functionalities have received little attention so far [START_REF] Sorensen | Massive mimo for crowd scenarios: A solution based on random access[END_REF]- [START_REF] Björnson | Random access protocol for massive MIMO: strongest-user collision resolution (SUCR)[END_REF].

Generally speaking, the network access functionalities refer to all the functions that a UE needs to go through in order to establish a radio link with the BS for data transmission and reception. To this end, the LTE standards specify a network entry procedure called random access (RA) by which uplink signals can arrive at the BS aligned in time and with approximately the same power level [START_REF] Sesia | LTE -the UMTS long term evolution: from theory to practice[END_REF]. In its basic form, the RA function is a contention-based procedure, which essentially develops through the same steps specified by the initial ranging process of the IEEE 802.16 family of standards.

This research has been supported by the ERC Starting Grant 305123 MORE, and by the research project 5GIOTTO funded by the University of Pisa.

Specifically, each UE trying to establish a communication link first acquires synchronization in DL and then accesses the network using the so-called random access channel (RACH), which is composed of a specified set of consecutive symbols and adjacent subcarriers. Each UE makes use of the RACH to transmit a randomly chosen code. As a consequence of the different UEs' positions within the cell, uplink signals are subject to UEs' specific propagation delays and arrive at the BS at different time instants. After detecting the selected codes and identifying how many UEs are using a given code (collision resolution), the BS proceeds extracting the corresponding timing and power information. Then, it will broadcast a response message informing the detected UEs that the random access procedure has been successfully completed and giving instructions for timing and power adjustment. The undetected UEs repeat the random access procedure until success notification. From the above discussion, it follows that code identification as well as multiuser timing and power estimation are the main tasks of the BS during the RA process.

The above procedure can be used for various purposes: initial access, handover, maintaining uplink synchronization, and scheduling request. Among the various purposes of the RA procedure and differently from [START_REF] Sorensen | Massive mimo for crowd scenarios: A solution based on random access[END_REF]- [START_REF] Björnson | Random access protocol for massive MIMO: strongest-user collision resolution (SUCR)[END_REF], we focus on initial access and handover in massive MIMO networks. Both problems have received great attention in conventional MIMO systems and different solutions are currently available [START_REF] Sesia | LTE -the UMTS long term evolution: from theory to practice[END_REF]- [START_REF] Miao | A ranging method for OFDMA uplink system[END_REF]. The methods illustrated in [START_REF] Zhou | OFDMA initial ranging for IEEE 802.16e based on time-domain and frequency-domain approaches[END_REF] performs code detection and timing recovery by correlating the received samples with time-shifted versions of a training sequence. A simple energy detector is employed in [START_REF] Mahmoud | Initial ranging for WiMAX (802.16e) OFDMA[END_REF] whereas a timing recovery scheme specifically devised for the LTE uplink is discussed in [START_REF] Sesia | LTE -the UMTS long term evolution: from theory to practice[END_REF]. Schemes for initial access based on subspace methods are proposed in [START_REF] Morelli | A robust ranging scheme for OFDMA-based networks[END_REF] and [START_REF] Sanguinetti | An ESPRIT-based approach for initial ranging in OFDMA systems[END_REF]. A solution based on the generalized likelihood ratio test is proposed in [START_REF] Sanguinetti | An initial ranging scheme for the IEEE 802.16 OFDMA uplink[END_REF], whereas [START_REF] Miao | A ranging method for OFDMA uplink system[END_REF] illustrates a RA algorithm that exploits a unique ranging symbol with a repetitive structure in the time-domain.

All the aforementioned schemes are developed by exploiting the degrees of freedom provided by the RACH in both time and frequency domains. In this work, we propose a novel RA protocol that aims at exploiting the spatial degrees of freedom provided by massive MIMO systems as well as the inherent different time instants of reception of UEs' signals (before the data transmission begins). The proposed procedure aims at resolving collisions and also at performing timing, channel, and power estimation. Specifically, it operates as follows. The number of active codes is first determined by using the minimum description length (MDL) principle [START_REF] Wax | Detection of signals by information theoretic criteria[END_REF], while the estimation of signal parameters via rotational invariance technique (ESPRIT) [START_REF] Roy | Esprit -direction-of-arrival estimation by subspace rotation methods[END_REF] is employed for timing recovery. The estimated timing delays are used to compute the leastsquare (LS) estimate of the channels of all identified UEs and also to get their own power levels. The estimated channels are eventually used by the BS in the DL to transmit a response message with "success" notification so as to inform the identified UEs that their random access procedure has been successfully completed.

II. SYSTEM DESCRIPTION

Consider a single-cell TDD massive MIMO system using OFDM as transmission technology. We denote by N FFT the number of subcarriers with frequency spacing ∆f . The BS is equipped with M antennas and K denotes the number of single-antenna UEs that would like to access the network. In doing this, each UE selects an available RACH and sends a request packet to the BS in order to notify the request of network entry. We assume that each RACH consists of Q consecutive OFDM symbols and N adjacent subcarriers such that Q and N are smaller than the (normalized) coherence time and bandwidth of all UEs, respectively. We denote h k ∈ C M the channel frequency response of UE k at the BS within the RACH and assume that

h k ∼ CN (0, β k I M ) (1)
where β k accounts for the large scale channel fading or path loss of UE k. Each UE transmits a randomly chosen code of length Q in parallel over N subcarriers (used for the random access procedure), which is taken from an orthogonal set C = {c 1 , c 2 , . . . , c Q }, with c H i c i = 1 ∀i. The signal transmitted by UE k arrives at the BS with a specific carrier frequency offset (CFO) ǫ k and a normalized (with respect to the sampling period) timing misalignment θ k . Observe that during the initial access, the CFOs are mainly due to Doppler shifts and downlink synchronization errors [START_REF] Morelli | Synchronization techniques for orthogonal frequency division multiple access (OFDMA): A tutorial review[END_REF]. Therefore, they are expected to be significantly smaller than ∆f . In this work, we assume that frequency estimation errors in the downlink are within 2% of ∆f as specified in [START_REF] Sanguinetti | An initial ranging scheme for the IEEE 802.16 OFDMA uplink[END_REF]. In such a case, the impact of CFOs on the random access process can reasonably be neglected. On the other hand, timing errors depend on the distances of UEs from the BS. Their maximum value depends on the propagation environment, and it can reasonably be approximated by the round trip propagation delay of a UE located at the cell boundary. This parameter is known and given by θ max = 2R/(cT s ), where R is the cell radius and T s = 1/(∆f N FFT ) the sampling period. A simple way to counteract the effects of timing errors relies on the use of sufficiently long cyclic prefixes (CPs) comprising N G ≥ θ max +L sampling intervals, with L being the maximum expected delay spread. This leads to a quasi-synchronous network in which timing errors do not produce any interblock interference and only appear as phase shifts at the output of the receive discrete Fourier transform (DFT) unit [START_REF] Morelli | Synchronization techniques for orthogonal frequency division multiple access (OFDMA): A tutorial review[END_REF]. Although such a solution is normally adopted during the random access procedure, the CP of data symbols should be made just greater than the channel length to minimize unnecessary overhead. It follows that accurate timing estimates must be obtained during the random access phase in order to avoid inter-block interference in the data section of the frame.

III. RANDOM ACCESS PROCEDURE

Under the above assumptions, the DFT output y m (n) ∈ C Q at antenna m over subcarrier n during the Q OFDM symbols takes the form for n = 1, . . . , N :

y m (n) = K k=1 √ ρ k h km e -j 2π N FFT (n-1)θ k c u k + w m (n) (2)
where w m (n) ∼ CN (0 Q , σ 2 I Q ) accounts for thermal noise, ρ k > 0 is the power level of the signal transmitted by UE k, and u k ∈ {1, 2, . . . , Q} denotes the index of the code selected by UE k for the random access. In writing the above expression, without loss of generality we have assumed that the first subcarrier of the RACH has index 0. The total power received from UE k in the random access slot is thus given by

p k = ρ k M m=1 |h km | 2 . ( 3
)
Since the number of available codes is limited to Q, collisions may occur between different UEs which are trying to access the network. These collisions need to be detected and resolved before any UE can establish a radio link with the BS for data exchanges. In the next section, we show how vectors {y m (n)} can be used to detect the active codes {c u k ; k = 1, . . . , K} and resolve possible collisions, by exploiting the large number M of antennas and the timing misalignments {θ k ; k = 1, . . . , K} of UEs' signals.

A. Code detection

The first problem is to determine the different code sequences which are being used, and the total number K of users which are transmitting over the RACH. Since at this stage the BS has no knowledge of which codes have been selected, the complete set C must be considered. Therefore, for each possible code c i , the following quantity is computed:

z im (n) = c H i y m (n). (4)
Let K i ⊆ {1, . . . , K} be the subset of UEs that have selected the same code c i and denote by K i its cardinality, i.e., K i = |K i |. Hence, from (2) and (4) we have

z im (n) = i l ∈Ki √ ρ i l h i l m e -j 2π
N FFT

(n-1)θi l + wim (n) [START_REF] Bjornson | Massive MIMO: ten myths and one critical question[END_REF] with wim (n) = c H i w m (n). The samples {z im (n); ∀n} are collected into the vector z im ∈ C N given by

z im = i l ∈Ki h ′ i l m a(θ i l ) + wim (6)
with a(θ i l ) = [1, . . . , e

-j 2π N FFT (N -1)θi l ] T ∈ C N and h ′ i l m = √ ρ i l h i l m .
The above equation indicates that z im has the same structure as the measurement model for a uniform linear array of passive sensors in the presence of multiple uncorrelated sources so that an estimate of K i can be obtained by performing an eigendecomposition (EVD) of the correlation matrix R zi = E{z im z H im }. To see how this comes about, let λ i1 ≥ λ i2 ≥ • • • ≥ λ iN be the eigenvalues of R zi arranged in non-increasing order. Then, from (6) it follows that

λ ik = µ ik + σ 2 k = 1, . . . , r i (7) λ ik = σ 2 k = r i + 1, . . . , N (8)
where r i ≤ K i and µ ik > 0 are, respectively, the rank and the non-zero eigenvalues of the matrix l∈Ki ρ l β l a(θ l )a(θ l ) H . Such a matrix is of rank r i = K i , iff θ i l = θ i k for l = k. Since the timing misalignments are continuous random variables, the probability that θ i l = θ i k for l = k is equal to zero, and hence r i = K i with probability 1. This means that if the true correlation matrix R zi were available all collisions could be in principle resolved provided that K i ≤ N -1. In practice, however, R zi is not available at the receiver and must be replaced with some suitable estimate. One common approach is to use the sample correlation matrix

Rzi = 1 M M m=1 z im z H im (9)
which provides an unbiased and consistent estimate of R zi when M is sufficiently large (as it occurs in massive MIMO systems). Performing the EVD of Rzi and arranging the corresponding eigenvalues λi1 ≥ . . . ≥ λiN in a non-increasing order, we can find an estimate of K i through informationtheoretic criteria. Two prominent solutions in this sense are based on the Akaike and MDL criteria. Here, we adopt the MDL approach, which looks for the minimum of the following objective function [START_REF] Wax | Detection of signals by information theoretic criteria[END_REF]:

Ki = arg min k 1 2 k (2N -k) ln M -M (N -k) ln ĝ(k) (10)
where ĝ(k) is the ratio between the geometric and arithmetic means of { λin ; n = k + 1, . . . , N }.

B. Timing estimation

The ESPRIT algorithm provides a means for estimating the frequencies of multiple complex sinusoidal signals embedded in white Gaussian noise [START_REF] Roy | Esprit -direction-of-arrival estimation by subspace rotation methods[END_REF]. This method belongs to the class of subspace decomposition schemes and exploits the shift-invariant structure of the received signal to get estimates of the unknown parameters in closed form [START_REF] Stoica | Introduction to spectral analysis[END_REF]. Compared to the widely used MUSIC estimator [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], ESPRIT exhibits similar accuracy with a remarkable reduction of complexity as it dispenses with computationally demanding peak search procedures. In the sequel, ESPRIT is applied to the model [START_REF] Sorensen | Massive mimo for crowd scenarios: A solution based on random access[END_REF] to find an estimate of {θ i l ; l = 1, . . . , K i } under the assumption that K i has been correctly estimated, i.e. Ki = K i . Observe that a fundamental assumption behind the ESPRIT estimator is that K i must not exceed N -1. This means that the maximum number of UEs that can simultaneously select a given code c i is limited by N -1.

We begin by arranging the eigenvectors of Rzi associated to the K i largest eigenvalues into the matrix

V i = [v 1 v 2 • • • v Ki ] ∈ C N ×Ki .
The timing offsets {θ i l ; l = 1, . . . , K i } of the K i UEs, which have selected code c i , are thus estimated in a decoupled fashion as

θi l = N FFT 2π arg{ψ il } l = 1, . . . , K i (11)
where {ψ i1 , . . . , ψ iKi } are the eigenvalues of

V i = V (1) H i V (1) i -1 V (1) H i V (2) i ( 12
)
and the matrices V

(1) i and V

(2) i are obtained by collecting the first and the last N -1 rows of V i , respectively. The estimates θi l are then used for identifying the active UEs within K i .

Clearly, the identifiability of the UEs associated to code c i is guaranteed as long as the timing misalignments are different. In other words, two or more UEs that pick the same code c i and are received by the BS with the same time delay are not distinguishable.

C. Channel and power level estimation

Once timing estimation has been performed, the BS must proceed acquiring information about the channels and power levels of the users in K i . From ( 5), the LS estimate of h ′ i l m is found to be

ĥ′ i l m = e T l ÂH i Âi -1 ÂH i z im i l ∈ K i ( 13
)
where e l denotes the lth component of the canonical basis for R Ki , and Âi ∈ C N ×Ki collects the vectors a( θi l ) with i l ∈ K i . Observe that Âi is a Vandermonde matrix. Therefore, the full-rank condition, needed for the computation of the inverse ÂH i Âi , is met if and only if θi k = θij ∀k = j. Assume that the timing offsets have been perfectly estimated, i.e. θi l = θ i l ∀l, so that Âi = A i . Plugging (6) into (13) we obtain

ĥ′ i l m = h ′ i l m + ǫ i l m (14)
where

ǫ i l m ∼ CN (0, σ 2 [A H i A i ] -1 l,l ). The above equation indicates that ĥ′ i l m is an unbiased estimate of h ′ i l m with variance σ 2 [A H i A i ] -1 l,l
. Taking (3) into account, the form of ( 14) suggests the following ad-hoc estimator of the total power:

pi l = M m=1 | ĥ′ i l m | 2 -M σ 2 [ ÂH i Âi ] -1 l,l + i l ∈ K i (15)
where (x) + denotes the maximum between 0 and x. The estimated channels { ĥ′ i l m } are eventually used by the BS in the DL to transmit a response message with a "success" notification to the identified UEs in K i . This message informs them that their random access procedure has been successfully completed, and contains instructions for timing and power adjustment so as to establish a data communication link over a set of specified subcarriers. Notice that pi l may occasionally be null for some i l ∈ K i . This could happen, for example, when the SNR of a given user is too low or when the number of users in K i has been erroneously estimated by the MDL algorithm. In this case, the BS recognizes the anomaly and does not transmit any message to the UEs for which pi l = 0.

IV. NUMERICAL RESULTS

Monte-Carlo simulations have been used to assess the performance of the proposed random access algorithm in terms of code detection capability as well as timing, channel, and power estimation accuracy. The results are obtained for 1000 different channel realizations and UE distributions. We assume that the UEs are uniformly distributed in a circular cell with radius normalized to one and minimum distance 0.1. The path loss function β k is obtained as β k = d -κ k where d k is the distance of UE k from the BS and κ = 3.7 is the path loss exponent. The DFT size is N FFT = 512 and the random access slot is composed of Q = 8 OFDM symbols and N = 16 adjacent subcarriers. Unless otherwise specified, the number of UEs that simultaneously access the network is K = 15, so that K i can be at maximum equal to N -1 = 15. Each UE randomly selects a code from a Walsh-Hadamard codebook of dimension Q = 8. The timing misalignments {θ k } are randomly taken in the interval [0, θ max ] with θ max = 256.

We begin by investigating the performance of the proposed algorithm in terms of probability of correctly identifying the UEs that are trying to get access to the network. Fig. 1 plots the performance of the MDL algorithm as a function of M when the SNR at the cell edge, SNR e , is 0, 5, or 10 dB. As it is seen, the MDL performs very well as the number of UEs is correctly estimated with a probability higher than 90% for all the considered SNR e values when M > 40. As expected, adding more antennas improves the system performance, but at a slow pace. Essentially, the curves in Fig. 1 show the probability that all the collisions within each set K i are correctly resolved by MDL when the true correlation matrix R zi is not available and its estimate Rzi is used instead. In this regard, it is interesting to compare the results of Fig. 1 with the probability of collision in conventional systems in which collision-avoidance entirely rests on the choice of different code sequences (see for example [START_REF] Zhou | OFDMA initial ranging for IEEE 802.16e based on time-domain and frequency-domain approaches[END_REF], [START_REF] Sanguinetti | An initial ranging scheme for the IEEE 802.16 OFDMA uplink[END_REF]). Assuming that QN codes are available (recall that QN is the total number of resources available for our system in the time-frequency domain), the probability that no collision occurs is given by:

p nc = (QN )! (QN -K)! (QN ) -K . (16)
For example, with QN = 128 and K = 15 we have p nc = 0.426. Contrasting this value with the results in Fig. 1, it is seen that, with the proposed procedure, the probability of identifying correctly all the UEs trying to access the network is more than double that of a conventional system. Notice that for the conventional system it has been assumed that two users with different codes can be distinguished independently of their relative power levels (and hence independently of the SNRs). On the other hand, the results in Fig. 1 are more realistic since they take into account the power levels of the different UEs trying to access the network. Fig. 2 illustrates the normalized mean-square-error (NMSE) of the timing and channel estimates versus M for SNR e = 0, 5, and 10 dB. For channel estimation, the normalized mean-square-error is defined as

NMSE h = E M m=1 |h ′ m -ĥ′ m | 2 /( M m=1 |h ′ m | 2 )
where the statistical expectation is taken with respect to the thermal noise, the channel realizations, and the UE distribution. The results of Fig. 2(a) show that the NMSE of the timing estimates decreases fast as M grows large and it is smaller than one sampling interval for M > 40. From Fig. 2(b), it is seen that the NMSE of channel estimates takes small values. Therefore, the BS can reliably make use of the estimated channels in the subsequent DL phase to transmit a response message to the identified UEs.

The results of Fig. 3 investigate the performance of the procedure as a function of K when the ratio M/K = 8 as envisioned in massive MIMO systems. As expected, increasing K deteriorates the performance of MDL since the number of collisions increases and tends to approach the maximum tolerable value given by N -1 = 15. However, a probability of correct estimation as large as 80% is guaranteed for K up to 32.

V. CONCLUSIONS

We proposed a random access algorithm for initial access and handover in the uplink of massive MIMO systems. By exploiting the spatial degrees of freedom offered by massive MIMO systems as well as the inherent different time instants of reception of uplink signals, the proposed solution resolved all the collisions with high probability and, at the same time, performed timing, channel and power estimation with high accuracy. Numerical results demonstrated the robustness of the proposed solution to overloaded situations.

Fig. 1 .

 1 Fig. 1. Probability of estimating correctly all the UEs in each set K i as a function of M when K = 15 and SNRe is 0, 5 or 10 dB.

Fig. 2 .

 2 Fig. 2. NMSE of timing and channel estimates as a function of M when K = 15 and SNRe is 0, 5 or 10 dB.

Fig. 3 .

 3 Fig.3. Probability of estimating correctly the number of UEs in each set K i as a function of K when M/K = 8 and SNRe is 0, 5 or 10 dB.