Annemarijn Nijmeijer

Pierre-Eric Lauri

Jean-Michel Harmand

Gregoire T Freschet

Jean-Daniel Essobo Nieboukaho

Patrick Kenfack Fogang

Seguy Enock

Stéphane Saj
email: stephane.saj@cirad.fr

Jean- Daniel Essobo Nieboukaho

Long-term dynamics of cocoa agroforestry systems established on lands previously occupied by savannah or forests

Keywords: afforestation, ecosystem service, land-use legacies, litterfall, plant diversity, soil quality, Theobroma cacao

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Ecosystem services and underlying functions depend on both current ecosystem ecological attributes and historical legacies. Disturbances, natural or anthropogenic, can lead to long-term fluctuations in ecosystem structure and functioning. The frequency, type, size, timing and intensity of a disturbance determines its impact which could eventually lead to land-use change [START_REF] Chapin | Principles of Terrestrial Ecosystem Ecology Springer[END_REF]. In general, the impact of human interference, for example in the form of agricultural activities, is larger than that of natural events (Foster et al., 2003). All organic components, above-and below-ground, seem to be highly influenced by past land-use [START_REF] Bellemare | Legacies of the agricultural past in the forested present: an assessment of historical land-use effects on rich mesic forests[END_REF][START_REF] Jangid | Land-use history has a stronger impact on soil microbial community composition than aboveground vegetation and soil properties[END_REF][START_REF] Perring | Global environmental change effects on ecosystems: the importance of land-use legacies[END_REF]. Such legacies are demonstrated for an array of ecosystem attributes such as nutrient cycling (Dupouey et al., 2002), carbon storage [START_REF] Freschet | Aboveground and belowground legacies of native Sami land use on boreal forest in northern Sweden 100 years after abandonment[END_REF], soil microbial community and heterotrophic respiration [START_REF] Kallenbach | Land-use legacies regulate decomposition dynamics following bioenergy crop conversion[END_REF]. Yet, land-use legacies might be overruled if current disturbances are strong and long enough, ultimately leading to novel ecosystems [START_REF] Hobbs | Novel ecosystems: theoretical and management aspects of the new ecological world order[END_REF].

Agricultural systems can be established from contrasting ecosystems. The conversion of a tropical forest into an agricultural suitable environment is mainly achieved by slash-and-burn technics, which are very similar all over the world [START_REF] Nepstad | Large-scale impoverishment of Amazonian forests by logging and fire[END_REF][START_REF] Achard | Determination of Deforestation Rates of the World's Humid Tropical Forests[END_REF][START_REF] Palm | Slash-and-burn agriculture: the search for alternatives[END_REF].

Tropical grasslands are also used for agricultural purposes after vegetation burning [START_REF] Kugbe | Annual vegetation burns across the northern savanna region of Ghana: period of occurrence, area burns, nutrient losses and emissions[END_REF]. Over the last century, half of the tropical grasslands and savannahs had been converted to agricultural land (UNDP, 2005). Often, the same crops are grown on differing past land-uses. This is the case of cocoa agroforestry systems (cAFS) from the forest-savannah transition zone of Central Cameroon which are either created on forest or on savannah land [START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF]. In the forestderived cAFS (F-cAFS), farmers operate a selective clearing of forest trees for the provision of shade along with food and/or additional income [START_REF] Jagoret | Farmers' assessment of the use value of agrobiodiversity in complex cocoa agroforestry systems in central Cameroon[END_REF]Saj et al., 2017a). In the savannahderived cAFS (S-cAFS), farmers use specific techniques to build up a tree canopy [START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF].

Central Cameroon farmers' practices are known to drive cAFS' structure, tree species composition and biomass on the long run [START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF]Saj et al., 2017a;Jagoret et al., 2018a).

Interestingly, full-grown S-cAFS and F-cAFS seem to exhibit comparable multi-strata structure, cocoa yields and C storage abilities despite differing previous land-use and a fortiori differing legacies [START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF][START_REF] Saj | Carbon storage and density dynamics of associated trees in three contrasting Theobroma cacao agroforests of Central Cameroon[END_REF][START_REF] Nijmeijer | Carbon dynamics in cocoa agroforestry systems in Central Cameroon: afforestation of savannah as a sequestration opportunity[END_REF].

Land-use change from savannah or forest into cAFS probably affects both above-and below-ground ecosystem functioning and related services. For instance, conversion of forest into cAFS was shown to alter litter inputs and soil properties, putatively altering nutrient cycling [START_REF] Beer | Shade management in coffee and cacao plantations[END_REF][START_REF] Schroth | Plant-soil interactions in multistrata agroforestry in the humid tropics[END_REF][START_REF] Hartemink | Nutrient stocks, nutrient cycling, and soil changes in cocoa ecosystems: A review[END_REF][START_REF] Adeniyi | Development of a composite soil degradation assessment index for cocoa agroecosystems in southwestern Nigeria[END_REF]. Similarly, transitioning from savannah to another (agro or eco)system or vice-versa was shown to alter both litter inputs and soil properties (Don et al., 2011;[START_REF] Nouvellon | Age-related changes in litter inputs explain annual trends in soil CO2 effluxes over a full Eucalyptus rotation after afforestation of a tropical savannah[END_REF][START_REF] Sugihara | Effect of vegetation on soil C, N, P and other minerals in Oxisols at the forest-savanna transition zone of central Africa[END_REF]. In this paper we explored the long-term legacies of past land-use in F-and S-cAFS from Central Cameroon.

First, we hypothesised that S-and F-cAFS ecosystem attributes were likely to progressively drift away from those of their initial ecosystem. We expected that the timespan needed to observe significant changes would depend both on the previous land-use (PLU) and the considered ecosystem attribute. Secondly, we hypothesised that the magnitude and the direction of the temporal evolution, hereafter called trajectories, of some ecosystem attributes were to differ between S-and F-cAFS. We expected these differential trajectories to attain comparable levels for some of the studied attributes. We tested both hypotheses using a selected number of attributes in cAFS, whose levels were studied through decade-long chronosequences. On the one hand, we focused on S-cAFS trajectories and their PLU, savannah, and on the other hand, we studied F-cAFS trajectories and their PLU, forest. Each of the 25 ecosystem attributes studied (plant biomass, production, litter production and storage, soil organic matter and nutrient contents and associated perennials diversity…) relates directly to at least one of the following six (dis)services: species conservation, carbon storage, crop production, nutrient cycling, soil quality and soil pollution.

Material and methods

Site characteristics and sampled plots

The study was carried out in the district of Bokito, in the villages Bakoa and Guéfigué (4⁰30 latitude N and 11⁰10 longitude E), located in a forest-savannah transition zone in Central Cameroon. The landscape consists of hills with gentle slopes at an altitude between 400 and 550 m a.s.l. and is characterised by a patchwork of forests, agroforests and herbaceous savannahs [START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF].

Annual rainfall ranges from 1300 to 1400 mm with a main dry season lasting from mid-November to the beginning of March [START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF]. The region is dominated by desaturated ferralitic soils (Elangwe, 1979). Here the cocoa plantations consist mainly of decade-old diversified agroforestry systems (cAFS) containing many associated trees [START_REF] Sonwa | Diversity of plants in cocoa agroforests in the humid forest zone of Southern Cameroon[END_REF][START_REF] Jagoret | Long-term dynamics of cocoa agroforests: a case study in central Cameroon[END_REF][START_REF] Saj | Carbon storage and density dynamics of associated trees in three contrasting Theobroma cacao agroforests of Central Cameroon[END_REF]. We selected 16 plots of cAFS created on savannah (S-cAFS) and 16 plots created on secondary forest (F-cAFS) whose age was distributed along a gradient from 1 (farmers just started the planting process of cocoa) to over 70 years old. Besides, five control plots of each past land-use (savannah and forest) were selected for comparison with the S-and F-cAFS. Savannahs in the area are annually burned and periodically cultivated, therefore we chose plots where no agricultural activity could be noticed for the last seven years, as savannah controls. The forest patches are degraded compared to other secondary forests that can be found in the region (Saj et al., 2017a). Plot size was of 40 m x 60 m (2400 m²), each containing a subplot of 20 m x 40 m (800 m 2 ; see below and [START_REF] Nijmeijer | Carbon dynamics in cocoa agroforestry systems in Central Cameroon: afforestation of savannah as a sequestration opportunity[END_REF].

Ecosystem attributes

Conservation

A floristic survey enabled us to compute two attributes which relate to species conservation, ie. perennial species richness and diversity. We also computed species rarefaction curves. We surveyed all the perennial plants found in the studied systems including both typical trees and "tree-like" perennials, i.e.: cocoa, associated trees, palms and bananas. All the perennials with a diameter at breast height (DBH) over 30 cm were surveyed in 2400 m² plots. Those with a DBH between 5 and 30 cm were surveyed in 800 m 2 subplots within the 2400 m² area already selected. Height, DBH and species were recorded as described in [START_REF] Nijmeijer | Carbon dynamics in cocoa agroforestry systems in Central Cameroon: afforestation of savannah as a sequestration opportunity[END_REF]. When the DBH couldn't be measured, we followed the recommendations from Weyerhaeuser and Tennigkeit (2000) to indirectly estimate it. In total, we counted 3036 cocoa trees and 670 associated perennials in the cAFS studied. In the forest and savannah control plots, we identified 304 and 92 perennials, respectively. Associated perennials were classified according to their species: (i) native species from the African continent and (ii) exotic species introduced from other continents. Among all the individuals counted in cAFS, forest and savannah plots, 88.4 % could be identified at the species level. The number of individuals per species in the plots were extrapolated to number of individuals per species per hectare. Abundance and diversity data of associated perennials were used to calculate species richness using Shannon-Wiener and Sørensen indices [START_REF] Peet | The Measurement of Species Diversity[END_REF][START_REF] Magurran | Measuring biological diversity[END_REF].

Carbon storage

Using the above-cited survey we calculated the total biomass of large associated trees (DBH > 30 cm) and the biomass of cocoa trees, two attributes directly relating to carbon storage (Saj et al;[START_REF] Colwell | EstimateS: Statistical estimation of species richness and shared species from samples. Version 9 and earlier[END_REF].

Aboveground biomass (AGB) of associated trees was estimated after [START_REF] Chave | Improved allometric models to estimate the aboveground biomass of tropical trees[END_REF]. Their belowground biomass (BGB) was estimated after [START_REF] Cairns | Root biomass allocation in the world's upland forests[END_REF]. The biomass of each tree was then calculated as the sum of AGB and BGB.

Crop production

Five attributes related to crop production were computed using the above-cited floristic survey and cocoa pod production counts, namely: cocoa tree stand density and basal area, cocoa tree basal area share (CTBAS), accessible cocoa yield and banana stem density. CTBAS was calculated as: [basal area of cocoa trees (at 1.30m height) : basal area of whole stand] (Saj et al. (2017b). Cocoa pod production and subsequent yield in cocoa beans were estimated for three consecutive production cycles in 2014, 2015 and 2016. Pods with a length > 10 cm were counted four to five times during the year (April, June, August, October and December) to get most of the annual production cycle. Such a count indicates the maximum number of pods per plot that can be produced and physiologically reach maturity within a production cycle. This has been defined as the "accessible" production of pods from which we can then derive the "accessible" yield in cocoa beans used in this study with calculations described in Saj et al. (2017b).

Nutrient cycling

Six attributes related to nutrient cycling were computed, namely: annual leaf and total litterfalls, leaf and total standing litters, leaf litter and total litter cycling indicators. Litterfall was monitored from November 2015 to October 2016. Three 0.5 m² collectors were randomly placed within the 2400 m² area of each plot, and once a month their positions were changed to account for spatial variability (12 different positions for each collector during the study; [START_REF] Schroth | Decomposition and nutrient supply from biomass[END_REF]. They were emptied every two weeks (111 collections in total). Once collected, the litter was air-dried and sorted to distinguish leaf litter from total (ie. leaf + branches + reproductive organs) litter. Subsamples were then dried at 60 °C for 72 hours and weighed for dry matter content calculation. Standing litter was measured at the end of the rainy season (October 2016), when decomposition started to slow down due to seasonal climate transition. In each plot, litter was collected in four randomly distributed 1 m 2 quadrants for fresh weight measurement. Litters from a 0.0625 m 2 (25 x 25 cm) sub-quadrant were then dried at 60 °C for 72 hours and weighed for dry matter content calculation. These subsamples were also sorted to distinguish leaf litterfall from total (ie. leaf + branches + reproductive organs) litterfall. Standing litter and litterfall were not measured in the savannah control plots due to the very low density of perennials in these plots. Two indicators of litter decomposition were calculated after [START_REF] Tripathi | Productivity and nutrient cycling in recently harvested and mature bamboo savannas in the dry tropics[END_REF] : (i) one including all the organs constituting the litter (leaves, branches and reproductive organs), which we called "total litter cycling": [total litterfall : (total litterfall + total standing litter)] ; (ii) one including only leaves, hereafter called "leaf litter cycling":

[leaf litterfall : (leaf litterfall + leaf standing litter)].

Soil quality and contamination

Ten attributes related to soil quality and possible contamination by pesticide-driven Cu were computed, namely: pH, organic carbon (C), total N, C:N, Inorganic P, exchangeable K, cation exchange capacity, as well as exchangeable Ca, Mg and Cu. A composite soil sample of 0-15 cm depth was built from eight subsamples (each of 135 cm 3 ) systematically distributed in each 800 m 2 subplot after brushing off the litter layer. Particle size distribution and soil chemical composition and particle size distribution were determined by the IITA soil laboratory in Yaoundé (November 2016; www.iita.org).

Soils were air-dried and sieved at 2 mm. Particle size distribution (three fractions) of sieved soil samples (0-2mm) was estimated after [START_REF] Bouyoucos | A recalibration of the hydrometer method for making mechanical analysis of soils[END_REF]; pH in water was determined in a 1:2.5 (10w:v) soil:water suspension. To determine soil organic C (SOC), soils were further ground to reach a particle size ≤ 0.5 mm before proceeding to chromic acid digestion and spectrophotometric analysis [START_REF] Heanes | Determination of total organic-C in soils by an improved chromic acid digestion and spectrophotometric procedure[END_REF]. Total nitrogen (N) was determined from a wet acid digest [START_REF] Buondonno | Comparing tests for soil fertility. II. The hydrogen peroxide/sulfuric acid treatment as an alternative to the copper/selenium catalyzed digestion process for routine determination of soil nitrogen-kjeldahl[END_REF] and analyzed by colorimetry [START_REF] Anderson | Tropical soil biology and fertility: A handbook of methods[END_REF]. Inorganic phosphorus (P) was extracted using Bray 1 extractant and analyzed using the molybdate blue procedure [START_REF] Murphy | A modified single solution method for the determination of phosphate in natural waters[END_REF].

Exchangeable cations (calcium (Ca), magnesium (Mg), potassium (K)) were extracted by ammonium acetate at pH 7 and analyzed by flame atomic absorption spectrophotometry [START_REF] David | The determination of exchangeable sodium, potassium, calcium and magnesium in soils by atomic-absorption spectrophotometry[END_REF].

Exchangeable copper (Cu) was extracted using the Mehlich-3 procedure and determined by atomic absorption spectrophotometry [START_REF] Mehlich | Mehlich 3 soil test extractant: A modification of Mehlich 2 extractant[END_REF].

Chronosequence validation

We defined four cAFS age categories [0-14], [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30], [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50] and [>50] years to represent the different temporal stages of cAFS while trying to balance the number of plots per category. Yet we noticed a significant soil clay content heterogeneity among the plots sampled [START_REF] Nijmeijer | Carbon dynamics in cocoa agroforestry systems in Central Cameroon: afforestation of savannah as a sequestration opportunity[END_REF] which challenged the validity of the chronosequence [START_REF] Pickett | Space-for-time substitution as an alternative to long-term Studies[END_REF]. To take into account a putative soil texture effect, a preliminary ANCOVA, including age categories, past land-use and soil clay content as a covariate, was run on the whole set of cAFS. It showed that the soil clay content still significantly impacted -and interacted with -cAFS age and/or previous land-use on numerous ecosystem attributes (data not shown). This was probably due to the fact that most of the S-cAFS in the sampled region were located at the foot of gentle slopes while F-cAFS were located slightly higher [START_REF] Nijmeijer | Legacies of past land-use in complex cocoa agroforestry systems in Bokito (Central Cameroon): long-term effects on ecosystem multifunctionality[END_REF]. Consequently, we split the studied plots into two according to their clay content: a group with low clay content soils (LCCS; containing between 9.8 to 16.2 % of clay), a group with high clay content soils (HCCS; containing more than 19.6 % of clay). The number of F-cAFS plots in the LCCS group and the number of S-cAFS plots in the HCCS group were too low to establish chronosequences. Hence, we ended up studying a S-cAFS chronosequence on LCCS, and a F-cAFS chronosequence on HCCS. The LCCS group was comprised of 20 plots (n S-AFS = 11, n savannah controls = 3, n F-cAFS = 3, n forest controls = 3). The HCCS group contained 18 plots (n S-AFS = 4, n savannah controls = 2, n F-cAFS = 10, n forest controls = 2). We further used ANOVAs to check, separately for each soil clay content group, that cAFS age categories did not differ in soil clay content.

Appendix 1 shows the mean values resulting from ANOVA and post-hoc (SNK) of the 25 studied attributes for the 4 land-use types regardless of plot age and soil texture (42 plots in total).

Data analyses

We performed all statistical analyses for the LCCS and HCCS groups separately. Two principal component analyses (PCA) followed by a Varimax rotation were performed on the whole set of ecosystem attributes studied to estimate cAFS plot positions in multivariate space of ecosystem attributes and identify major axes of variations. The scores of each S-or F-cAFS plot on the two axes of the PCA were then extracted and, for each axis, one-way ANOVAs with post-hoc SNK tests were run to determine which age categories significantly differed from their respective control (ie. past land-use). Each of the 25 ecosystem attributes (Tab 1) studied were then tested separately to check for differences between age categories and control plots with one-way ANOVAs followed by SNK as post-hoc (or with non-parametric Kruskal-Wallis tests followed by the Steel-Dwass-Critchlow-Fligner test when variables did not exhibit homogeneity of variance). These analyses were followed by the study of the temporal evolution of each attribute using linear, polynomial or logarithmic regressions -choosing for each the best fitted model (RMSE). The significance of each regression was then tested comparing the model's residuals against their predicted values (F-test). All statistical analyses were performed using XLStat (Addinsoft, 2015) and statistical significance was set at P < 0.05. Finally, (associated) perennial species rarefaction curves were computed to check for differences between the four land-use types. These curves were computed using EstimateS (version 9.1.0) from [START_REF] Colwell | EstimateS: Statistical estimation of species richness and shared species from samples. Version 9 and earlier[END_REF] and their visual comparison was used to evaluate the differences between land-uses [START_REF] Barlow | Quantifying the biodiversity value of tropical primary, secondary, and plantation forests[END_REF].

Results

Principal component analyses

The first two axes of the PCA explained 63% and 54.7% of the total variance of the two systems studied i.e. S-cAFS in the low clay content soil group (LCCS) and F-cAFS in the high clay content soil group (HCCS), respectively (Fig 1a,b). For both groups, the attributes related to soil quality, nutrient cycling and carbon storage were the main contributors to the first axis (D1), accounting respectively for 40.27% of the total variance in LCCS group and 36.58% in HCCS group. For both groups, the attributes related to cocoa stands were the main contributors to the second axis (D2), accounting respectively for 22.73% of the total variance in the LCCS group and 18.12% in the HCCS group (Fig

Ecosystem services levels and trajectories

Conservation

S-cAFS exhibited equivalent associated perennials diversity and species richness to F-cAFS and savannah but lower than those of forest (Tab 2a). While no significant differences in Shannon index and species richness were observed between local forest and young F-cAFS (<30 years), mature and old F-cAFS (more than 30 years old) exhibited lower associated perennials diversity values (Tab 2b).

Regressions showed that, around 30 years old, S-cAFS and F-cAFS reached comparable diversity levels (Fig 3a). S-cAFS associated perennial species richness appeared constant (ca. 7.8 sp ha -1 ) while it decreased with time in F-cAFS reaching that of S-cAFS ca. 90-100 years after establishment (Fig 3b).

Rarefaction curves were computed for African and total (associated) perennial species identified in each system type (Fig 4a,b). Forested systems (forest control, F-and S-cAFS) revealed comparable curve shapes while their levels underlined clear differences. These systems were found to exhibit very different species accumulation abilities than those exhibited by savannah (Fig 4 a,b).

Furthermore, the sampling effort needed in cAFS to reach forest control species richness levels appeared at least three times higher. S-cAFS exhibited slightly lower abilities to accumulate species than F-cAFS -notably due to their lower associated tree abundance [START_REF] Nijmeijer | Carbon dynamics in cocoa agroforestry systems in Central Cameroon: afforestation of savannah as a sequestration opportunity[END_REF]

Carbon storage

In S-cAFS, large (DBH >30cm) associated tree and cocoa tree biomass increased with age, especially between [15-30] and [31-50] years old categories (Tab 2a). Large associated tree biomass did not differ significantly between F-cAFS age categories, S-cAFS and forest plots (Tab 2b). Regressions showed that large associated tree biomass in S-cAFS reached about the same values (140 t ha -1 ) as F-cAFS ca. 70 years after establishment (Fig 3c). Regressions showed that cocoa tree biomass followed similar trends and reached ca. 20-25 t ha -1 50 years after establishment (Fig 3d).

Crop production

For the attributes related to cocoa production, no significant differences were found between F-cAFS age categories and S-cAFS (Tab 2b). For S-cAFS, basal area of cocoa trees, cocoa tree basal area share (CTBAS) and accessible yield were found to increase with time, especially between [0-14] and [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30] years old categories (Tab 2a). Regressions underlined that S-cAFS exhibited more cocoa trees during the first 20 years (> 1500 ha -1 ) (Fig 5a). Cocoa tree density reached a maximum when S-cAFS were around 20 years old and when F-cAFS were approx. 40 to 50 years old. They then slowly decreased and stabilised around 1000 trees ha -1 for both systems (Fig 5a). Cocoa tree basal area regressions followed the same trend for S-cAFS and F-cAFS up until 60-70 years of age. Afterwards, basal areas decreased in F-cAFS (Fig 5b). Cocoa tree basal area share (CTBAS) followed the same trend in S-and F-cAFS while it appeared consistently higher in S-cAFS. After 60 years, regressions showed that CTBAS reached a value of ca. 40% in S-cAFS versus ca. 30-32% in F-cAFS (Fig 5c). F-cAFS accessible yields did not evolve significantly with time and showed a mean of 693 kg ha -1 . S-cAFS accessible yield increased with time and tended to reach values of 800 kg ha -1 or higher 40 years after establishment (Fig 5d). Banana stems appeared more numerous in F-cAFS that were [15-30] years of age (Tab 2b).

Banana density did not however show any significant trend on the long-term (mean of 97.5 ind ha -1 ) in F-cAFS. Contrastingly, it quickly decreased after S-cAFS establishment as banana stems were then barely found in plots older than 5 years of age (Fig 5e).

Nutrient cycling

In S-cAFS, leaf litterfall increased the first 20 years after establishment and could not be distinguished from F-cAFS and forest levels from 20-30 years onwards (Tab 2a; Fig 6a). S-cAFS total litterfall increased with time but remained lower than that of local forest (Tab 2a,Fig 6d). In F-cAFS, leaf litterfall did not significantly evolve with time and showed a mean of 4.88 t year -1 ha -1 . Total litterfall was found lower than that of forest in old F-cAFS (> 50 years) (Tab 2b). Total litterfall increased significantly until ca. 40 years in S-cAFS, time at which it reached the F-cAFS total litterfall mean (7.93 t year -1 ha -1 ; Fig 6d). No significant differences could be detected for standing litters for both S and F-cAFS (Tab 2a, b). Yet, in F-cAFS, regressions showed that leaf standing litter increased with time and could reach values close to 2 t ha -1 for plots older than 80 years -a value corresponding to the mean obtained in S-cAFS at the same age (Fig 6b). In S-cAFS, total standing litter increased marginally with time (p < 0.1) to reach equivalent amounts to those in F-cAFS already 15-20 years after establishment (mean F-cAFS = 4.19 t ha -1 ; Fig 6e). In S-cAFS, total litter cycling significantly increased and could not be distinguished from F-CAFS or forest total litter cycling when over 30 years of age (Tab 2a). For F-cAFS, both litter cycling indicators decreased with time and ended up lower than forest indicators, after 30 years for leaf cycling and 50 years for total cycling (Tab 2b). Regressions showed that both leaf litter and total litter cycling indicators displayed significant opposite trends between S-cAFS (increasing) and F-cAFS (decreasing) to eventually reach comparable values at ca. 60 years of age (Fig 6c,e).

Soil quality and contamination

In S-cAFS, soil pH, cation exchange capacity (CEC) as well as soil organic C, total N, exchangeable Ca and Cu concentrations increased significantly with time -reaching levels comparable to F-cAFS and/or local forest after 30 or 50 years. Oppositely, the C:N ratio significantly decreased with time (Tab 2a). In F-cAFS, soil attributes didn't reveal much significant differences with time, inorganic P was found lower after 50 years old, C:N ratio tended to increase with time while pH exhibited an opposite trend (Tab 2b). Regressions showed that S-cAFS soil pH increased over time to reach that of F-cAFS (ca. 6.1-6.2) 50 years after establishment (Fig 7a). No significant evolution was found for soil organic C in F-cAFS, these values had large variability (mean = 2%). Soil organic C increased significantly in S-cAFS where the regression showed a ca. 50% increase 70 to 80 years after plot establishment (Fig 7b). Regressions showed that 65-70 years after establishment, S-cAFS and F-cAFS exhibited similar levels of soil total N i.e. approx. 0.13% (Fig 7c). Regressions for soil C:N ratio showed an increase for F-cAFS and a steady decrease for S-cAFS -similar levels (ca. 14) being reached around 50-55 years after plot establishment (Fig 7d). No significant regression was found for soil exchangeable P in S-cAFS, however a marginal decrease was detected in F-cAFS (mean P in S-cAFS= 4.8 ppm; Fig 7e). Regressions made on CEC were not significant in F-cAFS (mean= 7 cmol kg -1 ) whilst a significant increase was found in S-cAFS (Fig 7f). No significant regression was found for soil exchangeable K content in S-cAFS (mean= 0.16 cmol kg -1 ) while the decrease was significant in F-cAFS (Fig 7g). Regressions made on exchangeable Ca showed a marginal decrease in F-cAFS and a significant increase in S-cAFS -similar levels (ca. 5-5.5 cmol kg -1 ) being met around 65 years after plot establishment (Fig 7h). Regressions showed no significant evolution for exchangeable Mg, neither in S-cAFS nor in F-cAFS (mean: S-cAFS= 1 cmol kg -1 ; F-cAFS= 1.9 cmol kg -1 ; Fig 7i). Finally, regressions made on soil exchangeable Cu content showed a significant increase in S-cAFS while no trend could be established for F-cAFS (mean = 7 ppm, Fig 7j).

Discussion

S-cAFS and F-cAFS general trajectories

Our results demonstrated that F-and S-cAFS underwent highly contrasting trajectories after their establishment. Ecosystem attributes related to soil quality, plant diversity and nutrient cycling appeared to generally converge among S-cAFS and F-cAFS along the chronosequences studied, whereas attributes related to crop production varied in a more similar way. Overall our results revealed that, for the pool of attributes and services studied, agroforestry systems always evolved away from their past land-use state. Such differences became statistically apparent after 15 to 30 years for S-cAFS and after 30 to 50 years for F-cAFS. Hence, depending on attributes, past land-use legacies differed in magnitude and lasted from a few years to several decades. S-and F-cAFS displayed comparable levels of most attributes after a certain timespan, as previously suggested by other studies in the same region [START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF]Saj et al., 2017a;[START_REF] Nijmeijer | Carbon dynamics in cocoa agroforestry systems in Central Cameroon: afforestation of savannah as a sequestration opportunity[END_REF].

However, for other attributes the temporal projections did not seem to display a stabilization phase at later stages. It suggests that these attributes, after converging at intermediate time scales, may keep evolving in contrasting directions for S-and F-cAFS. This could be related to the different types of soils studied [START_REF] Adeniyi | Development of a composite soil degradation assessment index for cocoa agroecosystems in southwestern Nigeria[END_REF] as well as farmers' differential management within a generation and between generations (Jagoret et al., 2018a) -both putatively altering ecosystem's functioning.

Answering this question would need further investigations on a broader array of soil types and longer timespans.

Ecosystem services temporal trajectories

Species conservation

Our results underlined the relatively fast (ca. 20 years) convergence in terms of perennial species diversity between S-and F-cAFS. Interestingly, S-cAFS exhibited a large number of African species. This high number of native species in S-cAFS could be explained by the important gene bank (native trees, seeds and seedlings) present in their soils or in that of nearby forests or agroforests. Such a result also points at the active introduction, preservation or transfer of these species into S-cAFS (Jagoret et al., 2018b). Interestingly, rarefaction curves highlighted the similarity of species accumulation between S-and F-cAFS. Yet, despite a fair proportion of African species (over 50% of associated individuals), S-cAFS demonstrated an overall lower potential of conservation since the abundance of associated perennials was found lower than in F-cAFS [START_REF] Nijmeijer | Carbon dynamics in cocoa agroforestry systems in Central Cameroon: afforestation of savannah as a sequestration opportunity[END_REF]. Besides, some species occurred in S-or F-cAFS even they were not present in forest controls. This was the case of locally consumed Cola species but also of Entandrophragma cylindricum, a commercial timber species, which is currently classified 'vulnerable' in the IUCN red list (IUCN, 2018). All in all, F-cAFS still exhibited a forest legacy a century after establishment whilst both S-and F-cAFS showed an ability to enrich the local pool of species. Our study emphasised the important role that both agroforestry systems may currently play in the region in perennial species conservation, even though they do not preserve as many species as forests [START_REF] De Beenhouwer | A global meta-analysis of the biodiversity and ecosystem service benefits of coffee and cacao agroforestry[END_REF][START_REF] Mortimer | Supporting and regulating ecosystem services in cacao agroforestry systems[END_REF]. Furthermore, since associated diversity and abundance are largely driven by the practices and uses of associated flora [START_REF] Jagoret | Farmers' assessment of the use value of agrobiodiversity in complex cocoa agroforestry systems in central Cameroon[END_REF]Saj et al., 2017a), a particular attention shall be given to old systems which rejuvenation depends on farmer strategies and may eventually lead to less multifunctional and simplified systems (Saj and Jagoret, 2017;Jagoret and al, 2018b).

Carbon storage

In terms of C storage, the land-use legacy appeared very significant in young systems and tended to disappear on the long term. While the biomass of large associated trees of S-cAFS increased steadily over decades, it decreased in F-cAFS. These results underline the ability of S-cAFS to gain in C storage with ageing, at least up to a certain point [START_REF] Saj | Carbon storage and density dynamics of associated trees in three contrasting Theobroma cacao agroforests of Central Cameroon[END_REF]. The results further emphasised the slow, but steady, decrease in C storage abilities in F-cAFS. In a region where old cAFS are numerous [START_REF] Jagoret | Long-term dynamics of cocoa agroforests: a case study in central Cameroon[END_REF], such results point to the need to include selection/management (at the species and individual level) of potentially large trees in the renewal process of these systems (Jagoret et al., 2018b). Such management schemes would lead to the renewal of mature and old systems whilst keeping good C storage abilities. However, the timespans -and levels -of C storage of Cameroonian cAFS remain highly significant in comparison to other agricultural systems or cocoa monocultures [START_REF] Schroth | Contribution of agroforests to landscape carbon storage[END_REF][START_REF] Schneidewind | Carbon stocks, litterfall and pruning residues in monoculture and agroforestry cacao production systems[END_REF].

Crop production

The attributes chosen to indicate cocoa production abilities correspond to cocoa stands' or system attributes known for their contribution to yield in the studied systems (Saj et al., 2017b). The higher density of cocoa trees in S-cAFS compared to F-cAFS up to 20-30 years after establishment may reflect the insurance practices of local farmers which often "overseed" their S-cAFS to control grasses such as Imperata cylindrica [START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF]. These practices probably support successful establishments since accessible yields in S-cAFS appeared overall similar or slightly higher than that of F-cAFS plots. Except for very young S-cAFS, where banana stems provide beneficial shade to young cocoa seedlings and temporary income/food/feed [START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF], banana stems appeared more numerous in F-cAFS than in older S-cAFS. This discrepancy could underline farmer practices favouring cocoa in soils where competition for water can be harsh as the S-cAFS studied were established on relatively sandy soils [START_REF] Van Vliet | Chapter Five -Mineral Nutrition of Cocoa: A Review[END_REF]. It could also underline the differential role S-and F-cAFS may have in the provision of co-products [START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF][START_REF] Saj | Carbon storage and density dynamics of associated trees in three contrasting Theobroma cacao agroforests of Central Cameroon[END_REF]. Such differentiations may be interpreted as management legacies inherited from past land-use or adaptation to soil type. Furthermore, cocoa yields reached comparable levels at the latest 15 years after cAFS establishment. This confirms that, despite unfavourable conditions in savannah (such as low soil fertility, weed competition and almost no shade trees), afforestation is clearly achievable using shaded cocoa [START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF]. Finally, the constantly lower basal area share of cocoa trees in F-cAFS compared to S-cAFS may be considered as a legacy from previous land use (forest) as mature trees are already present at establishment. It can be hypothesised that the maintenance of these trees with an already developed canopy may alter light capture on the long term [START_REF] Blaser | Shade trees have limited benefits for soil fertility in cocoa agroforests[END_REF].

Nutrient cycling

Litterfall and its subsequent decomposition are considered to be effective in the improvement and conservation of soil quality and, as such, play an important role in the regulation of nutrient cycling [START_REF] Vitousek | Nutrient cycling in moist tropical forest[END_REF]de Carvalho et al., 2014). Litterfall in S-cAFS increased steadily up to 40-50 years after establishment, reaching the levels of F-cAFS litterfall. In F-cAFS the forest legacy maintained a leaf litterfall amount close to that of the initial forest. The litter cycling indicators clearly showed an increase in S-cAFS contrasting with a decrease in F-cAFS. These indicators show comparable levels around 50 years, suggesting a convergence for both S-and F-cAFS. The temperature and moisture levels in forests with a closed canopy generally favour decomposition of litter and could explain the lower standing litter of the forest control plots compared to the more open canopy of cAFS [START_REF] Prescott | The influence of the forest canopy on nutrient cycling[END_REF][START_REF] Steffan-Dewenter | Tradeoffs between income, biodiversity, and ecosystem functioning during tropical rainforest conversion and agroforestry intensification[END_REF]. The probably drier and warmer conditions in young S-cAFS, which are still lacking a closed canopy and producing mainly low quality cocoa litter [START_REF] Dawoe | Litterfall and litter nutrient dynamics under cocoa ecosystems in lowland humid Ghana[END_REF], could also explain the higher litter stock in young S-cAFS, despite their lower litter production. Overall, the similar cocoa productivity after savannah and forest conversions may suggest that, while young cocoa plantations established after savannah benefited less from the nutrients released by litter decomposition than these established on forest, the lower competition for light they encountered offset the low soil fertility and actually authorised an equivalent productivity level compared to young F-cAFS [START_REF] Blaser | Shade trees have limited benefits for soil fertility in cocoa agroforests[END_REF]. Besides, as cAFS mature and the cocoa stands develop, the poor quality of cocoa leaf litter and its subsequent lower decomposition rate [START_REF] Dawoe | Litterfall and litter nutrient dynamics under cocoa ecosystems in lowland humid Ghana[END_REF] may account for the slight increasing leaf standing litter noticed in F-cAFS chronosequence.

Soil quality and contamination

In S-cAFS established on low clay content soils many of the studied attributes exhibited a steady evolution with time underlining an overall improvement of soil quality and nutrient availability. This improvement can be explained by the increase of litter -and linked litter cycling -in S-cAFS over time when compared to the very low soil organic content (SOC) content and litter cycling of the initial savannahs. This was not the case in F-cAFS established in high clay content soils where a decrease in pH, K and Ca and an increase of the C:N ratio may point to a low but steady decline of soil quality and nutrient availability -results in line with the recent study of [START_REF] Adeniyi | Development of a composite soil degradation assessment index for cocoa agroecosystems in southwestern Nigeria[END_REF]. Noticeably, soil Cu concentrations seemed to increase with time at least in S-cAFS and could point to a putative soil contamination by pesticides, as farmers use copper oxide or hydroxide to treat against black pod disease [START_REF] Sonwa | Management of cocoa: Constraints during acquisition and application of pesticides in the humid forest zones of southern Cameroon[END_REF][START_REF] Jagoret | Long-term dynamics of cocoa agroforests: a case study in central Cameroon[END_REF]. As a result of the differences in previous land-uses (savannah vs. forest) -emphasised by the differences in soil texture (LCCS vs HCCS) -the convergence of soil attributes between S-and F-cAFS took a long the timespan: 40-60 years after plot establishment. This timespan may have been increased by the lower soil clay content of the initial savannah. Indeed, [START_REF] Nijmeijer | Carbon dynamics in cocoa agroforestry systems in Central Cameroon: afforestation of savannah as a sequestration opportunity[END_REF] found that the annual increase in SOC concentration of S-cAFS was lower in soils with low clay content (10-15%) than in soils with higher clay content (20-25%). In the same way, Feller et al. (2001) reported that the annual rate of soil carbon increase was lower in a degraded than in a non-degraded soils. Soil organic C concentration is generally positively related to fine silt and clay content [START_REF] Plante | Impact of Soil Texture on the Distribution of Soil Organic Matter in Physical and Chemical Fractions[END_REF]. As CEC is a positive function of soil clay and organic carbon concentrations, forest soils also may have had intrinsically larger nutrient contents than savannah soils.

The values found for macronutrients fell within the range already noticed in Cameroonian cAFS (Duguma et al., 2001). Yet, soil organic matter (represented here by the SOC) did not decrease in F-cAFS after conversion contrasting with other studies [START_REF] Yang | Long-term impacts of land-use change on dynamics of tropical soil carbon and nitrogen pools[END_REF][START_REF] Snoeck | Temporal changes in VAM fungi in the cocoa agroforestry systems of central Cameroon[END_REF]. The forest control plots in our study region were highly disturbed with low and rather similar litter production as in F-cAFS, which may explain the rather similar SOC contents of both systems. Furthermore, the production of cocoa was reported to be linked to the acidification of the soil [START_REF] Hartemink | Nutrient stocks, nutrient cycling, and soil changes in cocoa ecosystems: A review[END_REF][START_REF] Snoeck | Temporal changes in VAM fungi in the cocoa agroforestry systems of central Cameroon[END_REF][START_REF] Adeniyi | Assessing the relationship between soil quality parameters of Nigerian alfisols and cocoa yield[END_REF]. However, despite a clear decrease over the long term, pH in F-cAFS did not show significant differences with the forest control. Other studies revealed idiosyncratic results [START_REF] Dawoe | Litterfall and litter nutrient dynamics under cocoa ecosystems in lowland humid Ghana[END_REF][START_REF] Isaac | Diagnosis of Nutrient Imbalances with Vector Analysis in Agroforestry Systems[END_REF] and it may be hypothesised that the composition and abundance of associated trees and their diverse litter restitution to the soil could account for such results (Finzi et al., 1998). This could also partly explain the pH increase with time in S-cAFS. The overall lower pH and nutrient contents in savannah plots might be explained by annual burning and periodical cultivation of savannah with staple crops, practices which are known to lead to extensive carbon and nutrient losses resulting in soil desaturation [START_REF] Kugbe | Annual vegetation burns across the northern savanna region of Ghana: period of occurrence, area burns, nutrient losses and emissions[END_REF]Dubiez et al., 2018). The low levels of soil C and N in savannah plots are in accordance with previous studies [START_REF] Lal | The potential of world cropland soils to sequester C and mitigate the greenhouse effect[END_REF][START_REF] Jagoret | Afforestation of savannah with cocoa agroforestry systems: a small-farmer innovation in central Cameroon[END_REF]. However, the increase of organic C in older S-cAFS showed that the systems were able to build-up SOC over time and increase their soil fertility which may support cocoa yield [START_REF] Adeniyi | Assessing the relationship between soil quality parameters of Nigerian alfisols and cocoa yield[END_REF]. The decreasing trajectory exhibited by K in the soils of F-cAFS after conversion may be linked to lower recycling of K in cAFS than in forest and high absorption of K by cocoa trees and its exportation with cocoa-beans harvest (van Vliet and Giller, 2017). Yet, in S-cAFS, the relative stability of K at a minimum level over time and the trend towards these K levels in F-cAFS after around 80 years suggests that K input to the top soil through recycling (litter and throughfall) [START_REF] Hartemink | Nutrient stocks, nutrient cycling, and soil changes in cocoa ecosystems: A review[END_REF]) could, at least partly, compensate K loss through cocoa bean exports.

Conclusion

As a whole, the 25 attributes related to the six ecosystem (dis)services (ES) studied underwent distinct temporal trajectories and showed that S-and F-cAFS could be significantly distinguished from their past land-use, savannah and forest, after 30 and 50 years, respectively. Past land-use legacies were strong enough to last decades but, in the end, S-and F-cAFS agroforests revealed typical and largely convergent profiles despite different initial systems and soil types. Such results indirectly point to the role of farmers in the establishment and maintenance of decade-old agroforestry systems as well as the significant contribution of associated plant diversity to ES provision in these systems. The levels and timespans at which the studied complex cAFS provide multiple ES are worth comparing with simpler (or full sun) systems. In the current context of climate change, where the choice of species and densities of associated perennials in cocoa plantations are debated, such multifunctionality comparisons could help determine long-term strategies which would use plant biodiversity as an actionable lever to adapt and mitigate climate change.

de Carvalho, W.R., Vasconcelos, S.S., Kato, O.R., Capela, C.J.B., Castellani, D.C., 2014. Short-term changes in the soil carbon stocks of young oil palm-based agroforestry systems in the eastern Amazon. Agroforestry Systems 88, 357-368. Don, A., Schumacher, J., Freibauer, A., 2011. Impact of tropical land-use change on soil organic carbon stocks -a meta-analysis. Glob. Change Biol. 17 Foster, D., Swanson, F., Aber, J., Burke, I., Brokaw, N., Tilman, D., Knapp, A., 2003. The Importance of Land-Use Legacies to Ecology and Conservation. BioScience 53, 77-88. Table 1 Regressions coefficients (R 2 ) of the principal component analyses (PCA) of the studied attributes after Varimax rotation for the first two dimensions, D1 and D2, run for S-cAFS in low clay content soils and for F-cAFS for high clay content soils (in bold: R² > 0.5). Tables 2a and2b Results from ANOVA (F) or Kruskal-Wallis (H) and mean values for the 25 attributes in S-and F-cAFS, savannah and forest controls, from plots with low clay content soils (LCCS, 2a) and plots with high clay content soils (HCCS, 2b). In the case of LCCS plots, age categories were applied only on S-cAFS (2a is followed by an exponent stating for which dimension (D1 and/or D2) its mean was found significantly different from past land-use mean. In (a), forests and F-cAFS plots and means appear in grey. These plots contributed to the PCA but were not included in the ANOVA made on plot scores to distinguish between S-cAFS age categories and savannah. In (b), savannah and S-cAFS plots and means appear in grey. These plots contributed to the PCA but were not included in the ANOVA made on plot scores to distinguish between F-cAFS age categories and local forest. Arrows between savannah and S-cAFS (a, orange) and between forest and F-cAFS (b, green) were added manually on the biplot to highlight the "virtual" path followed by the agroforestry systems according to their past land-use. and grey for S-cAFS) correspond to the variables' means when no significant temporal evolution could be found. °: P < 0.1; *: P < 0.05; **: P < 0.01.

Tables

Service

Figure 4

African (a) and total (b) perennial plant species rarefaction curves for forest (green line with green squares), F-cAFS (green line with white squares), S-cAFS (orange line with white circles) and savannah (orange line with orange circles) plots. The 95% intervals (shaded regions) were obtained by a bootstrap method based on 100 replications. The numbers within the brackets next to the system type correspond to: the number of plots needed (first number) to reach the sample size, the number of individuals counted (second number) and the number of species (third number) obtained for the studied sample size.

Graphical abstract caption

Highlights of the study design and considered ecosystem services.

  Fig 2a). Agroforestry systems established on forest lands (F-cAFS) could not be distinguished from

  . The share of African individuals in S-cAFS and F-cAFS appeared close, being respectively of 44.1 and 47.6% (Fig 3 a,b). Forest control plots exhibited a fair number of exotic individuals (> 10%, Fig 4 a,b). Finally, Fand S-cAFS shared 29 and 19 African species with forest plots, respectively (Sørensen index: 0.68 and 0.48). S and F-cAFS shared 22 African species (Sørensen index = 0.61).

  , 1658-1670. Dubiez, E., Freycon, V., Marien, J.-N., Peltier, R., Harmand, J.-M., 2018. Long term impact of Acacia auriculiformis woodlots growing in rotation with cassava and maize on the carbon and nutrient contents of savannah sandy soils in the humid tropics (Democratic Republic of Congo). Agroforestry Systems, online, 1-12. Duguma, B., Gockowski, J., Bakala, J., 2001. Smallholder Cacao (Theobroma cacao Linn.) cultivation in agroforestry systems of West and Central Africa: challenges and opportunities. Agroforestry Systems 51, 177-188. Dupouey, J.L., Dambrine, E., Laffite, J.D., Moares, C., 2002. Irreversible impact of past land use on forest soils and biodiversity. Ecology 83, 2978-2984. Elangwe, 1979. Carte géologique de la République du Cameroun. Ministère des mines de l'eau et de l'énergie de la République du Cameroun, Cameroon. Feller, C., Albrecht, A., Blanchart, E., Cabidoche, Y.M., Chevallier, T., Hartmann, C., Eschenbrenner, V., Larré-Larrouy, M.C., Ndandou, J.F., 2001. Soil organic carbon sequestration in tropical areas. General considerations and analysis of some edaphic determinants for Lesser Antilles soils. Nutrient Cycling in Agroecosystems 61, 19-31. Finzi, A.C., Canham, C.D., Van Breemen, N., 1998. Canopy tree-soil interactions within temperate forests: species effects on pH and cations. Ecological Applications 8, 447-454.

Figure captions Figure 1

 1 Figure captions

Figure 2

 2 Figure 2

Figure 3

 3 Figure3

Figure 5

 5 Figure 5

Figure 6

 6 Figure 6

Figure 7

 7 Figure 7

Figure 1 (

 1 Figure 1

  Figure 4 809

  Figure 1 Click here to download Figure: Nijmeijer et al -R1 AGEE 21604 Fig 1.pdf

Table 1 Click here to download Tables: Nijmeijer et al -R1 AGEE 21604 Tab 1.pdf

 1

				S-cAFS		F-cAFS
				low clay content soils	high clay content soils
	Service	Attribute	Abbreviation	D1	D2	D1	D2
	Conservation	Perennials shannon index	Shannon	0.36	0.03	0.36	0.19
		Perennials species richness (nb ha -1 )	Species R	0.28	0.13	0.31	0.05
	C Storage	Large associated trees biomass (t ha -1 )	LT biomass	0.53	0.03	0.53	0.00
		Cocoa trees biomass (t ha -1 )	Cac biomass	0.05	0.87	0.00	0.16
	Crop production	Density of cocoa trees (nb ha -1 )	Cac dens	0.00	0.70	0.02	0.72
		Basal area of cocoa trees (m 2 ha -1 )	BA Cac	0.06	0.89	0.01	0.75
		Cocoa trees basal area share (%)	CTBAS	0.01	0.76	0.01	0.25
		Mean accessible yield (kg ha -1 )	Cac yield	0.06	0.88	0.00	0.00
		Density of banana stems (nb ha -1 )	Ban dens	0.02	0.08	0.07	0.12
	Nutrient cycling	Leaf litterfall (t ha -1 )	LLF	0.56	0.11	0.04	0.08
		Total litterfall (t ha -1 )	TLF	0.59	0.15	0.23	0.09
		Leaf standing litter (t ha -1 )	LSL	0.05	0.28	0.00	0.11
		Total standing litter (t ha -1 )	TSL	0.08	0.31	0.19	0.23
		Leaf litter cycling	L cycle	0.52	0.02	0.39	0.36
		Total litter cycling	T cycle	0.70	0.05	0.56	0.46
	Soil quality	pH	pH	0.72	0.02	0.81	0.76
		Organic C (%)	Org C	0.86	0.03	0.79	0.00
		Total N (%)	N	0.80	0.01	0.92	0.02
		C:N	C/N	0.60	0.00	0.54	0.07
		Inorganic P (ppm)	P	0.64	0.00	0.63	0.01
		CEC (cmol kg -1 )	CEC	0.83	0.03	0.62	0.00
		Exchangeable K (cmol kg -1 )	K	0.34	0.14	0.45	0.00
		Exchangeable Ca (cmol kg -1 )	Ca	0.76	0.06	0.80	0.05
		Exchangeable Mg (cmol kg -1 )	Mg	0.54	0.02	0.62	0.00
	Soil contamination	Exchangeable Cu (ppm)	Cu	0.10	0.07	0.24	0.03

Table 2 Click here to download Tables: Nijmeijer et al -R1 AGEE 21604 Tab 2.pdf

 2

	(2b)
							High clay content soils
					F-cAFS	F-cAFS	F-cAFS	F-cAFS
	Service	Attribute	F or H	forest	[0-14] years	[15-30] years	[31-50] years	[50>] years	S-cAFS	savannah
	Conservation	Perennials shannon index	4.237 *	2.571 a	2.276 ab	1.681 ab	1.500 b	1.583 b	1.775 ab	1.452 b
		Perennials species richness (nb ha -1 )	6.492 **	20.5 a	14.0 ab	13.0 ab	9.7 bc	8.3 bc	8.0 bc	2.500 c
	Carbon storage	Large associated trees biomass (t ha -1 )	2.113	284.4	248.6	267.5	182.2	124.9	191.2	14.3
		Cocoa trees biomass (t ha -1 )	4.463 *	-	5.3 b	5.1 b	25.9 a	18.2 a	19.2 a	-
	Crop production	Density of cocoa trees (nb ha -1 )	3.693 *	-	870.8	937.5	1633.3	895.8	1040.6	-
		Basal area of cocoa trees (m 2 ha -1 )	3.795 *	-	2.545	2.548	10.515	8.154	8.369	-
		Cocoa trees basal area share (%)	1.140	-	18.5	6.6	36.5	30.5	30.2	-
		Mean accessible yield (kg ha -1 )	0.080	-	630.9	669.6	723.5	732.6	801.2	-
		Density of banana stems (nb ha -1 )	4.842 *	-	62.5 b	175.0 a	112.5 ab	58.3 b	18.8 b	-
	Nutrient cycling	Leaf litterfall (t ha -1 )	3.094 °7.7 a	5.4 a	7.6 a	5.5 a	5.0 a	5.4 a	-
		Total litterfall (t ha -1 )	4.220 *	15.1 a	8.9 ab	15.4 a	9.1 ab	6.8 b	10.3 ab	-
		Leaf standing litter (t ha -1 )	3.121 °0.7 a	0.9 a	1.0 a	0.9 a	2.2 a	1.0 a	-
		Total standing litter (t ha -1 )	0.370	4.5	3.5	4.6	4.1	4.8	4.8	-
		Leaf litter cycling	10.405 **	0.836 a	0.674 b	0.812 a	0.622 b	0.573 b	0.638 b	-
		Total litter cycling	6.176 **	0.947 a	0.890 a	0.928 a	0.905 a	0.732 b	0.892 a	-
	Soil quality	pH	4.872 *	6.5 a	7.1 a	7.0 a	6.3 a	5.7 ab	6.6 a	4.9 b
		Organic C (%)	1.598	2.18	2.10	2.96	1.76	1.78	2.49	1.43
		Total N (%)	2.233	0.20	0.17	0.24	0.13	0.12	0.19	0.08
		C:N	3.832 *	11.2 b	12.1 b	12.4 b	13.8 ab	15.1 ab	13.4 ab	18.8 a
		Inorganic P (ppm)	4.800 *	7.215 ab	5.716 ab	-	3.977 ab	3.332 b	9.110 a	2.270 b
		CEC (cmol kg -1 )	1.500	8.0	7.3	11.3	5.5	6.9	9.5	5.5
		Exchangeable K (cmol kg -1 )	2.246	0.29	0.33	0.50	0.19	0.19	0.34	0.18
		Exchangeable Ca (cmol kg -1 )	2.356	8.56	10.23	15.86	5.09	4.69	10.16	1.61
		Exchangeable Mg (cmol kg -1 )	1.076	2.14	2.15	2.86	1.38	1.92	2.17	1.11
	Soil contamination	Exchangeable Cu (ppm)	0.928	3.44	2.78	4.39	3.48	15.79	14.83	1.93

Acknowledgments

This study was supported by the STRADIV (n°1405-018; Agropolis Fondation) and AFS4FOOD (EuropeAid/130-741/D/ACT/ACP) projects as well as by CIRAD (French Agricultural Research Centre for International Development) and IRAD (Institute of Agricultural Research for Development of Cameroon). This research was conducted within the Research and Training Platform "DP Agroforestry

Cameroon". We thank A. Agoume and J.P. Bidias, our field assistants in Bokito, and R. Ndango (IITA).

We also thank L. Defaye for her kind revision of the English language of this manuscript.

Appendix 1

Results from ANOVA (F) and mean values of the 25 studied ecosystem attributes in S-and F-cAFS, savannah and local forest controls, irrespective to their soil clay content ad age. ° : P < 0.1; *: P < 0.05; **: P < 0.01; ***: P < 0.001. Differences at P < 0.05, after SNK test, between land-uses are indicated by different letters after the mean.


Kevin Godineau
email: kevin.godineau@ens-paris-saclay.fr

Sylvain Lavernhe

Christophe Tournier

Opto-mechanical modelling of an additive manufacturing laser scanning head including assembly defects

Keywords: Accuracy, Assembly, Defect, Geometric modelling

Additive manufacturing (AM) machines used in SLM or SLA are composed of galvanometric scanning system in order to focus and steer a laser beam to melt raw materials. A major problematic in these AM processes is to master the laser spot position which essentially depends on the machine geometry and the galvanometers angular positions. In literature, most of existing models make strong assumptions concerning the geometry of the laser scanning system. Moreover, the position of the laser spot is often obtained by interpolating a table of correspondence, experimentally determined, between the angular positions of the galvanometers and the cartesian coordinates in the working plane. All these approximations induce deviations of the laser spot compared to the desired position and affect machine performance. This article presents two kinematic models of the galvanometric laser scanning system in an AM machine: a nominal model of the system and a model with assembly defects consideration. These kinematic models, often used for machine tools and robots, are here applied to create a virtual AM machine. Thereby, the laser spot position can be simulated knowing the geometry of the machine, the possible assembly defects, and the orientation of the galvanometers. The second kinematic model is then used to extract the influence of assembly defects on the laser spot position. The work described in this paper allows us to highlight and quantify the theoretical impact of an assembly defect on the precision of the laser spot position in an AM machine.

Introduction

In order to increase machine tools and robot performances, mathematical models have been developed to express the relationships between the work space (endeffector position) and the joint space (axes position). These models, which are now largely integrated within the NC controls, are enhanced by numerous behaviours (assembly defects, deformation of the structure due to cutting forces...). Thus, it's possible to take into account all these parameters to model the machine in order to remove as many undesirable behaviours as possible during machining.

The work presented in this article focuses on the kinematic modelling of SLM and SLA additive manufacturing machines. For these processes, opto-mechanical systems used to position precisely the laser spot on the material can be divided into two categories [START_REF] Ehrmann | Optics for vector scanning[END_REF]: pre-objective laser scanning systems and post-objective laser scanning systems.

This article is only devoted to post-objective laser scanning systems. The laser beam passes through a focalization device and hits successively two rotary mirrors to be deviated towards the working plane to melt raw material. The rotary mirrors are actuated by two galvanometers and are used to precisely position the laser spot. In this system, the work space is defined by the laser spot position, and the joint space is defined by the galvanometers angular positions.

Well-known mathematical models establish the relationships between the laser spot position and the angular positions of the mirrors [START_REF] Cui | Calibration of a laser galvanometric scanning system by adapting a camera model[END_REF][START_REF] Eberle | Building a vector model representation of a two-axis laser scanhead using numerical analysis for simulation purposes[END_REF][START_REF] Li | Beam deflection and scanning by two-mirror and two-axis systems of different architectures: a unified approach[END_REF][START_REF] Xie | Correction of the image distortion for laser galvanometric scanning system[END_REF][START_REF] Manakov | A Mathematical Model and Calibration Procedure for Galvanometric Laser Scanning Systems[END_REF]. However, strong assumptions are made. For instance, the rotary axis and the reflection surface of the mirror are often mixed [START_REF] Cui | Calibration of a laser galvanometric scanning system by adapting a camera model[END_REF][START_REF] Xie | Correction of the image distortion for laser galvanometric scanning system[END_REF][START_REF] Manakov | A Mathematical Model and Calibration Procedure for Galvanometric Laser Scanning Systems[END_REF] which may simplify the geometry. Furthermore, from a control point of view, to reduce the computation time, correspondence table is preferred to determine the position setpoints of the galvanometers from the laser positions. To ensure a high precision of the machines, the geometrical deviations of the laser on the working plane have to be compensated by the command [START_REF] Delgado | Reducing field distortion for galvanometer scanning system using a vision system[END_REF]. But the identification of this compensation directly on the correspondence table can be a long and not optimal experimental process. It should be possible to take advantage of the continuous formulation of the geometry with defects to improve the building of the final correspondence table.

In this article, two kinematic models are proposed to model the actual behaviour of AM machines equipped with optomechanical systems. The novelty is to express directly the assembly errors of each subcomponent that make the laser spot position deviate from its nominal position. The aim of the study is to be able to visualize, to quantify and to understand the impact of these assembly defects onto the working plane.

The paper is then organized as following: Section 2 is dedicated to the nominal kinematic model of the full optomechanical system. Assembly defects are included inside the second model in section 3. Finally, in section 4, assembly defects and their influences on the working plane are analysed.

Nominal model

In order to define a mathematical formulation between joint space and work space, a nominal kinematic model is developed. This model describes the galvanometric scanning system operation. Parameters used in this model are named X, Y and Z for translations (A, B and C for rotations).

(

) (3) (4) (1) 2

Assumptions and parameterization

Due to particles ejection and smoke production in selective laser melting context, two windows w1 and w2 are necessary to dissociate the production environment from the laser source environment. The nominal model of the system represented in figure 1 Therefore, the laser beam is modelled by a light ray. 8. Distance between the mirror rotary axis ax (ay) and the mirror reflection surface mx (my) is named 𝑒 𝑚 and is not equal to zero. The nominal kinematic model is expressed considering the optical path and its successive transformations caused by each component. In order to establish more easily the kinematic model of the full device, the system is decomposed into local sub-models. Three sub-models are established. The first one, determine, whatever the system axis, the orientation of the light ray reflected by a mirror. The second sub-model gives the coordinates of the intersection point between the light ray and a flat surface. The last sub-model describes the deviation due to a window on light ray orientation and position. These three optical models render the forward kinematic model non-linear and do not allow us to use the work developed in the robotic system literature or machine tools literature.

The mathematical formulations are expressed with homogeneous formalism except for the three previous models. Transformation matrices are defined using extrinsic Tait-Bryan parameterization (xyz). Equation 1represents the homogeneous transformation matrix to express a vector from the 𝑅 𝑎𝑦0 frame into the 𝑅 𝑠 frame for the mirror y (figure 2 left):

𝑇 𝑠/𝑎𝑦0 = ( 1 0 0 0 0 cos (𝐴 𝑎𝑦0 ) sin (𝐴 𝑎𝑦0 ) 0 0 -sin (𝐴 𝑎𝑦0 ) cos (𝐴 𝑎𝑦0 ) 0 0 𝑌 𝑎𝑦0 𝑍 𝑎𝑦0 1
)

Forward Kinematic Model

The Forward Kinematic Model (FKM) of the system (𝑥, 𝑦) = ℎ(𝜃 𝑥 , 𝜃 𝑦 ) is determined by following the light ray from the source to the work plane. By using the three sub-models and the transition matrix previously described, the FKM is obtained. Variables 𝐾 1 , 𝐾 2 and 𝐾 3 , defined by equations 2, 3 and 4 are used to reduce the FKM equations size in equations 5 and 6. Many models of the literature are regrouped within this model. For example: the suppression of the windows (𝑒 𝑣 = 0), the suppression of the mirror thickness (𝑒 𝑚 = 0), and some assumptions about mirror positions allows this model

to represent the complex model described by [START_REF] Manakov | A Mathematical Model and Calibration Procedure for Galvanometric Laser Scanning Systems[END_REF]. Furthermore, as the nominal model takes into account many geometrical parameters, it allows us to model a wide variety of different two axis galvanometric laser scanning systems. Thanks to its low mathematical complexity, lots of real-time applications can embed such native equations. Hence, this model represents a virtual machine. Coupled with the Inverse Kinematic Model (IKM), it allows us to simulate the opto-mechanical chain of an additive manufacturing process.

Inverse Kinematic Model

The inverse kinematic model (𝜃 𝑥 , 𝜃 𝑦 ) = ℎ -1 (𝑥, 𝑦) is used to convert work space point coordinates (𝑥, 𝑦) provided by the slicer to joint space setpoints (𝜃 𝑥 , 𝜃 𝑦 ). The Inverse Kinematic Model cannot be computed analytically considering the FKM equations and the non-linearities. Therefore, a Newton Raphson method is used to determine the solutions numerically. This method significantly increases the computation time compared to the FKM, which makes it impossible to be used for real time computations for scan speed interpolation. As mentioned previously, correspondence tables based on this kinematic model between the two spaces (joint space and work space) can be favorably used [START_REF] Delgado | Reducing field distortion for galvanometer scanning system using a vision system[END_REF].

Model with regard to assembly defects

In this section, assembly defects are introduced inside the previous model in order to simulate the laser spot position knowing the geometry of the machine and the possible assembly defects.

List of assembly defects

There are numerous defects which disrupt the nominal operation of the laser scanner device: optical defects, surfaces geometrical defects, thermal defects which introduce deformations of all system parts, dynamic defects… The rest of the article focuses on one of the major issues: the assembly defects, i.e., the non-nominal assembly of subcomponents considering that each one is without defect.

The previous model assumptions remain valid except the third assumption and the sixth assumption mentioned above. Due to the parameterization of the nominal model, thirty assembly defects that significantly influence the laser spot position have been identified. Those assembly defects are collected in table 1. Each defect is noted 𝛿, the indices x, y and z (or a, b and c) correspond to the position defects along the axes x, y and z (or rotation defects around the axes x, y and z).

Integration of assembly defects

The modelling of the thirty assembly defects inside the opto-mechanical chain is represented in figure 3. Those defects are integrated in the nominal model by adding new non-linearized homogeneous transformation matrices. For example, in the nominal model, in order to pass from the mirror frame to the source frame, three transformation matrices are necessary (𝑇 𝑎𝑥/𝑚 , 𝑇 𝑎𝑥0/𝑎𝑥 , 𝑇 𝑠/𝑎𝑥0 ). For the model with defects, six transformation matrices are now necessary (𝑇 𝑚/𝑚𝑑𝑒𝑓 , 𝑇 𝑎𝑥𝑑𝑒𝑓/𝑚 , 𝑇 𝑎𝑥/𝑎𝑥𝑑𝑒𝑓 , 𝑇 𝑎𝑥0/𝑎𝑥 , 𝑇 𝑠/𝑎𝑥0 , 𝑇 𝑠𝑑𝑒𝑓/𝑠 ).

This new FKM with assembly defects consideration (𝑥, 𝑦) = ℎ 𝜹 (𝛿 𝑦 𝑠 , … , 𝛿 𝑐 𝑝 , 𝜃 𝑥 , 𝜃 𝑦 ) allows to simulate a greater number of possible behaviors of the system to represent the actual one. Moreover, the assembly defects inserted in the final model are not linearized, therefore it is possible to completely modify the geometry of the system by interpreting geometric defects as new geometrical parameters.

The IKM with defects is also solved with a Newton Raphson method. However, increasing the number of variables increases the mathematical complexity of the model leading to longer computation time of the jacobian matrices.

Influence of assembly defects

In this section, the two previous models are used to visualize and characterize the influence of assembly defects on the laser spot position in the working plane.

Visualization of assembly defects

The subtraction between the second model (model with all assembly defects) and the first model (nominal assembly model) allows us to quantify the theoretical impact of each assembly defect on the laser spot position (cf. equation 7). Indeed, by setting successively all defects parameters to zero except 𝛿 𝑖 which is analysed in the function ℎ 𝛿 𝑖 , its influence 𝐼 𝛿 𝑖 is directly obtained for all mirrors' positions.

𝐼 𝛿 𝑖 (𝜃 𝑥 , 𝜃 𝑦 ) = ℎ 𝛿 𝑖 (0, … ,0, 𝛿 𝑖 , 0, … ,0, 𝜃 𝑥 , 𝜃 𝑦 ) -ℎ(𝜃 𝑥 , 𝜃 𝑦 ) The displacement of the laser spot position induced by the laser source assembly defect 𝛿 𝑦 𝑠 can be divided into two components. The first one, represents the global displacement of the laser spot positions. This displacement is called linear displacement because it is due to the 𝑂 𝑝 point displacement for the joint configuration (𝜃 𝑥 , 𝜃 𝑦 ) = (0, 0), when the assembly defect is added. In figure 4 this global linear displacement is equal to 1 mm. The second component of the displacement represents the non-linearity of the two mirrors. The shape of this displacement is mainly due to the position where the laser spot hits the mirror x compared to the position of the mirror rotary axis x. This displacement, called non-linear displacement, is about 0.1 mm in figure 4.

Impact and characterization of each assembly defect

Figure 5 gathers the deviation magnitudes caused on the spot position by a unitary influence (𝛿 ? = 1 𝑚𝑚) for the eleven positioning defects. Each assembly defect is represented by a bar, the linear component of the displacement due to assembly defects is represented by the first value of the colorbar (white area) while the non-linear component is represented by the colorbar itself (blue to yellow area). The first bar corresponds to the assembly defect 𝛿 𝑦 𝑠 represented in figure 4. Figure 5 shows that 𝛿 𝑥 𝑝 and 𝛿 𝑦 𝑝 have only the linear displacement component. The assembly defects 𝛿 𝑦 𝑎𝑥 , 𝛿 𝑦 𝑎𝑦 𝑎𝑛𝑑 𝛿 𝑧 𝑝 are only composed by the non-linear component of the displacement. The nineteen orientation defects can also be represented in a similar figure. This figure, not represented in this article, would show that windows orientation defects have a negligible impact compared to other defects (relative influence equal to 1000). In general terms, the more the spot laser is distant to the origin 𝑂 𝑝 , the larger are the defects influences on the laser spot position and so on the machine precision.

Using the developed models, a virtual machine is created to quantify the displacement error of the laser spot. Figure 6 represents this error for an AM machine which the assembly defects include in Table 2. Assembly defects are randomly generated between ± 0.1 mm for position defects and ± 0.1 mrad for orientation defects.

Conclusion and perspectives

In order to study the impact of assembly defects on the laser spot position, two models have been created. The first model represents the nominal system. The second model introduces thirty assembly defects to represent more precisely the behavior of a real machine. These two models allow to highlight and quantify the theoretical impact of an assembly defect on the precision of the laser spot position in an AM machine.

One of the major perspectives of this study consists in identifying the values of the assembly defects from experimental measurements of the laser spot position. Hence, the use of these models would save a lot of time during the calibration phase of the AM machine.

  and figure 2 is submitted to the following assumptions: 1. The intersection of the laser beam and the work plane p is not influenced by the focalization device (Dynamic Focus Module or DFM). Therefore, only the laser beam positioning device is studied. 2. The origin of the laser source 𝑂 𝑠 is located at the centre of the DFM equivalent focus lens. 3. The laser beam hits the work plane perpendicularly in the following joint configuration: (𝜃 𝑥 , 𝜃 𝑦 ) = (0,0). The intersection point is named 𝑂 𝑝 . For this configuration mirrors rotary axes are named ax0 and ay0. 4. Flat surfaces (mirror surfaces, work plane, windows) are considered as perfect planes. 5. Windows thickness is constant and equal to 𝑒 𝑣 . Therefore, windows surfaces are parallel to each other. 6. Windows surfaces and work plane are also parallel to each other. 7. General geometrical optics properties are assumed.

Figure 1 .

 1 Figure 1. Nominal model

Figure 2 .

 2 Figure 2. Mirror parameterization -mirror y (left), mirror x (right)

  𝐾 1 = 𝑒 𝑚 + 𝑍 𝑎𝑦0 cos(𝐴 𝑎𝑦0 + 𝜃 𝑦 ) -𝑌 𝑎𝑦0 sin(𝐴 𝑎𝑦0 + 𝜃 𝑦 ) cos(𝐴 𝑎𝑦0 -𝐵 𝑎𝑥0 + 𝜃 𝑦 ) 𝐾 2 = 𝑍 𝑝 + 2𝑒 𝑣 -𝐾 1 cos(𝐵 𝑎𝑥0 ) 𝐾 3 = √cos(2𝜃 𝑥 ) 2 cos(2𝜃 𝑦 ) 2

Figure 3 .

 3 Figure 3. Representation of assembly defects in the model

Figure 4 represents

 4 Figure 4 represents 𝐼 𝛿 𝑦 𝑠 (𝜃 𝑥 , 𝜃 𝑦 ), the displacement of the laser spot for all joint configurations when a unitary assembly

Figure 4 .

 4 Figure 4. Impact of the assembly defect (𝛿 𝑦 𝑠 = 1 𝑚𝑚) on the laser spot position.

Figure 5 .

 5 Figure 5. Macroscopic characterization of each position assembly defect influence.

Figure 6 .

 6 Figure 6. Impact of 30 different assembly defects (cf.Table 2) on the laser spot position.

Table 1 .

 1 List of all assembly defects included in the kinematic models. 𝛿 𝑦 𝑎𝑥 𝛿 𝑧 𝑎𝑥 𝛿 𝑎 𝑎𝑥 𝛿 𝑏 𝑎𝑥 𝛿 𝑐 𝑎𝑥 Mirror x 𝛿 𝑧 𝑚𝑥 𝛿 𝑎 𝑚𝑥 𝛿 𝑏 𝑚𝑥 Rotary axis y 𝛿 𝑦 𝑎𝑦 𝛿 𝑧 𝑎𝑦 𝛿 𝑎 𝑎𝑦 𝛿 𝑏 𝑎𝑦 𝛿 𝑐 𝑎𝑦

		x	y	z	a	b	c
	Laser source		𝛿 𝑦 𝑠	𝛿 𝑧 𝑠		𝛿 𝑏 𝑠 𝛿 𝑐 𝑠
	Rotary axis x
	Mirror y			𝛿 𝑧 𝑚𝑦 𝛿 𝑎 𝑚𝑦 𝛿 𝑏 𝑚𝑦
	Windows 1				𝛿 𝑎 𝑤1 𝛿 𝑏 𝑤1
	Windows 2				𝛿 𝑎 𝑤2 𝛿 𝑏 𝑤2
	Work plane	𝛿 𝑥 𝑝 𝛿 𝑦 𝑝	𝛿 𝑧 𝑝	𝛿 𝑎 𝑝	𝛿 𝑏 𝑝 𝛿 𝑐 𝑝

Table 2 )

 2 on the laser spot position.

Table 2 .

 2 Values of assembly defects included in the kinematic model which generates the Figure6.𝛿 𝑧 𝑚𝑥 -0.0252 𝛿 𝑎 𝑚𝑥 -0.0621 𝛿 𝑦 𝑎𝑦

	defect mrad	defect mm
	𝛿 𝑏 𝑠	-0.0693	𝛿 𝑦 𝑠	-0.0633
	𝛿 𝑐 𝑠	-0.0654	𝛿 𝑧 𝑠	-0.0846
	𝛿 𝑎 𝑎𝑥	0.0302	𝛿 𝑦 𝑎𝑥	-0.0398
	𝛿 𝑏 𝑎𝑥	0.0635	𝛿 𝑧 𝑎𝑥	-0.0232
	𝛿 𝑐 𝑎𝑥	0.0533
				-0.0993
	𝛿 𝑏 𝑚𝑥 0.0293	𝛿 𝑧 𝑎𝑦	-0.0434
	𝛿 𝑎 𝑎𝑦	0.0277	𝛿 𝑧 𝑚𝑦 0
	𝛿 𝑏 𝑎𝑦	0.0184	𝛿 𝑥 𝑝	-0.0848
	𝛿 𝑐 𝑎𝑦	-0.0349	𝛿 𝑦 𝑝	0.0114
	𝛿 𝑎 𝑚𝑦 -0.0754	𝛿 𝑧 𝑝	-0.0452
	𝛿 𝑏 𝑚𝑦 0.0472
	𝛿 𝑎 𝑤1 -0.0687
	𝛿 𝑏 𝑤1	-0.0131
	𝛿 𝑎 𝑤2 0.0665
	𝛿 𝑏 𝑤2	-0.0280
	𝛿 𝑎 𝑝	-0.0736
	𝛿 𝑏 𝑝	0.0399
	𝛿 𝑐 𝑝	-0.0028