Growth, Productivity, Biomass and Carbon Stock in Eucalyptus saligna and Grevillea robusta Plantations in North Kivu, Democratic Republic of the Congo

: Initiated by the World Wildlife Fund (WWF) more than a decade ago in North Kivu, single-species plantations of Eucalyptus saligna and Grevillea robusta constitute, with other village plantations, the current legal source of wood-energy for the communities bordering the Virunga National Park (PNVi). This study assesses the growth and productivity of these plantations in two sites with different soil and climatic conditions to predict their production over time. The study also assesses the carbon stock and long-term CO 2 ﬁxation in the biomass of the studied plantations to deduce their contribution to climate change mitigation. Non-destructive inventories were carried out during three consecutive years in 20 E. saligna and 12 G. robusta plantations in Sake and Kirumba. Analysis of the data revealed that both species have similar diametric growth while height growth and productivity were signiﬁcantly higher in the E. saligna plantations. The productivity of E. saligna was also higher in Kirumba than in Sake, while that of G. robusta was higher in Sake than in Kirumba. The differences observed were mainly related to species, silviculture, altitude and concentration of bioavailable elements in the soils. The analysis of productivity evolution over time allowed us to determine optimal rotations at 8 and 12 years, respectively, for E. saligna and G. robusta plantations. The relationships between biomass or carbon stock and tree diameter were not different between the studied species but were signiﬁcantly different at the stand level. If silviculture was standardized and plantations carefully monitored, carbon stock and long-term CO 2 ﬁxation would be higher in G. robusta plantations than in E. saligna plantations. These results indicate that while for productivity reasons E. saligna is the favoured species in wood-energy plantations to quickly meet the demand of the growing and disadvantaged population living in the vicinity of PNVi, carefully monitored G. robusta plantations could be more interesting in terms of carbon credits. To simultaneously optimise wood-energy production and carbon storage in the plantations initiated in North Kivu, E. saligna and G. robusta should be planted in mixture. In addition, species and site characteristics adapted silvicultural management practices must be applied to these plantations, which are very important for the region, its population and its park. Finally, the economic proﬁtability as well as the sustainability of the plantations should be assessed in the longer term in North Kivu.


Introduction
Forest plantations provide a variety of ecosystem services [1,2], contribute to improving people's livelihoods [3,4], provide raw materials for the wood and paper industry [5][6][7], are a source of wood for domestic energy [8,9] and store carbon in their above-and below-ground biomass [10,11]. Currently, large-scale tree planting on non-forested and/or deforested land is seen as one of the ways to restore degraded ecosystems [12][13][14] and to combat global warming [15,16]. High biomass producing forests such as young plantations in the tropics and subtropics can indeed increase the terrestrial carbon sink and thus slow down the accumulation of CO 2 in the atmosphere [10,11,15,16].
Forest plantations are generally composed of one or two species planted at regular spacings and now involve mainly exotic taxa that grow rapidly [5,17]. In the tropics and subtropics, the genera Acacia, Casuarina, Eucalyptus, Gmelina, Grevillea, Leucaena, Pinus or Tectona are generally the most used [9]. Among these species, the genera Eucalyptus and Acacia are amongst those planted in short rotation [18,19]. While the genus Eucalyptus is mainly used in forestry for the production of essential oils, tannins, pulp, timber and woodenergy [20], the genus Acacia is used in agroforestry as well as in forestry for the restoration of degraded land, soil fertility improvement, timber and wood-energy production [19,21,22]. As for the genus Grevillea, it is used in agroforestry as a shade tree for coffee and tea trees, but also in forestry to produce wood-energy over rotations of 10 to 20 years and timber over rotations of 25 to 40 years [23,24].
In the Democratic Republic of Congo (DRC), plantations of fast-growing exotic species have existed since the 1980s but remain poorly documented. These plantations are mainly developed to meet the demand for wood-energy, notably in Kinshasa [25], North Kivu [26], Kwilu [27], Central Kongo [28], Upper Katanga [29] and recently in the Tshopo [30]. Most of these plantations are based on Acacia auriculiformis A. Cunn. ex Benth. except in Kivu where the genus Eucalyptus is mainly used. Other species including Grevillea robusta A. Cunn. and Acacia mearnsii (De Wild.) Pedley are also planted in Kivu [31]. These plantations are generally monospecific, sometimes in association with food crops. Apart from the A. auriculiformis plantations in Mampu on the Bateke plateau, the other plantations in the DRC have not yet been the subject of scientific studies, particularly with regard to production and carbon sequestration potential, economic profitability or the sustainability of their management. Among these plantations are those initiated since 2007 as part of the EcoMakala project around the Virunga National Park (PNVi) by the World Wildlife Fund (WWF). These monospecific plantations based on Eucalyptus spp., G. robusta or A. mearnsii currently cover nearly 11,200 ha [32] and are expected to play a major role in the protection of the PNVi by providing wood-energy to the population living near the park, including the town of Goma, and by providing an alternative financial income to the farmer-planters. In 2011, the EcoMakala project was approved and registered as a geographically integrated REDD+ pilot project in eastern DRC [33]. Between 2020 and 2021, EcoMakala+ received USD 1.3 million, making North Kivu the first province to benefit from the sale of carbon credits in DRC [34].
To establish recommendations for management adapted to the objectives and expected production of these plantations, it is important to characterise and model the determinants of their productivity as suggested by Rondeux [35]. To this end, growth and productivity models must be fitted at the tree or stand level. Vanclay [36] defines a growth model as a synthesis of inventory data to better characterise the growth and evolution of forests. Growth models are generally equations or systems of equations that predict the growth and yield of a forest stand for given ecological conditions. These empirical models are usually Forests 2022, 13, 1508 3 of 24 calibrated by combining inventory data with environmental data (sometimes also remote sensing data) and are tools for forest researchers and managers to predict production and explore silvicultural options [35,36]. The validity of the models obtained is often restricted to the conditions under which the data were collected. Given the spatial and temporal variability of production factors (species, age, seed source, ontogenic stage, silviculture, climate, soil type, topography), it is indeed risky to use such models to extrapolate to conditions other than those inventoried.
Given the importance of the E. saligna and G. robusta plantations initiated around the PNVi, particularly in terms of wood-energy production and expected carbon remuneration, the main purpose of this research is to quantify their productivity and carbon stock so as to propose adapted management measures. Specifically, in the case of such plantations installed in two sites with contrasting soil and climate conditions (Sake and Kirumba), the study aims to provide answers to the following research questions: (i) What are their growth and productivity? (ii) Can we define/design/establish a silviculture adapted to the biotic and abiotic factors that significantly influence their growth and productivity? (iii) What should be the rotation length as well as the age and intensity of thinning to maximise their production? (iv) To which extent are those plantations contributing in the mitigation of climate change and what is their potential income on the carbon markets?
Based on the literature review summarised in Table 1, the basic hypotheses tested in this study were that productivity and carbon stock are higher in plantations (i) of E. saligna rather than G. robusta, and (ii) on fertile volcanic soils of Sake (see [37,38]) rather than non-volcanic ones of Kirumba.

Study Sites
The study was conducted in the province of North Kivu (DRC), specifically in Sake and Kirumba (Figure 1). The relief is essentially mountainous in both sites with an altitude varying between 1470 and 1760 m in Sake and between 1650 and 2200 m in Kirumba (Table S1). The soils of Sake are derived from ancient (notably at Kimoka, Kirotshe and Luhonga) and recent (notably at Mubambiro) basaltic lava [37,39]. These soils are slightly acidic or neutral (pH-water = 6.8 ± 0.4; Table S2). Soils developed on recent lavas differ from those derived from ancient lavas in their low cation exchange capacity (CEC = 3.5 ± 0.4 cmolc/kg at Mubambiro compared to 14.0 ± 9.9 cmolc/kg at Kimoka, 42.7 ± 15.3 cmolc/kg at Kirotshe and 15.5 ± 0.9 cmolc/kg at Luhonga) and bioavailable calcium concentration ([Ca]: 2200 ± 420 µg/g in Mubambiro compared to 4420 ± 1960 µg/g in Kimoka, 7720 ± 2080 µg/g in Kirotshe and 3180 ± 80 µg/g in Luhonga). In Kirumba, the soils are essentially derived from granitic rocks [39]. These soils are acidic with high and locally sometimes toxic aluminium content (pH-water = 4.6 ± 0.4 and [Al] = 620 ± 370 µg/g), low cation exchange capacity and low bioavailable calcium concentration (respectively, CEC = 3.5 ± 2.1 cmolc/kg and [Ca] = 280 ± 290 µg/g). According to the FAO WRB classification, Sake soils are Haplic Acrisols while Kirumba soils are Aluandic Andosols [40]. The Sake plantations benefit from a tropical monsoon climate (Am) with an average temperature of 19.9 • C and 2716 mm/year of rainfall. In Kirumba, the climate is equatorial (Af) with an average temperature of 18.4 • C and rainfall of 3750 mm/year [41]. In both sites, the vegetation consists of a mosaic of forest and savannah [42]. Eucalyptus spp. and G. robusta plantations are established in small clumps on private concessions. average temperature of 18.4 °C and rainfall of 3750 mm/year [41]. In both sites, the vegetation consists of a mosaic of forest and savannah [42]. Eucalyptus spp. and G. robusta plantations are established in small clumps on private concessions.

Study Material
The study was conducted in 32 single-species plantations, of which 20 were of E. saligna (Photo 1a) and 12 of G. robusta (Photo 1b). Literature data related to the botanical, ecological, silvicultural and cultural characteristics of these two species originating from Australia are summarised in Table 1.

Study Material
The study was conducted in 32 single-species plantations, of which 20 were of E. saligna ( Figure 2a) and 12 of G. robusta (Figure 2b). Literature data related to the botanical, ecological, silvicultural and cultural characteristics of these two species originating from Australia are summarised in Table 1.   The plantations studied were selected to provide a sample of the different ages, planting spacings and topographical situations per species and site. These plantations were established between 2010 and 2016 by transplanting 15-20 cm tall seedlings at spacings of 3 m × 3 m, 2.5 m × 2.5 m or 2 m × 2 m for E. saligna and 4 m × 4 m or 3 m × 3 m for G. robusta. During the first three years, the trees are regularly hoe-weeded (twice or thrice a year) and food crops are usually grown in association with them. From the fourth year onwards, maintenance is done with a machete (once or twice a year). Agricultural practices as well as extensive grazing (cattle, goats and sheep) have been observed in some plantations. In most of the plantations studied, stems were not evenly distributed. This irregularity in the spatial distribution of stems was due to a lack of replanting following, depending on the plantation, malicious acts (e.g., theft of seedlings), natural mortality, the passage of a bush fire and/or the early cutting of some stems by thieves. The success rate (see Equation (5)) varied between 30 and 100% in E. saligna plantations and between 18 and 94% in G. robusta plantations (Table S1). The lowest success rates were observed in the Sake plantations (18%-70%) and the highest in the Kirumba plantations (52%-100%).
Thinning and/or pruning operations were rare in these plantations. They are mainly exploited by clear-cutting to produce wood-energy (mainly charcoal) and construction wood (planks, rafters and sleepers). This cutting generally takes place between 6 and 10 years in E. saligna plantations and between 10 and 15 years in G. robusta plantations. However, some E. saligna plantations dedicated to timber production remain beyond 10 years and are often thinned between 8 and 10 years. Thinned wood is generally used for energy purposes (fuelwood and charcoal). Table 1. Botanical, ecological, silvicultural and cultural characteristics of the tree species used in the studied plantations.

Trunk and bark
The bole is generally straight and free of branches for half or two-thirds of the height. Its bark is rough and persistent, brownish or greyish and often peels off in long strips The bole is generally without buttresses. The bark is grey to light brown, becoming furrowed with age

Inventory Design and Data Collection
Non-destructive inventories were carried out periodically in November-December 2018, 2019 and 2020 in 20 E. saligna plantations and 12 G. robusta plantations. Particularly in the G. robusta plantations, a fourth data collection campaign took place in May 2022 because the three-yearly data were judged not sufficient to accurately monitor tree growth. The area (in hectares, ha) and centre of each plantation were determined with a GPS (Garmin GPSMAP 64stc, Olathe, KS, USA). This area varied between 0.49 and 3.0 ha. In each sampled plantation, the inventory was carried out on a constant area of 500 m 2 distributed over three circular plots, each with a radius of 7.30 m measured horizontally to correct for slope. The layout of the plots was adapted to the shape of each plantation. The first plot was always placed in the centre of the plantation (Figure 3). The other two plots were then positioned at approximately half the distance between the centre and the boundary of the plantation along the long side. Each plot was centred on a tree. All trees in the plot were measured and marked and their polar coordinates were determined with reference to the centre tree. As the area of the plantations varied between 0.49 and 3.0 ha, the sampling rate was therefore between 1.67 and 10.2%.  The diameter at 1.30 m from the ground (DBH in cm) and the total height (Ht in m) of each tree in the plot were measured with a dendrometric tape and a laser rangefinder (Forestry Pro), respectively. The DBH measurement level was marked with sticky paper which was then removed after marking the tree with oil paint [49,50]. A total of 1044 trees were measured, of which 763 were in E. saligna plantations and 281 in G. robusta plantations.
The topographical situation of each plot was characterised using the classification of Méthot et al. [49]. The slope of the land was measured with a Suunto clinometer. Due to the absence of weather stations at both sites, climate variables (temperature and precipitation) were acquired from global online data [41] using the geographical coordinates of The diameter at 1.30 m from the ground (DBH in cm) and the total height (Ht in m) of each tree in the plot were measured with a dendrometric tape and a laser rangefinder (Forestry Pro; Nikon, Tokyo, Japan), respectively. The DBH measurement level was marked with sticky paper which was then removed after marking the tree with oil paint [49,50]. A total of 1044 trees were measured, of which 763 were in E. saligna plantations and 281 in G. robusta plantations. The topographical situation of each plot was characterised using the classification of Méthot et al. [49]. The slope of the land was measured with a Suunto clinometer (Vantaa, Finland). Due to the absence of weather stations at both sites, climate variables (temperature and precipitation) were acquired from global online data [41] using the geographical coordinates of each plantation.
Soil samples were taken with a soil auger to a depth of 30 cm. For each plantation, two composite samples of 400 g each were taken, one from samples taken in the plots and the other from samples taken in neighbouring fields and/or fallows at a minimum distance of 30 m from the plantation boundaries. These samples were each packed in a plastic bag, marked and sent to the Laboratory of Plant Ecology and Biogeochemistry of the Université Libre de Bruxelles (ULB) for physico-chemical analyses (Table S2). The granulometry (proportion of sand, silt and clay) was determined by sieving and sedimentation after destruction of the organic matter. The pH-water was measured with glass electrodes in a soil suspension according to NF-ISO 10390. The cation exchange capacity (CEC, in cmolc/kg) was evaluated on soil extracts with cobaltihexamine trichloride according to ISO 23470:2007 [51]. The bioavailable elements (Ca, Mg, K, Al, Cu, Fe, Mn, P and Zn, in µg/g) were extracted with 0.5 M EDTA 0.02 M ammonium acetate solution at pH 4.65 [52]. Comparison of the values of these parameters showed that the soils collected from the plantations did not differ significantly (Student's t-test, p-value > 0.05) from those of the surrounding fields. At this stage, the soil variables do not seem to be affected by the plantations.
Wood cores for the infradensity measurements required for biomass estimation were collected in 2019 from 10 plantations ranging in age from 6 to 10 years, with 5 plantations per species. In each selected plantation, four sample trees were selected along the four cardinal points (E-W-N-S) around the central plot. Wood cores were extracted at 1.30 m from the ground on either side of the tree trunk (180 • apart) using a Pressler auger that was driven into the tree perpendicular to the trunk axis [53]. Each of the 80 wood cores collected was stored in a paper envelope and marked. All wood cores were then sent to the Wood Biology Laboratory in Yangambi (DRC) for infradensity measurements. Each core of diameter D (in cm) and length L (in cm) was placed in an oven for drying at 105 • C for 48 h [54]. Once dried, the core was weighed on a 0.001 g precision balance. The infradensity (ρ in g/cm 3 , Equation (1)) was calculated by dividing the dry mass of the core (DM in g) by its wet volume (V c in cm 3 , Equation (2)).

Data Analysis
The data were manipulated with Excel (version 2013, Microsoft Corporation, Redmond, WC, USA) and analysed with R (version 4.2.0, R Core Team, Vienna, Austria) [55].

Estimation of Dendrometric Characteristics
The characteristics recorded for each plantation are: age (in years); stand density (N in stems/ha, Equation (3)); basal area (g i and G in m 2 /ha, Equations (6) and (7)); mean diameter (DBH in cm); mean total height (Ht in m); dominant height (Hdo in m); total volume (V in m 3 /ha, Equation (9)); quadratic diameter (d g in cm, Equation (10)); mean tree volume (v m in m 3 , Equation (11)); and mean basal area (g m in m 2 , Equation (12)). From the observed density (N) and the theoretical density (N , Equation (4)) derived from the spacing used at planting (L and e in m, cf. Table S1), the success rate was calculated for each plantation (SR in %, Equation (5)). The average diameter and the average total height were respectively calculated as the arithmetic mean of the diameters and total heights of all trees measured in the plots. The dominant height is the average height of the 100 largest trees in one hectare. It was calculated as the arithmetic mean of the total heights of the five largest trees in the inventory area (500 m 2 ). The total volume was calculated as the sum of the individual volumes of the trees in the plot. The individual volume (Vi in m 3 ) was estimated with the allometric equation proposed by Deleuze et al. [56] where 0.496 is the shape coefficient (Equation (8)). The quadratic diameter corresponds to the diameter of the tree with mean basal area g m and mean volume v m .

Estimation of Growth and Productivity
Growth and productivity were assessed by calculating the average annual increments in diameter (ADBH in cm/year), total height (AHt in m/year) and volume (AMV in m 3 /ha/year). Each of these increments corresponds to the ratio between the size observed at a given time and the age of the plantation at that time [35].

Comparison of Dendrometric Characteristics
The averages of the dendrometric characteristics were calculated for each species/site and compared between them ("unpaired Student's parametric t-test") with a threshold α of 0.05 (Table 2; Figures 4-6). Plantations were previously classified into three categories based on the average age calculated over the three inventories and the species (Figure S1) using the hclust function of the vegan package [57]. The first category includes plantations aged ≥ 8 years; they are referred to as "older plantations". The second category consists of plantations aged between 5.5 and 7.9 years (Intermediate-aged plantations). The third category gathers plantations ≤ 5.4 years old (younger plantations).

Predictions of Increment, Basal Area and Carbon Stock
To define suitable silvicultural itineraries and to identify the biotic and abiotic factors influencing the growth and productivity of the studied plantations, models based on multiple linear regressions (MLR) were developed in two steps. In the first stage, increments were modelled as a function of stand age and density to allow managers to easily predict stand growth and productivity. These two variables were also used in the predictions of basal area and carbon stock, so as to set the age of thinning and estimate the carbon stock over the entire rotation, respectively. In order to identify the variables that significantly influence the growth and productivity of the two species in the two sites, the increments were, in the second step, modelled by combining all measured biotic and abiotic variables (site, stand age and density, altitude, slope, mean temperature and total annual precipitation, soil grain size and bioavailable elements concentrations in soils). These variables were all previously standardised according to Baillargeon [58]. The selection of variables was carried out step by step by minimising the Akaike information criterion (AIC), the root mean square error (RMSE) and the variance inflation factor (VIF) [59,60]. The model selected was one in which the explanatory variables were not redundant or collinear, with normally distributed (shapiro.test, p-value > 0.05) and homoscedastic (ncvTest, p-value > 0.05) residuals [61][62][63]. The collinearity between the variables was tested using the vif function of the car package [64]. Variables with a VIF greater than 5 were progressively eliminated from the model starting with the one with the highest VIF [63][64][65]. The quality of the predictions was tested by comparing the variance of the observed data with the variances of the predicted data using the fligner.test function in the stats package [55].

Expected Productivity over the Rotation, Age and Length of the Thinning Interval
According to Lanier [66], the rotation is the number of years between the establishment of a forest stand and its final felling. Thinning is the elimination of part of a forest stand in order to allow the best trees to develop. The length of the thinning interval (hereafter "thinning interval") is the time (in years) between two successive thinning passes on the same plot. In this study, the rotation was determined in such a way as to maximise the average increment in volume (Figures 4c, 5c and 6c). The expected productivities on this rotation were estimated from Models 3 and 8 (Table 3) for E. saligna and G. robusta, respectively. Given the temperament of these two species (both heliophilous) and the age of the stands (between two and 12 years), thinning was considered from the moment when the total basal area reached 30 m 2 /ha in accordance with Prégent [65]. This choice is justified by the fact that basal area is a good indicator of the degree of competition between stems and takes into account both the number and size of stems and the quality of the site [65,67]. The age of first thinning for each site and species was determined from the threshold basal area (30 m 2 /ha) and stand density using Models 4 and 9 (Table 3) for E. saligna and G. robusta, respectively. Thinning intensity was estimated from Equation (13) [65] assuming a thinning that removes 10 m 2 /ha of basal area. Under these conditions, the residual basal area (after thinning) is set at 20 m 2 /ha. The thinning interval (Table 4) was calculated from the age difference between two consecutive thinnings. The number of stems to be removed during the thinning was determined by multiplying the thinning intensity by the stand density before the thinning. The variability of the age of first thinning as a function of site, species and stand density was tested by MLR.
Thinning intensity = (threshold basal area − residual basal area) threshold basal area × 100 To enable decision-makers to define the area to be planted, projections were made (Equation (14)) based on the expected productivity and volume of wood-energy needed to satisfy the demand of an arbitrary sample of one million inhabitants (case of the city of Goma) with reference to FAO estimates for which the average annual individual consumption in developing countries is one cubic meter of wood per inhabitant [9].
Area to be planted (in ha) = volume of wood to be produced m 3 /an expected productivity (m 3 /ha/an)

Biomass, Carbon Stock and Financial Potential of Carbon Stored by the Plantations Studied over the Rotation Period
Above-ground biomass was estimated at the tree and stand level. At the tree scale, biomass was estimated from Equation (15) proposed by Chave et al. [68]. (15) where AGB est = above-ground biomass (in kg), D = DBH (in cm), H = Ht (in m) and ρ = wood infradensity (in g/cm 3 ).
At the stand level, the above-ground biomass (AGB in Mg/ha) is the sum of the individual biomasses. The below-ground biomass (BGB in Mg/ha) was calculated by multiplying the above-ground biomass by 0.2 as suggested by Ponce-Hernandez [69]. The total biomass was obtained by summing the AGB and BGB. The carbon stock (StC, in Mg/ha) was estimated by multiplying the total biomass by 0.5 as recommended by the IPCC [70]. The total carbon stock over the rotation was estimated from models 5 and 10 (Table 3) for E. saligna and G. robusta, respectively. The CO 2 equivalent was calculated by multiplying the carbon stock by 3.67 [70]. The long-term fixed CO 2 was estimated by dividing the CO 2 equivalent by the rotation time [71]. The financial potential due to the CO 2 fixation by the studied plantations was estimated by multiplying the long-term fixed CO 2 by the international market price of forest carbon, i.e., USD 18 per Megagram of CO 2 on 30 June 2022 according to the Gold Standard [72]. The biomass (or carbon stocks) of E. saligna and G. robusta plantations were compared with each other using Student's t-test. The variability of biomass or carbon stock with site, species, age and stand density was tested by MLR.

Dendrometric Characteristics
All dendrometric characteristics assessed differed significantly by species in the intermediate-aged plantations ( Table 2). In the younger and older plantations, only the height (Ht, Hdo and AHt) and volume (V and AMV) characteristics were significantly species dependent. For all age groups combined, the average density assessed in E. saligna plantations (724 ± 303 stems/ha) was significantly higher (p < 0.001) than that found in G. robusta plantations (461 ± 216 stems/ha). The ADBH of E. saligna (2.6 ± 0.6 cm/year) was similar (p = 0.121) to that of G. robusta (2.4 ± 0.5 cm/year). In contrast, AHt and AMV of E. saligna (3.0 ± 0.6 m/year and 30 ± 13 m 3 /ha/year, respectively) were significantly higher (p < 0.001) than those of G. robusta (2.1 ± 0.3 m/year and 15 ± 7 m 3 /ha/year, respectively). ADBH and AHt of both species decreased as the stands matured (Figure 4a,b), with maximum values observed between 5 and 7 years for E. saligna and between 3 and 5 years for G. robusta. The maximum value of AMV was observed between 7 and 8 years for E. saligna and between 10 and 12 years for G. robusta (Figures 4c, 5c and 6c). Based on these observations, the rotation that maximises production is 8 years for E. saligna plantations and 12 years for G. robusta plantations. In addition, the average increments of E. saligna were higher in Kirumba than in Sake ( Figure 5), while those of G. robusta were higher in Sake than in Kirumba ( Figure 6).

Variation Factors for Growths
The growth and productivity of the plantations studied were significantly dependent on age, stand density and site characteristics ( Table 3). The fitted models showed that the combined effect of age and stand density explained 28, 30 and 29% of the variability in diameter, total height and volume increments for E. saligna (Models 1-3) and 69, 50 and 59% for G. robusta (Models 6-8), respectively. The addition of environmental predictors to the models refined the predictions by reducing the AIC and RMSE and explaining 48, 52 and 76% of the variability in diameter, total height and volume increments for E. saligna (Models 11-13) and 74, 75 and 87% for G. robusta (Models 14-16), respectively. * We preferred to present models fitted with raw data (with age expressed in 'years' and stand density in 'stems/ha') to allow forest managers to easily predict stand growth and productivity without resorting to data transformation.

Thinning and Expected Productivity over the Rotation
Two scenarios were tested for the 8-year rotation for E. saligna plantations and the 12-year rotation for G. robusta plantations. The first scenario considers the production of wood-energy in plantations with different silviculture depending on the species and generally little monitoring as observed in the field (see Section 2.2). The second scenario aims to produce timber from plantations of identical silviculture species in which thinning removes 10 m 2 /ha when the basal area of the plantation reaches 30 m 2 /ha, i.e., a thinning intensity of 33%. In both scenarios, the rotation was set at 8 years in E. saligna plantations and 12 years in G. robusta plantations. The application of Models 4 and 9 revealed that in scenario 1, the age of the first thinning is well beyond the respective rotations (Table 4). In scenario 2, one or two thinnings would be required over the 8-and 12-year rotations, respectively, in E. saligna and G. robusta plantations installed at 2 × 2 m and/or 2.5 × 2.5 m spacings. The thinning interval would be 3 to 4 years in E. saligna plantations and 4 to 5 years in G. robusta plantations. The age of first thinning varied significantly (MLR: adjusted R 2 = 0.96; p < 0.001) by site (p = 0.002), stand density (p < 0.001) and site-species interaction (p < 0.001). Table 4. Age, interval and intensity of thinning in E. saligna and G. robusta plantations in North Kivu. The age of first thinning was determined for each site from the threshold basal area (30 m 2 /ha) and stand density using Models 4 (for E. saligna) and 9 (for G. robusta). The asterisks indicate that the stand density after thinning is that of a stand: unthinned '*', thinned once '**', thinned twice '***'.  (Table 5) over the entire rotation confirmed that the productivity of E. saligna plantations was higher than that of G. robusta plantations, regardless of the site and scenario considered. In Kirumba, the ratio of E. saligna to G. robusta productivity was 3.0 in the first scenario and 2.0 in the second. In Sake, the ratio was 1.5 in the first scenario and 1.2 in the second. Based on these forecasts, meeting the annual wood-energy demand of an estimated one million inhabitants would require between 20,000 and 32,000 ha of monospecific E. saligna plantations or between 27,000 and 72,000 ha of monospecific G. robusta plantations.

Wood Infradensity, Biomass, Carbon Stock and Economic Potential of Carbon Sequestered on Rotation
In plantations ranging in age from 6 to 10 years, the average infradensity was 0.56 ± 0.06 g/cm 3 for E. saligna wood and 0.59 ± 0.01 g/cm 3 for G. robusta wood. The biomass (and thus the carbon stock) of the two species did not differ significantly (Student's test: p = 0.981) for trees of the same size (DBH, Figure 7a) although their heights differed. At the stand level (Figure 7b,c), biomass or carbon stock varied significantly (MLR: Adjusted R 2 = 0.73; p < 0.001) by species (p = 0.029), site (p < 0.001), age (p = 0.028), stand density (p = 0.003) and site-species interactions (p < 0.001) on the one hand, and stand age and density on the other (p < 0.001). Table 5. Expected productivity over the rotation and area to be dedicated to plantations. Expected productivities were calculated from Models 3 and 8 ( Table 3). The area to be dedicated to plantations was estimated from Equation (14). The rotation stand density corresponds to the stand density after thinning (see Table 4).   (Table  6). Long-term CO2 fixation over the respective rotations is estimated to be between 50 and 65 Mg CO2-equivalent per hectare in the E. saligna plantations and between 30 and 100 Mg CO2-equivalent per hectare in the G. robusta plantations. The financial potential of carbon credits over the entire revolution is estimated to be between USD 900 and USD 1200 per hectare in E. saligna plantations and between USD 500 and USD 1900 per hectare in G. robusta plantations. Table 6. Carbon stock, long-term CO2 fixation and financial potential of carbon credits in E. saligna and G. robusta plantations in North Kivu. The carbon stock was calculated from models 5 (for E. saligna) and 10 (for G. robusta). Long-term fixed CO2 and the financial potential of carbon credits were deducted from the carbon stock following the methodology described in Section 2.4.6. The rotation stand density corresponds to the stand density after thinning (see Table 4).  (Table 6). Long-term CO 2 fixation over the respective rotations is estimated to be between 50 and 65 Mg CO 2 -equivalent per hectare in the E. saligna plantations and between 30 and 100 Mg CO 2 -equivalent per hectare in the G. robusta plantations. The financial potential of carbon credits over the entire revolution is estimated to be between USD 900 and USD 1200 per hectare in E. saligna plantations and between USD 500 and USD 1900 per hectare in G. robusta plantations. Table 6. Carbon stock, long-term CO 2 fixation and financial potential of carbon credits in E. saligna and G. robusta plantations in North Kivu. The carbon stock was calculated from Models 5 (for E. saligna) and 10 (for G. robusta). Long-term fixed CO 2 and the financial potential of carbon credits were deducted from the carbon stock following the methodology described in Section 2.4.6. The rotation stand density corresponds to the stand density after thinning (see Table 4).

Different Dendrometric Characteristics According to Species, Silviculture and Site Characteristics
The results of this study showed that in Sake as in Kirumba, E. saligna is a much more productive species than G. robusta. In addition, the productivity of E. saligna was higher in Kirumba than in Sake ( Figure 5), while that of G. robusta was higher in Sake than in Kirumba ( Figure 6). Despite the differences, the average productivities observed (30 ± 13 m 3 /ha/year at the age of 6.2 ± 2.1 years for E. saligna (N = 724 ± 303 stems/ha) and 15 ± 7 m 3 /ha/year at the age of 7.7 ± 2.3 years for G. robusta (N = 461 ± 216 stems/ha)) are comparable to the averages reported in other sites. Indeed, for E. saligna, Dyson [73] reports productivity varying between 17 and 39 m 3 /ha/year over 10 years in plantations in Kenya (N = 509 stems/ha). For the same species, Walters [74] reports productivity of between 33 and 46 m 3 /ha/year over 15 years in plantations in Hawaii (N = 247 to 400 stems/ha). Concerning G. robusta, Pandey [75] reported a productivity of 10 to 12 m 3 /ha/year over 10 to 15 years (N = 800 to 1200 stems/ha) in India, while Muchiri et al. [76] found productivities varying between 8 and 24 m 3 /ha/year over 30 years (N = 200 stems/ha) in some Kenyan agroforestry systems.
As the two species have a similar diameter increment (2.6 ± 0.6 cm/year for E. saligna and 2.4 ± 0.5 cm/year for G. robusta), the observed difference in interspecific productivity is mainly the result of a greater increment in height for E. saligna  [11][12][13][14][15][16]. Climatic conditions (temperature and rainfall) can also influence growth and productivity [77,78], but these two variables were not determinant in our study where rainfall (3750 mm/year in Kirumba and 2716 mm/year in Sake) was well above the water requirements of both species (800-1800 mm/year for E. saligna and 600-1700 mm/year for G. robusta according to Orwa et al. [45,48]). These results support findings that forest stand growth or productivity varies with species, silviculture and ecological site conditions (e.g., [36,[79][80][81][82][83]).
The difference in productivity of E. saligna between Sake and Kirumba was, according to Model 13 (Table 3), mainly related to the positive effects of altitude and stand density, which were higher in Kirumba (respectively, 1580 ± 70 m in Sake versus 2010 ± 190 m in Kirumba, and 685 ± 343 stems/ha in Sake versus 786 ± 217 stems/ha in Kirumba). In addition, E. saligna was more productive on acidic and toxic soils in Kirumba than on slightly acidic or neutral soils in Sake because the species is well adapted to both acidic and basic well-drained soils in tropical and subtropical regions [43,45]. The decrease in productivity of E. saligna in Sake can finally be explained by the high calcium concentrations in the volcanic soils of the region ([Ca] = 4400 ± 1970 µg/g in Sake compared to 900 ± 2010 µg/g in Kirumba). Indeed, high calcium concentrations lead to the retrogradation of phosphorus and the insolubilisation of micronutrients, including manganese, and prevent their assimilation by plants [84,85]. As for G. robusta, its productivity in Sake is, according to Rojas-Sandoval [24], favoured by the presence of volcanic soils, rich in exchangeable cations with a pH favourable to their assimilation by the trees. The low productivity of this species in Kirumba was due to the increase in altitude and the aluminium toxicity of the soils (Model 16).
These results support the conclusion that tree plantations initiated in North Kivu require management measures adapted to the edaphic conditions of the sites [86]. While the soils of Sake have a pH and CEC favourable to tree growth and productivity, the acidic and aluminium-rich soils of Kirumba theoretically require mineral or organic amendments. Given the negative effect of calcium on micronutrient availability in E. saligna, liming is not recommended. Furthermore, given the socio-economic, logistical and technical contexts of the region, farmer-planters do not have sufficient access to manures and/or composts for organic amendments. Itinerant livestock rearing could be beneficial in the region, but it is only practised by a minority in the Sake plantations. To this end, the introduction of soil-improving species in mixtures should be favoured, especially nitrogen-fixing species such as those of the Acacia, Albizia, Calliandra, Casuarina and Leucaena genera. The gradual replacement of monocultures would eventually make up for some of the weaknesses of this type of planting [87,88]. These species can also be used in mixtures in Sake plantations to combat apatitic retrogradation linked to the high presence of calcium in E. saligna. However, experimental plantations based on these species are first needed in the region to: (i) ensure their adaptability or not to local ecological conditions; (ii) check that the benefits expected according to Kelty [87] are actually observed; and (iii) serve as controls with farmer-planters and/or donors. In addition, an awareness and information campaign on the benefits of these species and/or mixed plantations should be conducted among farmers, managers and donors to ensure their acceptability as suggested by Messier et al. [89]. Possible concerns about the impact of these species on soil water resources are discounted [18,[90][91][92][93][94] as the region benefits from sufficient rainfall (>2500 mm/year).

Empirical Models to Predict Growth and Productivity
Our study developed predictive models for growth and productivity of E. saligna and G. robusta plantations based on age, stand density [95,96], but also on the combination of selected biotic and abiotic factors [36]. The analysis of the residuals confirmed the good fit of the models to the data (Table 3). Furthermore, the increases predicted from stand age and density did not differ significantly (Fligner and Killeen test, p > 0.05) from those observed or from those predicted by combining biotic and abiotic factors (Figures S2 and  S3). Nevertheless, the models presented were developed under the conditions of the study: first rotation, monospecific even-aged stands, agricultural precedent, soil fertility inherited from the original ecosystem (confirmed by comparison with surrounding soils), etc. If these models have been shown to be effective under these conditions, it would be interesting to verify their extent/scope of application. Unfortunately, in addition to the cost and time involved in such a large survey, local security conditions, the geographical area to be covered, the extreme fragmentation of existing plantations and the poverty in the region do not favour the periodic collection of data, by disrupting the work schedule and obliging planters to cut their plantations prematurely. However, as data on site characteristics are still not readily available, managers can now predict the productivities of even-aged stands of E. saligna and G. robusta established in Sake and Kirumba from stand age and density using Models 3 and 8, respectively.

Results That Suggest Ages of the First Thinning and Rotations Adapted to the Species, Site and Silviculture
Analysis of the evolution of stand productivity over time allowed us to set the rotation at 8 years for E. saligna plantations and at 12 years for G. robusta plantations. These rotations are comparable to those used in other sites for both species (e.g., [45,48,97,98]) and are within the range of cutting ages intuitively used by planters in North Kivu (see Section 2.2). Furthermore, the rotation set for E. saligna is similar to that applied in the A. auriculiformis agroforestry systems of Mampu (8-10 years) on the Bateke Plateau in the DRC [99,100].
Concerning thinning, our results showed that the age and periodicity of intervention varies with species, site and spacing (Table 4) as indicated by Prégent [66]. In carefully monitored E. saligna plantations, thinning to remove 33% of the stems should be carried out from the age of 4, 5 or 6 years, depending on the site and spacing, and could continue at intervals of 3 to 4 years if the stand remains standing. These thinning practices are comparable to those applied in plantations of the same species in South Africa where timber stands established at densities of 1330 stems/ha are generally thinned by 25% of stems from the age of 6 years, and then these thinnings are continued at intervals of 3 to 5 years until the final harvest of the wood [43]. In G. robusta plantations, thinning of 33% of stems should be applied with intervals of 4 or 5 years from the age of 6 years if stands are established at 2 × 2 m spacings or somewhat later from 10 years or more for stands established at spacings ≥ 2.5 m. These ages and thinning intervals are comparable to those applied in Brazil [23] and Rwanda [101] where plantations of the same species are thinned by about 30% of stems between 10 and 15 years, with an interval of 5 years.
Referring to the productivity forecasts (Table 5), it can be seen that in addition to the wood that would result from thinning operations, plantations established at close spacings (2 to 2.5 m) and scrupulously monitored could produce a volume of raw wood greater than or equal to that of unthinned plantations established at wide spacings. This result indicates that to optimise wood-energy production on small plots in North Kivu, E. saligna and G. robusta should be planted at spacings between 2 and 2.5 m. Spacing of less than 2 m would not be promoted in the two study sites for both species as they would require early thinning with very low individual woody biomass production. However, for timber production purposes, spacings of 3 to 4 m would ideally be of interest.

Similar Interspecific Biomass and Carbon Stocks at the Tree Level but Significantly Different at the Stand Level
For the same DBH, the above-ground biomass (and therefore the carbon stock) of an individual E. saligna was similar to that of an individual G. robusta (Figure 7a). This similarity in individual biomass between the two species is explained by the compensation of the difference in height by the difference in wood infradensity. Indeed, in E. saligna and G. robusta plantations of similar age and DBH, E. saligna individuals were generally taller with a slightly low wood infradensity (ρ = 0.56 g/cm 3 ) while G. robusta individuals were generally less tall with a slightly high wood infradensity (ρ = 0.59 g/cm 3 ). The wood infradensities found for both species are comparable to those reported between 5 and 10 years in other sites for the same species (e.g., [44][45][46]48]). At the stand level, biomass and carbon stock in age-matched plantations were significantly higher for E. saligna than for G. robusta (Figure 7b,c). This result is justified by a generally higher productivity and stand density in E. saligna plantations than in G. robusta plantations. Considering the 8-year rotation for E. saligna plantations and the 12-year rotation for G. robusta plantations, the total biomass was, however, similar (ca. 230 Mg/ha) for both species. This average biomass production is higher than that reported over 10 years in A. auriculiformis plantations (145 Mg/ha) on the Bateke plateau in the DRC [99].

What Choice Should Be Made between E. saligna and G. robusta?
In the two scenarios analysed, the average productivity of E. saligna plantations over the 8-year rotation was significantly higher than that of G. robusta plantations over the 12-year rotation, with differences of between 12 and 28 m 3 /ha/year in the first scenario and between 6 and 22 m 3 /ha/year in the second ( Table 5). As for the carbon stock (Table 6), the first scenario showed that it was similar in plantations of both species (ca. 117 Mg/ha). If silviculture were standardised and plantations scrupulously followed (second scenario), this carbon stock would be better optimised in G. robusta plantations (ca. 206 Mg/ha) than in E. saligna plantations (ca. 126 Mg/ha). These carbon estimates are within the range of stocks (12-228 Mg C/ha) reported in tropical agroforestry systems [102]. Long-term CO 2 fixation was estimated to be between 50 and 65 Mg/ha in E. saligna plantations, with higher values in Kirumba, and between 30 and 100 Mg/ha in G. robusta plantations, with higher values in Sake (Table 6). While E. saligna is the species to be favoured for productivity reasons in wood-energy production plantations to try to meet the demand of the growing and disadvantaged population living in the vicinity of the PNVi, scrupulously monitored G. robusta plantations could be more interesting in terms of carbon credits. Thus, to optimise simultaneously the production of wood-energy and the carbon stock in the plantations initiated in North Kivu, E. saligna and G. robusta should be planted in a mixture. Indeed, mixed plantations are thought to be more productive, able to store more carbon, better for restoring biodiversity, and reduce the impact of disturbances more quickly than single-species plantations [11,89,103].

Conclusions and Research Perspectives
The E. saligna and G. robusta plantations in North Kivu should not be considered as mere providers of wood-energy but also as real carbon sinks, providing the areas of forest lost through logging would be compensated by new, well-managed plantations. The assessment of productivity showed differences that could lead planters to favour E. saligna. The performance of this species, including on slopes and acidic soils (Kirumba), is an asset for valorising deforested land on less favourable or uncultivated soils in the region. Estimates of carbon stock and long-term CO 2 fixation revealed that G. robusta is also a species of choice. Thus, to promote more resilient and diversified production systems, E. saligna and G. robusta should be planted in mixture. Empirical models have been developed locally to predict the productivity of monospecific even-aged plantations of the two most commonly used species in the region. The validity of these models should be extended to other growing environments and plantation types. The impact of stand density (and therefore spacing) and site characteristics on the growth and productivity of the plantations studied led us to propose pragmatic approaches for silviculture of the two species in the two target sites, including the addition of natural fertilisers where possible (itinerant breeding). The economic profitability of the plantations in North Kivu should be assessed so as to deduce their socio-economic impact in the region, including possible income from "carbon credits". To do this, rapid and efficient methods of inventorying plantations should be developed, using, for example, "remote" (satellites, drones) and/or innovative (machine learning, LIDAR) approaches. Finally, a multi-criteria analysis of their management system/mode should be carried out to discuss their sustainability.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/f13091508/s1. Table S1. Presentation of the plantations studied. ID: plantation identifier; Es: Eucalyptus saligna; Gr: Grevillea robusta; Alt.: altitude; Topo: topographical situation; Temp.: temperature; Af: equatorial climate; Am: tropical monsoon climate. Table S2. Physico-chemical characteristics of soils in the plantations studied. Es: Eucalyptus saligna; Gr: Grevillea robusta. Figure S1. Hierarchical classification of plantations based on average age (in years) and species. The average age was estimated over three years and is equivalent, for each plantation, to the arithmetic mean of the ages calculated in 2018, 2019 and 2020. Es: Eucalyptus saligna; Gr: Grevillea robusta. Figure S2. Projection of observed and predicted increases for E. saligna. Observed: observed values; Predicted1: predicted values from stand age and density; Predicted2: predicted