Rodolphe Dubois

Alexandre Eudes

Vincent Frémont

On Data Sharing Strategy for Decentralized Collaborative Visual-Inertial Simultaneous Localization And Mapping

This article introduces and evaluates two decentralized data sharing algorithms for multi-robot visualinertial simultaneous localization and mapping (VI-SLAM): Factor Sparsification for Visual-Inertial Packets (FS-VIP) and Min-K-Cover Selection for Visual-Inertial Packets (MKCS-VIP).

Both methods make robots regularly build and exchange data packets which describe the successive portions of their map, but rely on distinct paradigms. While FS-VIP builds on consistent marginalization and sparsification techniques, MKCS-VIP selects raw visual and inertial information which can best help to perform a faithful and consistent re-estimation while reducing the communication cost. Performances in terms of accuracy and communication loads are evaluated on multirobot scenarios built on both available (EUROC) and custom datasets (SOTTEVILLE).

I. INTRODUCTION A. Considered problem

The use of a fleet of robots is increasingly being considered for applications like the exploration or the inspection of large-scale areas. For such tasks, Simultaneous Localization And Mapping (SLAM) techniques are core algorithmic blocks which support higher-level decision and control algorithms. When performing SLAM, a robot relies on the measurements acquired by its embedded sensors to estimate in real time a model of its environment and its trajectory within it.

While single-robot SLAM has been extensively studied [START_REF] Cadena | Past, present, and future of simultaneous localization and mapping: Toward the robust-perception age[END_REF], its multi-robot counterpart has introduced new kinds of challenges to cope with when exchanging data between robots. Robots share knowledge by exchanging data packets designed to provide adequate information: i) to enrich the recipient's map by correctly registering the received data w.r.t. it; ii) to refine it by correlating the provided information with its underlying estimation model, such that subsequent re-estimations should enhance the accuracy of the estimates without compromising their consistency i.e. the correct assessment of the covariance matrices on the estimates.

Additionally to real-time requirements, robots must face various constraints when building packets, the most of which being communication constraints such as limited communication range and bounded bandwidth. Thus, communication policies should enhance the autonomy of each agent, and the amount of transmitted data should be as flexible as possible depending on those constraints.

B. Related works

1) Centralized multi-robot SLAM: Most of the data sharing methods that have been introduced so far are centralized: they only allow data exchanges between single agents and a central server unit. Forster et al. [START_REF] Forster | Collaborative monocular slam with multiple micro aerial vehicles[END_REF] developed a centralized monocular visual SLAM framework which assigns low-level tasks such as keyframe selection and processing to the agents while high-level tasks for Collaborative Structurefrom-Motion like mapping, loop closing and map merging are performed by a ground station. Karrer et al. [START_REF] Karrer | Cvi-slam collaborative visualinertial slam[END_REF] have recently introduced a similar system for VI-SLAM called CVI-SLAM. However, those methods often require to exchange significant amounts of data that might exceed the available communication bandwidth. This is the case when the whole set of 2D keypoints and visual descriptors of each frame is communicated. Furthermore, each agent must remain within the server's communication range, so it reduces its autonomy. Finally, centralized architectures are critically vulnerable to failures of the central server.

2) Decentralized multi-robot SLAM: Because of the previously cited drawbacks, many authors have focused on decentralized SLAM methods. Schuster et al. [START_REF] Schuster | Multi-robot 6d graph slam connecting decoupled local reference filters[END_REF] proposed a submap-based approach for visual-stereo SLAM that makes each robot share some local submaps with neighboring robots, which are then aligned via ICP point-cloud registration. Contreras et al. [START_REF] Contreras | Efficient decentralized collaborative mapping for outdoor environments[END_REF] also proposed a decentralized method for LiDAR SLAM. To cope with limited bandwidth constraints, other works make robots exchange summarized representations of their map to meet the bandwidth requirements and ease their alignment within the recipient robot's map. Such methods often rely on marginalization and sparsification techniques [START_REF] Mazuran | Nonlinear factor recovery for long-term slam[END_REF]. For instance, Cunningham et al. [START_REF] Cunningham | Ddf-sam 2.0: Consistent distributed smoothing and mapping[END_REF] proposed DDF-SAM in which one robots exchange marginalized maps over some shared variables (a.k.a. separators) like commonly observed landmarks. In a similar way, Lazaro et al. [START_REF] Lazaro | Multi-robot slam using condensed measurements[END_REF] make robots exchange condensed measurements computed through marginalization over separators. Paull et al. [START_REF] Paull | Communicationconstrained multi-auv cooperative slam[END_REF] proposed a more robust framework based on the bookkeeping of inter-robot exchanges for acoustic underwater decentralized SLAM. Nonetheless, few methods of that kind have been specifically designed for VI-SLAM.

C. Contributions

In this paper, we propose and evaluate Factor Sparsification for Visual-Inertial Packets (FSVIP) and Min-K-Cover Selection Packets (MKCS-VIP), two decentralized data sharing algorithms for VI-SLAM. The rest of this paper is organized as follows. In section II, the classical landmark-based VI-SLAM problem is formulated within a Maximum Likelihood (ML) estimation scheme. Sections III and IV respectively introduce FSVIP and MKCS-VIP, whose performances are evaluated in section V. Concluding remarks are finally proposed in the last section of the paper.

II. MONOCULAR VISUAL-INERTIAL SLAM

In VI-SLAM, each robot relies on the images acquired by its monocular camera, and on the specific acceleration and angular velocity measurements supplied by its Inertial Measurement Unit (IMU) to estimate its own trajectory and to build a map of its environment. A Front-End module first interprets those measurements to build online a probabilistic graphical model known as a factor graph, as represented in the figure [START_REF] Cadena | Past, present, and future of simultaneous localization and mapping: Toward the robust-perception age[END_REF], which encodes the measurements likelihood w.r.t. hidden variables purposely introduced, based on the observations, to model the trajectory and the environment.

Keyframe states

Landmark positions

Inertial factors

Visual factors

Loop closure factor Fig. 1: SLAM graphical probabilistic model

The trajectory is modelled as a discrete set X of keyframes. Each keyframe x k is described by the pose of the robot T WB,k ∈ SE(3) -where the subscript k denotes time t k , W the gravity-aligned world frame (whose origin coincides with the position of the first keyframe) and B the body frame -and its inertial state

ξ k = [ W v k B b ak B b ωk ]
whose components respectively denote the velocity, the accelerometer and the gyrometer biases. The environment is sparsely represented as a set L of triangulated 3D landmarks whose observations have first been tracked through the acquired images. As illustrated in the figure [START_REF] Cadena | Past, present, and future of simultaneous localization and mapping: Toward the robust-perception age[END_REF], all those variables are mutually constrained by Gaussian factors, which are stochastic constraints derived from the measurements. Each IMU measurement u k yields the following inertial factor:

p (x k+1 |x k , u k ) ∝ exp - 1 2 x k+1 f k (x k , u k ) 2 Σ -1 u
(1) where x (•) ∈ X , f k is the integrated discrete-time dynamics of the robot and Σ u the covariance of the residual propagated from the IMU noise model. The minus operator returns the residual in a tangent space of the ambient composite Lie group in which one x (•) evolves [START_REF] Solà | A micro lie theory for state estimation in robotics[END_REF]. Equally, each landmark observation z i k yields the following factor:

p z i k |T WB,k , W l i ∝ exp - 1 2 z i k -π C (T WB,k , W l i ) 2 Σ -1 v
(2) where W l i ∈ L, Σ v is the covariance of the visual measurement, and π C is the projection of the landmark to the camera frame C. Finally, the robot can refine its factor graph by closing loops. A robot closes a loop when it knowingly observes an area it has already mapped. It thus derives a relative pose constraint from 2D-3D correspondences, merges the associated landmarks and corrects the inertial drift it had accumulated in the meanwhile.

A Back-End module finally infers the map and the trajectory based on the resulting factor graph. Popular approaches rely on Maximum Likelihood (ML) estimators which solve the following graph optimization problem:

{X * ML , L * ML } = arg min X ,L -log p(Z|X , L) (3)
where Z is the set of visual and inertial measurements. p(Z|X , L) is proportional to the product of the factors defined by ( 1) and [START_REF] Forster | Collaborative monocular slam with multiple micro aerial vehicles[END_REF]. Equation (3) yields a nonlinear leastsquare optimization problem. Note that in VI-SLAM, the visual observations allow to estimate the trajectory and the 3D structure of the environment up to a scale factor which is made observable by the IMU measurements.

III. FACTOR SPARSIFICATION FOR VISUAL-INERTIAL PACKETS (FS-VIP)

The method described in this section is transposed from the one introduced in [START_REF] Paull | Communicationconstrained multi-auv cooperative slam[END_REF] for acoustic underwater SLAM, that we adapted to the specific constraints of VI-SLAM (e.g. the need to properly condense the visual-inertial factors and to provide enough visual information to spot inter-robot loop closures).

A. Overview of the proposed method

An overview of all the processes involved in FS-VIP is given in the figure [START_REF] Forster | Collaborative monocular slam with multiple micro aerial vehicles[END_REF] and detailed in the next paragraphs. While the map is being built as an output of the SLAM algorithm, three additional threads are running to i) precompute the packets that will be sent to the other robots; ii) handle the packet exchange and iii) integrate the received packets. In parallel to the mapping process, each robot precomputes the summarized representations of the successive portions of its map that it will send to the other robots. are two user-defined thresholds, then the current keyframe is selected. The aim is to select a minimal set of keyframes which covers the area observed from the sub-trajectory keyframes.

x 0 x 1 x 2 x 3 x 4 B0 l 1 B0 l 2 B0 l 3 B0 l 4 B0 l 5 (
2) Subgraph extraction: We compute a new packet as soon as a new keyframe is selected and an estimated distance d min has been run since the last packet computation. The corresponding factor graph is then extracted as showed in the figure (3a). It holds all the keyframes and observed landmarks states as well as the visual and inertial factors linked to the extracted sub-trajectory.

3) Local linearization point computation: All the keyframe and landmark states included in the subgraph are then locally inferred by solely optimizing over the extracted factors and fixing the pose of the first keyframe. We thus get a local linearization point which only depends on the information held within the extracted subgraph. Note that the landmark positions should be estimated relatively to the first keyframe, while the keyframe poses can be estimated relatively to any gravity-aligned reference frame L. Another remark is that the subgraph extraction may make some previously well-constrained landmarks badly-constrained (according to the criteria defined in [START_REF] Schneider | maplab: An open framework for research in visual-inertial mapping and localization[END_REF] 1 ). Those should be identified and removed after running the optimization. [START_REF] Paull | Communicationconstrained multi-auv cooperative slam[END_REF] which relies on nearest landmark association to compute interrobot loop closures, minimal visual information should be communicated to spot them in a visual-inertial framework. We choose to transmit local visual information on the keyframes selected in III-B.1. For each one, we communicate the n keypoints keypoints associated with the most observed landmarks in the map along with their full visual descriptor, and among those observed landmarks, we select the n landmarks (≤ n keypoints ) most observed landmarks within the 1 A landmark is considered as well-constrained if it has enough observer keyframes, if is is not too close and too far from its closest observer keyframe and it the disparity angle of its observation bearing vectors is sufficient extracted subgraph.

4) Visual information selection: Contrary to

5) Subgraph marginalization:

As a fourth step, the local estimate's information matrix is computed as the Fisher Information Matrix (FIM) of the subgraph, evaluated at the inferred local estimate:

I Z (X * ML , L * ML ) = z∈Z J z {X ,L} Σ -1 z J z {X ,L} (4)
where J z {X ,L} is the measurement Jacobian matrix of measurement z w.r.t. the estimated states. The computed information matrix is then marginalized to keep only the pose states TWI,k of the selected keyframes and the positions W li of the selected landmarks. The marginalized information matrix is the Schur's complement of the marginalized variables:

I D = I Z (T S , L S ) = I SS -I SM I -1 MM I MS (5)
where the subscripts S and M respectively denote the selected and the marginalized variables, and T and L the sets of keyframe poses and landmark positions.

6) Subgraph consistent sparsification:

As suggested by the figure (3b), the resulting information matrix I D is dense since it has been filled in with cross-correlation blocks between the remaining variables. The associated distribution N ( mS , I -1 D ) is far too complex to be communicated as it is. Therefore, we take advantage on sparsification techniques to compute a set of uncorrelated relative pose and relative position Gaussian factors (that we call virtual measurements) connecting pose and landmark states in a simpler topology (as represented in figure (3c)), and which yield the same local estimate and the closest state information matrix in terms of information metrics. The virtual measurement information matrix Ω S is computed as the solution of an optimization problem which minimizes the loss of information under a consistency constraint:

ΩS = arg min ΩS∈D+ D KL N ( mS , (J S Ω S J S ) -1 )||N ( mS , I -1 D ) subjected to I S = J S Ω S J S ≤ I D ( 6
)
where D KL denotes the Kullback-Leibler divergence, D + is the set of block diagonal positive definite matrices which match to the desired sparsified topology and I S is the state information matrix derived from the computed factors. J S is the virtual measurement Jacobian matrix:

J S = ∂ ({T Bi-1 Bi } i∈T S , { B0 l i } i∈L S ) ∂ ({T LBi } i∈T S , { B0 l i } i∈L S ) (7)
which is a sparse block-diagonal matrix whose pattern is represented in the figure [START_REF] Schuster | Multi-robot 6d graph slam connecting decoupled local reference filters[END_REF]. The right (or local) relative pose Jacobian matrices are: where Ad (•) is the adjoint matrix on SE(3). The consistency constraint ensures that the resulting sparsified distribution does not add any artificial information ie does not artificially reduces the joint entropy over any subset of variables. As proved in [START_REF] Mazuran | Nonlinear factor recovery for long-term slam[END_REF], since J S is invertible by construction, a solution to the unconstrained version of ( 6) is:

J TB i B j {TLB i ,TLB j } = -Ad -1 TB i B j I 6×6 (8)
TB i B i+1 B 0 l i TLB i+1 B 0 l i
[Ω * S ] i = ({J S I -1 D J S } i ) -1 (9)
where the subscript i denotes the i th diagonal block corresponding to the i th factor. We finally need to impose the consistency constraint from the solution of the relaxed problem.

Classical methods [START_REF] Vandenberghe | Determinant maximization with linear matrix inequality constraints[END_REF] to solve problem [START_REF] Mazuran | Nonlinear factor recovery for long-term slam[END_REF], based on the solving of Semi-Definite Programming problems by interior point methods, which model the consistency constraint with log-det barriers, are intractable in real-time. We therefore use a non-optimal method to enforce this constraint with less complexity. Let λ min be the minimum eigenvalue of I D .

Each eigenvalue λ i of the computed blocks is compared to λ min and replaced by it if it exceeds it.

C. Packet exchange

Packet exchanges are handled in a separate thread. As described in figure [START_REF] Forster | Collaborative monocular slam with multiple micro aerial vehicles[END_REF], each robot periodically checks which robots are within communication range. For each spotted robot, it checks in its packet sending history timestamp of its last successful packet sending to that robot. If new packets have been precomputed in the meanwhile, then they are sequentially sent to the other robot, and the communication history is updated as soon as a reception acknowledgement signal is received back from the recipient robot. Figure (3c) shows the content of each packet: it holds the means and the covariance matrices of the virtual measurements, plus the selected visual information (2D keypoints and descriptors) of each communicated vertex. Additionally, every measurement of the recipient robot associated to one of the selected keyframes (e.g. relative pose factors estimated from AprilTag [START_REF] Wang | AprilTag 2: Efficient and robust fiducial detection[END_REF] detection) is also transmitted.

D. Packet integration

Packet integration is handled by a third thread. Each received packet is first pre-aligned by connecting it to the previously received packet and to the robot trajectory thanks to the communicated inter-robot observations. Finally, loop closure algorithms are run in order to spot 2D-3D correspondences, using the communicated visual information on the selected keyframes. The found matches are then used to associate some received 2D keypoints and merge some communicated landmarks with the map landmarks.

E. Conclusion on the proposed method

As a summary, the FS-VIP method make robots regularly exchange data packets which summarize the successive portions of their map. It is inspired from the method developed by Paull et al. [START_REF] Paull | Communicationconstrained multi-auv cooperative slam[END_REF] for acoustic underwater SLAM which we adapted for the visual-inertial framework. The main differences we have introduced are the selection of the landmarks and keyframes for which ones local visual information is appended. In [START_REF] Paull | Communicationconstrained multi-auv cooperative slam[END_REF], landmarks are scarce and merged on a distance criteria while visual-inertial landmarks are much more numerous and should be merged based on visual matching. This also makes the sparsification process heavier. To bound its computation load, we pre-compute standard packets, while in [START_REF] Paull | Communicationconstrained multi-auv cooperative slam[END_REF], the sparsity of the underlying posegraph allows to compute them specifically for each single exchange. Finally, another difference is that in [START_REF] Paull | Communicationconstrained multi-auv cooperative slam[END_REF], each packet comes with an inter-robot relative range measurement computed from the time of flight.

IV. MIN-K-COVER BASED-SELECTION FOR VISUAL-INERTIAL PACKETS (MKCS-VIP)

A. Overview of the proposed method

The method proposed in this section is based on the selection of landmarks and raw measurements to compute data packets that will be exchanged between robots. An overview of the MKCS-VIP method is given by the figure [START_REF] Contreras | Efficient decentralized collaborative mapping for outdoor environments[END_REF]. All the process mentioned in this figure are detailed below.

B. Exchange handler thread

We make each robot R X hold a reception history (recall that it was a sending history in FS-VIP). For each robot R Y , it stores the timestamp t Y→X of the most recent information it holds from robot R Y , may it have been received directly from robot R Y or relayed by a third robot R Z . We assume that the whole data of R Y anterior to t Y→X which had to be communicated has been successfully transmitted to R X and that the clocks of both robots are synchronized.

We use the following method to keep the history up-todate: each robot R X periodically sends an update request to the other robots R Y it can reach. R Y returns the vector (t Z→Y ) RZ ∈R where R is the set of all robots. Based on this information, it checks if it can send some data about trajectory of any robot R Z (including itself) such that t Z→X ≥ t Z→Y + t min where t X→X is confounded with the current time and t min is user-defined. If it does, an exchange is triggered about the data of R Y posterior to t Z→Y . If the exchange boils down to relaying received data, then such data can be sent as it is. However, if R X must send some information about its own trajectory, it is additionally required than a minimal distance of d min has been run since the last exchange, and then the packet is computed as described below.

C. Packet computation

Contrary to FS-VIP, packets do not summarize the successive portions of the trajectory by computing some virtual measurement factors, but by selecting the best landmarks and visual-inertial factors which can best help to faithfully and consistently re-estimate it, as well as some localized visual information to spot inter-robot loop closures.

1) Inertial factors: First, we must provide some inertial information to make the scale factor observable. As a first option, we can communicate all the bunches of raw IMU measurements between the successive keyframes over the subtrajectory (6 doubles per measurement). A possibly lighter option regarding the communication load is to send some pre-integrated IMU measurements [START_REF] Forster | Imu preintegration on manifold for efficient visual-inertial maximum-a-posteriori estimation[END_REF]. Such measurements are relative orientation, velocity and position factors derived from the propagation of the raw IMU measurements between vertices x i and x j . They are respectively denoted:

∆R ij , ∆v ij , ∆p ij ∼ N (∆ Rij , ∆v ij , ∆ pij ), Σ ∆ (10)
where Σ ∆ is the covariance matrix resulting from the IMU noise propagation. Those factors are derived assuming constant accelerometer and gyrometer biases B bai and B bgi between t i and t j . The pre-integrated measurements can incorporate bias updates, projected as first order corrections via the Jacobian matrices. Therefore, each transmitted pre-integrated measurement should include the above estimates, covariance and Jacobian matrices One complete preintegrated IMU measurement needs 87 doubles. Thus, sending pre-integrated IMU measurements is relevant as soon as there are on average more than about 15 IMU measurements between successive keyframes.

2) Visual factors: The inertial measurements need to be completed with some visual factors so that the trajectory is properly constrained. We thus add the visual factors associated to observations of a selected subset of landmarks. Those landmarks are selected using a heuristic approximation the solution to a MIN-K-COVER problem in order to mutually constraint the most keyframes with the minimum of landmark observations. The aim is to select the minimal subset of landmarks such that each keyframe observes at least n mkc landmarks i.e. such that each keyframe is covered. The greedy heuristic involves to select the landmark which covers the most not fully covered keyframes until all keyframes are fully covered. At the end, some keyframes may be overcovered. Therefore, we iterate over the keyframes and remove the extra-observations associated to the most observed selected landmarks, paying attention not to make it badlyconstrained. The associated 2D keypoints are transmitted without any visual descriptor.

3) Visual information: The selected visual-inertial factors allow to re-estimate the considered sub-trajectory. As in the FS-VIP, we complete them with some visual information localized on selected keyframes. To take loop closures into account and avoid to communicate redundant visual information, we build and update a covisibility graph to support their selection. An edge is added between two keyframes if their commonly observed landmark count and ratio exceed the respective user-defined thresholds n lim c and r lim c . When a new keyframe is added to the map, the covisibility graph is updated, and the keyframe is selected if it is not adjacent to a previously selected keyframe. For each selected keyframe, we select on its visual frame the n keypoints keypoints associated to the most observed and constrained landmarks within the map, with their visual descriptor, as described in III-B.4.

4) Loop closures and inter-robot observations:

The relative pose factors derived from loop closures and interrobot observations are valuable information to share, and are included in the communicated packet if their detection timestamp is relevant w.r.t. the reception history of the recipient robot. In order to preserve the map consistency, it is crucial that their associated covariance matrices are properly estimated. A loop closure factor between two keyframes x i and x j is estimated from the spotted 2D-3D correspondences. It involves solving two Perspective-n-Point (PNP) problems within a RANSAC scheme to estimate the global poses of both keyframes w.r.t. to their spotted commonly observed landmarks and thus derive their relative pose. Note that since the 3D structure is needed to spot loop closures, each robot is the only one to be able to detect loop closures against its trajectory from received 2D visual information.

D. Packet integration

As a summary, the final packet includes all the raw (or preintegrated) inertial measurements along the communicated sub-trajectory, the relative poses between successive vertices (needed for initialization), the visual observations associated to the selected track-lines (communicated without visual descriptors), some visual information localized on some selected keyframes (with the visual descriptors) as well as the loop closure and the inter-robot observation factors spotted in the meantime. The received packet is chained to the previously received ones, connected to other trajectories through inter-robot observations and loop closures are searched between the received packet and the recipient robot's trajectory. Note that the IMU and camera models

E. Conclusion on the proposed method

Compared to the FS-VIP, MKCS-VIP make robots exchange more data but brings multiple advantages. First, the communicated packets are far less complex to build since they only involve the selection and the extraction of relevant data. A second advantage is that by sending raw measurements or factors directly derived from their propagation, the resulting packet is consistent and is not likely to be as conservative as the ones built in FS-VIP. Note also that all the communicated factors are relinearizable. Finally, a third advantage is that it is possible to relay data between robots and complete them with loop closure and inter-robot observation factors without any risk of double counting.

V. PERFORMANCE EVALUATION

A. Considered test scenarios 1) Used datasets: Both presented methods are evaluated on multi-robot scenarios built on the 5 aerial sequences from the EUROC dataset [START_REF] Burri | The euroc micro aerial vehicle datasets[END_REF] and from 3 terrestrial sequences of a homemade dataset we acquired, which we call the SOT-TEVILLE dataset. 2) Multi-robot scenarios: The trajectories represented on figures (6a) and (6b) are estimated from the raw visualinertial measurements with VINS-MONO [START_REF] Qin | Vins-mono: A robust and versatile monocular visual-inertial state estimator[END_REF] and carried into the MAPLAB framework [START_REF] Schneider | maplab: An open framework for research in visual-inertial mapping and localization[END_REF] to be merged together (i.e. loop-closed and optimized). The resulting trajectory estimates are the best achievable when knowing the full set of visual-inertial measurements of all the involved robots. We consider them as groundtruth estimates and rely on them to simulate the inter-robot observations in the testing scenarios. We consider 4 multi-robot scenarios built using some subset of sequences taken from the EUROC and the SOTTEVILLE datasets, as summarized in table I. Scenario 1 confronts three disparate trajectories from the EUROC dataset, while the trajectories involved in scenarios 2 and 3 are similar and prone to inter-robot loop closures and observations. Scenario 4 involves the 3 trajectories from the SOTTEVILLE dataset. Note that while trajectory S3 is prone to loop closures and inter-robot observations with trajectories S1 and S2, this is not the case between those last two.

B. Method parameters

Table II gives the values of the user-defined parameters introduced throughout sections III and IV and used in the evaluation multi-robot scenarios.

C. Performance metrics

On each scenario, each tested method is evaluated regarding i) the accuracy of the estimation of its own trajectory and of the robots from which ones it has received data; ii) the communication load and iii) its matching performances regarding inter-robot loop closure detection. We classically use Absolute Root Mean Squared Errors (RMSE) for translation (ATE), rotation (ARE) and scale computed using [START_REF] Zhang | A tutorial on quantitative trajectory evaluation for visual(-inertial) odometry[END_REF]) to assess the estimation accuracy. For each robot, we evaluate

D. Evaluation of both methods

The tables IIIa and IIIb respectively give the a posteriori absolute RMSE for each robot estimating its own trajectory after having merged the information received from the other robots. The table IIIc shows for both methods, and for every robot in each scenario the a posteriori RMSE for each robot estimating the trajectories of the other robots from which it has received data. The table IIId displays the RMSE of the relative pose trajectories, estimated by the recipient robot, between sender and recipient robots. The tables IIIe and IIIf respectively detail the average amount of items each robot has communicated over the evaluated scenarios for both tested methods. Finally, the tables IIIg and IIIh quantify the number of inter-robot matchings (number of spotted loop closures or merged landmarks). In those tables, i ← j means that the related metric is associated to the trajectory of robot j estimated in the map of robot i.

First, the tables IIIa and IIIb show that fusing the information provided by the other robots can help the recipient robot to improve accuracy of the estimation of its own trajectory. When looking at tables IIIc and IIId, it comes that MKCS-VIP performs better than FS-VIP in terms of estimation accuracy of the trajectories of the other robots. Actually, MKCS-VIP brings multiple advantages over FS-VIP. First, each robot can relay the information it got from other robots. As a consequence, in scenario 4, as shown on the figure 6c, while robots associated to sequences S1 and S2 never observe each other and never observe the same scene, the robot S3 is able to relay the information of S1 to S2 and vice versa with MKCS-VIP. Secondly, it allows to exchange the loop closure factors each robot has spotted against its own map 3D structure. As a consequence, the resulting maps are more constrained than with FS-VIP for which one the marginalization process prevents any subsequent update on the exchange data. Additionally, it turns out that FS-VIP yields very conservative sparsified factors, such that many landmarks merges are needed (which is hardly reached as shown in table IIIg) so that the spotted inter-robot loop closures sufficiently constraint the trajectories of the robots.

VI. CONCLUSION

In this article, we introduced two decentralized data sharing methods based on the exchange of visual inertial packets built: from factor sparsification techniques (FS-VIP) and from landmark and raw measurement selection (MKCS-VIP). The performances of both methods have been investigated on two datasets. MKCS-VIP presents better results and takes advantages of exchanged and relayed data to build a complete map of the environment. Even if FS-VIP presents some limitations mainly due to the used consistent sparsification technique, it still has interesting results and could be useful for asymmetric computer power situation as this method is less heavy on the receiver side. As a perspective, we plan to conduct a more extensive testing on a wider dataset, possibly with conjoint aerial and terrestrial robots. This should bring more confidence on method parameters tuning and the analysis of their influence on the estimation accuracy and the communication load. Finally, both methods should be validated on real embedded systems.

Fig. 2 :

 2 Fig. 2: Overview of FS-VIP B. Packet precomputation

Fig. 3 :

 3 Fig. 3: Overview of the subgraph extraction, marginalization and sparsification process 1) Keyframe selection: When opening a new packet, the first step is to select online some keyframes along the trajectory, based on covisibility criteria which are looser than the ones used for regular keyframe selection. For each keyframe along the trajectory, we compute the number n c and the ratio r c of commonly observed landmarks with the lastly selected keyframe. If either n c ≤ n lim c or r c ≤ r lim c

Fig. 4 :

 4 Fig. 4: Pattern of the Jacobian matrix J S

Fig. 5 :

 5 Fig. 5: Overview of MKCS-VIP

Fig. 6 :

 6 Fig. 6: Groundtruth trajectories of the datasets and example of map merging using MKCS-VIP of the sender robot should be communicated at their first encounter. Each robot can consistently estimate the trajectory of the other robots, and such estimation can provide a base to register further data, such as local point-clouds asked on request.

  Table I gives some descriptive parameters of each robot trajectory, represented on figures (6a) and (6b).

	Seq.	Length Duration [m] [s]	ATE [m]	ARE [deg]	Scale	Scenarios 1 2 3 4
				EUROC dataset
	MH1	76	149	0.225 0.023 0.127
	MH2	69	183	0.240 0.020 0.122
	MH3	132	87	0.216 0.032 0.064
	MH4	92	87	0.364 0.023 0.122
	MH5	100	98	0.253 0.021 0.509
			SOTTEVILLE dataset
	S1	145	350	0.318 0.052 0.005
	S2	122	284	1.013 0.083 0.159
	S3	158	323	0.937 0.073 0.148

TABLE I :

 I Parameters of the sequences from the EUROC and SOTTEVILLE datasets and estimation accuracy of the trajectory estimation with VINS-MONO[START_REF] Qin | Vins-mono: A robust and versatile monocular visual-inertial state estimator[END_REF]

TABLE III :

 III Tables of results from the simulation on the EUROC and Sotteville datasets the absolute RMSE of its own trajectory estimation, the absolute RMSE of the trajectories of the other robots it reconstructed, and also the absolute RMSE of the trajectories

of inter-robot relative poses to assess the registration of each reconstructed map w.r.t. to the recipient robot's trajectory.

VII. ACKNOWLEDGEMENTS

We would like to thank the Société Nationale des Chemins de Fer (SNCF) for granting us the access to the warehouse where the above mentioned SOTTEVILLE dataset was acquired, as well as the Direction Générale de l'Armement (DGA) for financing this work.


Johanna Rüegg

Juan Carlos Quezada

Mathieu Santonja

Jaboury Ghazoul

Yakov Kuzyakov

Alexandre Buttler

Thomas Guillaume
email: thomas.guillaume@epfl.ch

Addresses

Drivers of soil carbon stabilization in oil palm plantations RUNNING TITLE: Soil carbon stabilization in oil palm plantations

Keywords:

Increasing soil organic carbon (SOC) in agroecosystems is necessary to mitigate climate change and soil degradation. Management practices designed to reach this goal call for a deeper understanding of the processes and drivers of soil carbon input stabilization. We identified main drivers of SOC stabilization in oil palm plantations using the well-defined spatial patterns of nutrients and litter application resulting from the usual management scheme. The stabilization of oil palm-derived SOC (OP-SOC) was quantified by δ 13 C from a shift of C4 (savanna) to C3 (oil palm) vegetations. Soil organic carbon stocks under frond piles were 20 and 22 % higher compared to harvest paths and interzones, respectively.

Fertilization and frond stacking did not influence the decomposition of savanna-derived SOC.

Depending on management zones, net OP-SOC stabilization equalled 16-27% of the fine root biomass accumulated for 9 years. This fraction was similar between frond piles and litter-free interzones, where mineral NPK fertilization is identical, indicating that carbon inputs from dead fronds did not stabilize in SOC. A path analysis confirmed that the OP-SOC distribution was largely explained by the distribution of oil palm fine roots, which itself depended on management practices. SOC mineralization was proportional to SOC content and was independent on phosphorus availability. We conclude that SOC stabilization was driven by C inputs from fine roots and was independent of alteration of SOC mineralization due to management. Practices favouring root growth of oil palms would increase carbon sequestration in soils without necessarily relying on the limited supply of organic residues.

INTRODUCTION

Soil organic carbon (SOC) depletion in agroecosystems is a major source of greenhouse gas emissions, resulting in losses of soil fertility and ecosystem stability [START_REF] Amundson | Soil and human security in the 21st century[END_REF].

Policymakers have recently acknowledged the promotion of soil C sequestration in agroecosystems as a promising strategy to simultaneously mitigate climate change and enhance food security [START_REF] Lal | Beyond COP 21: Potential and challenges of the "4 per Thousand" initiative[END_REF]. Nonetheless, management practices favouring soil C sequestration lead to highly variable outcomes, calling for a deeper understanding of processes and factors controlling C stabilization, especially after land-use changes [START_REF] Ghimire | Tillage, crop residue, and nutrient management effects on soil organic carbon in rice-based cropping systems: A review[END_REF][START_REF] Haddaway | How does tillage intensity affect soil organic carbon? A systematic review[END_REF]. This is particularly needed in regions that are undergoing rapid and substantial land degradation following conversion from largely forested landscapes to intensive agricultural systems, as has been the case in tropical landscapes dominated by oil palm (Elaeis guineensis) plantations. In Sumatra, a region with a long history of oil palm cultivation, a significant proportion of oil palm plantations has already reached a critical low level of SOC content (< 1%) in the topsoil [START_REF] Guillaume | Soil degradation in oil palm and rubber plantations under land resource scarcity[END_REF].

Stocks of SOC depend on the balance between soil C inputs from vegetation and outputs from SOC mineralization, erosion and leaching [START_REF] Lorenz | Soil Carbon Stock[END_REF]. This balance is strongly affected by the conversion of natural ecosystems to intensive agricultural land [START_REF] Guillaume | Carbon Costs and Benefits of Indonesian Rainforest Conversion to Plantations[END_REF]. Predicting impacts of land-use change and management on SOC dynamics and its stabilization faces major difficulties as many factors affect both litter input and mineralization processes that determine the fraction of C input stabilized in SOC. While higher plant biomass inputs might lead to higher SOC accumulation, an increase of fresh organic matter inputs may enhance the mineralization of more recalcitrant SOC that ultimately reduces the gain in SOC, a process known as priming [START_REF] Kuzyakov | Review of mechanisms and quantification of priming effects[END_REF]. Nutrient application modifies the stoichiometry of organic matter inputs and of soil organic matter, which in turn affects microbial processes controlling SOC stabilization [START_REF] Qiao | Carbon and nitrogen additions induce distinct priming effects along an organic-matter decay continuum[END_REF][START_REF] Zang | N fertilization decreases soil organic matter decomposition in the rhizosphere[END_REF]. For instance, altered microbial carbon use efficiency (i.e. the ratio of C incorporated into microbial biomass to the added C) or mining for nutrients from recalcitrant SOC result in either SOC gains or losses [START_REF] Finn | Ecological stoichiometry controls the transformation and retention of plant-derived organic matter to humus in response to nitrogen fertilisation[END_REF][START_REF] Kirkby | Nutrient availability limits carbon sequestration in arable soils[END_REF]. Quantifying soil C inputs remains a methodological challenge and few data are available, for instance, on root turnover, rhizodeposition, and the fraction of aboveground litter C stabilized in SOC [START_REF] Pausch | Carbon input by roots into the soil: Quantification of rhizodeposition from root to ecosystem scale[END_REF]. Experimental data are especially limited for perennial plants in tropics and subtropics.

The fast SOC turnover in the tropics makes soils particularly sensitive to land-use change [START_REF] Guillaume | Losses of soil carbon by converting tropical forest to plantations: Erosion and decomposition estimated by δ 13 C[END_REF][START_REF] Pabst | Land Use and Precipitation Affect Organic and Microbial Carbon Stocks and the Specific Metabolic Quotient in Soils of Eleven Ecosystems of Mt. Kilimanjaro, Tanzania[END_REF][START_REF] Zech | Factors controlling humification and mineralization of soil organic matter in the tropics[END_REF]. For example, soil C inputs decrease up to 90% when rainforests are converted to oil palm plantations, resulting in a rapid drop of SOC [START_REF] Guillaume | Carbon Costs and Benefits of Indonesian Rainforest Conversion to Plantations[END_REF].

Soil organic C losses are not, however, uniform within plantations, and specific management zones within the plantation may even exhibit a gain in SOC [START_REF] Khasanah | Carbon neutral? No change in mineral soil carbon stock under oil palm plantations derived from forest or non-forest in Indonesia[END_REF][START_REF] Rahman | Changes in soil organic carbon stocks after conversion from forest to oil palm plantations in Malaysian Borneo[END_REF]. Soil C inputs in mature oil palm plantations without cover crops (the most common practice) occurs mostly belowground through rhizodeposition because the understorey is frequently cleared. Significant aboveground C inputs occur only under frond piles, i.e. zones where dead fronds are stacked. A gradient of fertilizer application is superimposed on the gradient of soil C input: most fertilizers are applied around the trunks and, in some cases, additionally to the whole surface area at lower rate, including the frond piles but excluding the harvest path. These management practises lead to characteristic management zones with specific factors affecting SOC dynamics (Fig 1a). With plantation ageing, the heterogeneity of SOC distribution increases, depending on the distance to the tree, the presence of frond piles and fertilizer applications [START_REF] Frazão | Soil carbon stocks and changes after oil palm introduction in the Brazilian Amazon[END_REF][START_REF] Goodrick | Tree-scale spatial variability of soil carbon cycling in a mature oil palm plantation[END_REF]. Carbon and nutrients cycling are, thus, highly heterogeneous leading to a variable SOC equilibrium depending on management practices. This highlights the possibility to promote the increase of overall SOC stocks by redesigning management practices to reach this goal.

In Colombia, oil palm plantations are often established on native savanna grasslands dominated by grasses with a C4 photosynthesis pathway. The shift from C4 (grasses) to C3 vegetation (oil palm) allows source determination in soil organic matter using its 𝛿 13 C signature [START_REF] Balesdent | Natural 13C abundance as a tracer for studies of soil organic matter dynamics[END_REF]. The aim of the present study is to disentangle the effects of soil C and nutrient inputs on newly accumulated SOC (oil palm-derived) and the decomposition of old SOC (savanna-derived), taking advantage of the specific patterns of fertilization and soil C inputs in the four management zones. We hypothesized that soil C inputs and consequent SOC accumulation increase with fine root density, which itself depends on plantation age and distance to trees, as well as on the presence of frond piles. Fertilizer application is, however, expected to decrease SOC stabilization rates where soil C inputs are low. Hence, the specific aims of the study are to i) quantify the new oil palm-derived (C3 signature) and the old savanna-derived SOC stocks (C4) in a mature oil palm plantation established on native savanna grassland, ii) assess oil palm root development and its impact on SOC accumulation with plantation age, iii) determine the impacts of management practices on root development and soil microbial activity and iv) identify the main factors (root density, nutrient availability, microbial activity) driving SOC stabilization.

MATERIALS AND METHODS

Study area

The study was conducted in the Eastern Plains (Llanos Orientales), Department Meta, Colombia (4°05'7.0"N, 71°53'59.0"W). The region experiences a tropical climate (mean annual temperature of 26 °C and yearly precipitation of 2200 mm yr -1 ) with a distinct dry season from December to March and 95% of the yearly rain falling between April and November [START_REF] Lavelle | Soil ecosystem services and land use in the rapidly changing orinoco river basin of colombia[END_REF][START_REF] Rippstein | Condiciones naturales de la sabana[END_REF]. The study site lies in the slightly undulating well-drained high plains (Altillanura plana) dominated by Plinthosols and Ferralsols (IUSS Working Group WRB, 2014). These soils have a low fertility, high acidity and high aluminium saturation limiting agricultural production [START_REF] Lavelle | Soil ecosystem services and land use in the rapidly changing orinoco river basin of colombia[END_REF][START_REF] Rippstein | Condiciones naturales de la sabana[END_REF]. The natural vegetation is an herbaceous savanna with scarce bushes, which is drained by many small rivers. Gallery forests (morichales) grow in the depressions along these rivers [START_REF] Rippstein | Condiciones naturales de la sabana[END_REF].

Two unmanaged native savanna sites and three oil palm plantations with increasing age (2-, 4-, and 9-year old) were selected within an area of approximately 8 × 8 km and sampled in July and August 2016 (Fig. S1). Sites were carefully chosen with the help of plantation agronomists to ensure that the investigated plantations had been established on unmanaged native savannas, i.e. no cattle grazing in the past. Soils in the five selected sites were classified as sandy-loam Ferralsols (clay content and pH ranging from 11 to 16% and from 4.5 to 5.3 in the top 30 cm, respectively) with compacted top soil (bulk density ranging from 1.31 to 1.46 g cm -3 in the top 10 cm). Drivers of SOC stabilization were assessed in the 9-year old plantation. The 2-and 4-year old plantations were selected to assess the development of oil palm rooting system and its relationship with SOC accumulation.

Plantation management

Oil palm plantations had been established on native savannas whose soils were first loosened with a chisel plough (to a depth of 40 cm) and tilled with an overturning plough (10 cm).

Liming and phosphate rocks were applied before planting palms. Oil palm trees had been planted in a triangular grid pattern with a distance of 9 m between trees, leading to 143 palms ha -1 (Fig. 1b). In the young plantations (2-and 4-year old), a mixture of Kudzu (Pueraria phaseoloides) and Desmodium (Desmodium heterocarpon subsp. ovalifolium) cover crops (C3 vegetation) have been implemented after planting. Cover crops were never used in the mature plantation (9 years), and in this case the soil was always kept bare after planting with oil palms. Accordingly, C3-derived SOC in the 9-year old plantation are attributed solely to oil palm-derived SOC.

Oil palm management leads to four well-defined management zones in productive plantations (about 4 years after establishment). In the study region, the weeded circle (WC) around the oil palm trunk is always kept free of vegetation (Fig. 1). At young ages (2-and 4year old), all fertilizers are applied in WC. Associated with the beginning of harvest after 4 years, pruning starts and fronds are piled up in between palm trees (frond pile, FP). Each second avenue between palm lines becomes a harvest path (HP), where machines circulate.

In productive plantations (after 4 years), fertilizers are evenly spread from the harvest path by machines, i.e. all management zone receive the same amount of fertilizers except the harvest path that receive none. The remaining area, especially the alternating avenue, represents the fourth zone, where the soil is kept bare and fertilization starts after 4 years (interzone, IZ).

The relative surface area of the four management zones were 60% (IZ), 18% (HP), 12% (WC) and 10% (FP). Fertilization depends on oil palm stand age. The 9-year old plantation received per hectare during the first year about 240 kg of NPK, 60 kg of kieserite (MgSO4), and 25 kg of KCl, of zinc and of boron. The amount of NPK increased over time while the use of other types of fertilizer varied from year to year. The year before sampling, 600 kg of NPK, 200 kg of phosphate rocks and 160 kg of KCl-MgO per hectare were spread in the plantation.

Soil and roots sampling

Plots of 1 ha were established in all plantations and savannas in areas with homogeneous soil, far from the influence of roads, rivers or groundwater, free of laterite formations (arecife) and without former amendment of compost or residues from processing oil mills. Five trees as replicates for management zones were selected randomly in each plantation. For each of the five replicate of trees, the adjacent management zones (FP, HP, IZ, and WC) were sampled on a systematic grid (Fig. 1). Frond piles, harvest paths and interzones were sampled at the same distance from trees (4.5 m) to assess the effects of management starting 4 years after establishment. Two additional points in the interzone were sampled at 2 and 3 m away from the trunk to assess the horizontal expansion of oil palm roots (IZ2, IZ3, Fig. 1). In the 2-and 4-year old plantations, only the weeded circles and interzones were sampled because of the absence of frond piles and harvest path in young plantations. In the 2-year old plantation, the point IZ3 was not sampled, as root densities were already low in IZ2 (Fig. S2). In the savanna plots, five sampling points were selected at regular distances along a 100 m transect.

Roots and soils were sampled with a cylindrical corer of 5 cm diameter at three depth intervals (0-10, 10-20 and 20-30 cm). Roots were separated from soil by sieving at 2 mm and rinsed to remove attached mineral particles. Fine roots that passed through the sieve were manually picked. For the 2-and 4-year old plantations outside of weeded circles, cover crop roots were removed from oil palm roots. Oil palm roots were divided into coarse roots (> 2 mm, corresponding to primary and secondary roots of oil palms) and fine roots (< 2 mm, corresponding to tertiary and quaternary roots of oil palms). Dry root biomass was determined after drying at 60 °C for 48 h. Soil samples were air-dried and sieved at 2 mm directly after collection and further oven dried at 40 °C for 48 h prior to laboratory analyses.

Soil analysis

Total C and nitrogen (N) contents in soil, as well as δ 13 C signature were determined at the University of Göttingen with an isotope ratio mass spectrometer (Delta Plus, Finnigan MAT, Bremen, Germany). Because of the absence of carbonates in acidic soils, total C represents organic C. Residual water content was assessed by drying soil samples at 105 °C for 24 h. Bulk density was measured by inserting horizontally two cylinders of 100 cm 3 per depth at 0-5, 5-10, 10-20 and 20-30 cm depth in a soil pit located at the centre of each sampling plot.

Bulk density in each cylinder was determined after drying at 105 °C and averaged between four cylinders for 0-10 cm depth interval and 2 cylinders for 10-20 and 20-30 cm depth intervals. Carbon stocks were calculated multiplying C contents with bulk density and the layer thickness and the respective fractions of C4 (savanna-derived C) and C3 (oil palmderived) SOC. Soil available phosphorous (P) was determined using Bray II extraction method [START_REF] Bray | Determination of total, organic, and available forms of phosphorous in soils[END_REF]. Three grams of soil were extracted by shaking for 15 min with 20 ml of 0.03 N NH4F and 0.025 N HCl. Filtered extracts were mixed with a colorimetric reagent ((NH4)6Mo7O24 -SnCl2) and absorbance was measured with a UV/VIS spectrometer at 660 nm (Lambda 35, Perkin Elmer, Buckinghamshire, United Kingdom).

C3 and C4 derived carbon

Relative portion of oil palm-and savanna-derived C in SOC were calculated based on the differences in δ 13 C signature of biomass between savanna grassland, dominated by C4 photosynthetic pathway, and oil palms, C3 photosynthetic pathway, using two sources linear isotopic mixing model [START_REF] Balesdent | Natural 13C abundance as a tracer for studies of soil organic matter dynamics[END_REF]:

𝑓 𝑂𝑃_𝑑 = 𝛿 13 𝐶 𝑆_𝑑 -𝛿 13 𝐶 𝑁𝑆_𝑑 𝛿 13 𝐶 𝑂𝑃_𝑏 -𝛿 13 𝐶 𝑁𝑆_𝑏
where fOP_d is the fraction of oil palm (C3)-derived SOC, δ 13 CS_d is the isotopic signature measured in a soil sample at depth d, δ 13 CNS_d is the average isotopic signature at the corresponding depth in natural savanna reference sites, δ 13 COP_b is the averaged signature of oil palm fine roots in the mature oil palm plantation (mean = -28.2‰, standard deviation (SD) = 0.18, n = 9) and δ 13 CNS_b is the average isotopic signature of the aboveground and belowground savanna biomass determined on representative subsamples of aboveground biomass collected on 1 m 2 and fine roots collected with soil cores at each sampling point (mean = -13. 6 ‰, SD = 1.1, n = 19). This approach assumes that the 13 C fractionation occurring during the integration of biomass into SOC at each soil depth is the same for savanna and oil palm biomass [START_REF] Pausch | Soil organic carbon decomposition from recently added and older sources estimated by 13C values of CO2 and organic matter[END_REF]. No 13 C fractionation was observed between savanna plant biomass and savanna SOC in 0-10 cm layer (-13.7 ‰, SD = 0.5, n = 10).

Carbon stabilization per cumulative standing root biomass

Fine root stocks served as a proxy for the C input from oil palm roots. To account for differences in root stocks between management zones and root development time depending on the distance to the palm tree, cumulative standing fine root biomass in each management zone was estimated by fitting a linear model on root biomass measured in various plantation ages. The model was integrated starting from plantation establishment for WC, but only starting from 4 years for IZ, FP and HP, i.e. when oil palm rooting system reached 4.5 m away from the tree (Fig. S2). Assuming constant belowground C input per unit of fine roots [START_REF] Pausch | Estimation of rhizodeposition at field scale: Upscaling of a14C labeling study[END_REF], oil palm-derived SOC stocks were normalized per unit of cumulated fine root biomass to assess the effect of management zones on net C3 stabilization efficiency.

Soil incubation

Microbial biomass was analysed after incubation of topsoils (0-10 cm) from the 9-year old plantation (zones FP, HP, IZ and WC). Twenty grams of dry soil were rewetted to 60% WHC and incubated at 25 °C for 31 days. Jars were ventilated and weighed every week and rewetted once to compensate for the evaporated water. Carbon and nitrogen (N) in microbial biomass were measured at the end of the incubation by the fumigation-extraction method [START_REF] Vance | An extraction method for measuring soil microbial biomass C[END_REF]. For C and N contents, 5 g of incubated soil were fumigated for 24 h with ethanolfree CHCl3 in a desiccator. Soils were extracted by shaking 1 h in 25 ml solution of 0.5 M K2SO4 and then filtered. Non-fumigated samples were processed in parallel. Total extractable organic C and N were analysed with a TOC-N analyser (Shimadzu, Kyoto, Japan). Extractable C in the non-fumigated samples was assumed to represent dissolved organic carbon (DOC).

Microbial biomass C and N were calculated as the difference between fumigated and nonfumigated samples, which were used also to calculate microbial C:N ratio. Microbial biomass C was corrected by dividing extractable C with a factor of 0.45 [START_REF] Beck | An inter-laboratory comparison of ten different ways of measuring soil microbial biomass C[END_REF], N with a factor of 0.54 [START_REF] Brookes | Chloroform fumigation and the release of soil-nitrogen -a rapid direct extraction method to measure microbial biomass nitrogen in soil[END_REF].

Basal respiration

Basal respiration of the rewetted samples was measured with the MicroResp TM kit [START_REF] Campbell | A rapid microtiter plate method to measure carbon dioxide evolved from carbon substrate amendments so as to determine the physiological profiles of soil microbial communities by using whole soil[END_REF]. Three analytical replicates of 0.5 g for each field replicate were taken from the incubation jars after rewetting and incubated in 96 deepwell plate in parallel to the incubation in jars. Soil was kept moist by a moist paper towel fixed on the plate. Respiration was measured 1, 2, 4, 8, 15, 18, 24 and 31 days after rewetting using a MicroResp TM kit [START_REF] Campbell | A rapid microtiter plate method to measure carbon dioxide evolved from carbon substrate amendments so as to determine the physiological profiles of soil microbial communities by using whole soil[END_REF]. To remove residual CO2 in the wells, the plate was aerated with a fan before incubating for 6 h with the indicator plate on top (at 27 °C). The indicator plate was read before and after incubation with a spectrophotometer (Microplate reader BioTek SynergyMX) at 570 nm. Absorption calibration was done by dissolving a known amount of NaCO3 with 1 M HCl in excess in closed jars with eight microwells of the indicator plate for 6 h [START_REF] Campbell | A rapid microtiter plate method to measure carbon dioxide evolved from carbon substrate amendments so as to determine the physiological profiles of soil microbial communities by using whole soil[END_REF].

Soil respiration was partitioned using a two-pool mixed-model to describe SOC mineralization kinetics [START_REF] Bonde | Nitrogen Mineralization Kinetics in Soil During Longterm Aerobic Laboratory Incubations: A Case Study[END_REF]. The first pool follows a first-order decomposition kinetics, while the second follows a zero-order kinetics, corresponding to the stabilized basal respiration:

𝐶𝑚𝑖𝑛 = 𝐶 𝑙 (1 -𝑒 -𝑡𝑘 𝑙 ) + 𝐵𝑅 𝑡
where t is the time, Cmin is the cumulative CO2 mineralized to time t, Cl is the labile C pool released from sample preparation and re-wetting, kl is the decomposition constant of the labile pool and BR is the basal respiration. The basal respiration of each field replicate was determined by fitting the model on the three analytical replicates. The metabolic quotient is the ratio of basal respiration over microbial biomass C (Cmic).

Statistical analyses

All statistical analyses were performed using the open source software R version 3.2.1 (R Core [START_REF] Team | R: A language and environment for statistical computing[END_REF]. One sampling point in WC and one in IZ were removed from all analyses due to very low δ 13 C signatures, probably resulting from the former presence of C3 bush at these exact sampling locations. We used a linear mixed-effects model approach ("lme4" stocks vs. fine roots) were assessed by linear regressions. Average C stocks at plot scale down to 30 cm depth were calculated using the relative surface area of each management zone.

Associations among parameters were calculated using Pearson correlation. P-values were determined using the function cor.test.

A path analysis was performed to disentangle the direct and indirect effects of fine roots, soil microorganisms and nutrient application on the accumulation of oil palm-derived SOC by using the "lavaan" package [START_REF] Rosseel | lavaan: An R package for structural equation modelling[END_REF]. Based on the priori knowledge, we developed an initial conceptual model that was both consistent with our data and which made biological sense (Fig. S3). We first created five conceptual groups of measured variables, which represented i) nutrient application (measure of available P), ii) fine roots C inputs (measure of fine roots biomass), iii) SOM quality (measure of C:N ratio of SOM), iv) SOC accumulation (measure of oil palm-derived SOC) and v) soil microorganism effects. As a proxy for the microbial effect, we used the score of each sample on the first axis (PC1 = 68 %) of a principal component analysis (PCA) including all microbial related variables (soil basal respiration, metabolic quotient, Nmic, Cmic, microbial C:N ratio) (Fig. S4). The conceptual model hypothesized that fine roots biomass and microbial activity have a direct impact on oil-palm SOC accumulation. Both might be affected by nutrient application.

Additionally, microbial activity would be directly affected by root density through the amount of rhizodeposition and indirectly by changing organic matter quality. Finally, P availability could have a direct effect on microorganisms and SOM quality by changing the resources' stoichiometry in the soil. Pedoclimatic factors affecting SOC stabilization (e.g. soil mineralogy) were not specified in the model because they are identical between zones and cannot be modified by management. The adequacy of the model was determined by nonsignificant differences between the predicted and observed covariance matrices (χ2 tests, p > 0.05), low root mean squared error of approximation index (RMSEA < 0.1), high Tucker-Lewis index (TLI > 0.90) and high comparative fit index (CFI > 0.90) [START_REF] Grace | Structural equation modeling and natural systems[END_REF][START_REF] Rosseel | lavaan: An R package for structural equation modelling[END_REF].

RESULTS

Soil organic carbon stocks and origin

Observations done on the 9-year old plantation showed lower soil C stocks compared to the native savanna grasslands (Fig. 2). On average, 1.0 ± 0.2 kg C m -2 was lost down to 30 cm depth, considering the relative area of each management zone. However, spatial SOC distribution in the mature oil palm plantation depended on management zones. Soil C losses during that period of time in the top 10 cm under frond piles and, to a lesser extent, under weeded circles were lower than C losses in harvest paths and interzones (Fig. 2). This trend was similar down to 30 cm depth but the differences were not significant below 10 cm depth.

Differences in SOC stocks between management zones arose mainly from a higher accumulation of oil palm-derived SOC (C3-derived) under frond piles and weeded circles (Fig. 2). Oil palm-derived SOC after 9 years already accounted for between 27% (IZ) and 45% (FP and WC) of the total SOC stock in the top 10 cm, where differences between zones were highest. Below 20 cm depth, the contribution of oil palm-derived SOC dropped, accounting for 2% (IZ) and maximum 12% (WC) of the total SOC stocks. The amount of savanna-derived SOC (C4-derived) remaining after 9 years was lower in the top 10 cm than between 10-20 cm depth, indicating a faster decomposition of this C pool in the top soil compared to deeper soil layers. Nonetheless, management zones had little influence on the decomposition rates of savanna-derived SOC, except in 10-20 cm depth under weeded circles where more C was lost than under frond piles. A similar trend was observed in the top 10 cm that explains why total SOC under weeded circles was intermediate as compared to frond piles and harvest paths, despite the high amount of oil palm-derived SOC stabilized under weeded circles.

Root development and C inputs

Oil palm fine roots were first observed at 4.5 m away from palm trees in two out of the five investigated palms in the 4-year old plantations (Fig. S2). This indicates that oil palm rooting systems from adjacent palms started to overlap at that age, but root biomass at 4.5 m was still very low. Fine root biomass under weeded circles increased constantly during 9 years (Fig. 3a). At the age of 9 years, roots were observed in all management zones but root growth was strongly enhanced under frond piles, reaching the same fine root biomass as under weeded circles in only 5 years, well above fine root biomass in the interzones and harvest paths (Fig. 3a).

Oil palm-derived SOC under weeded circles was highly correlated with fine root biomass (Fig. 3b). Oil palm-derived SOC stocks corresponded to 70  5 % of fine root biomass stocks (slope = 0.70, R 2 = 0.80) at the time of measurement, independently of soil depth and plantation age. This percentage was relatively constant between the three soil depths; from 62  10% in the top 10 cm to 50  11% between 20-30 cm depth (Table 1). Oil palm-derived SOC accumulation per year was 3 to 4 times faster in the top soil 10 cm than in the underlying layers (Table 1).

While oil palm-derived SOC in the weeded circle was accumulating from the beginning of the plantation, root biomass was measured at fixed time points. Oil palm fine root biomass was integrated over the whole duration of the plantation to calculate the cumulative fine roots biomass for each year and soil depth under weeded circles. Accordingly, oil palm-derived SOC accumulation corresponded to 14  2%, 11 2% and 11  1% of the cumulated fine roots biomass stocks under weeded circles at 0-10, 10-20 and 20-30 cm depth, respectively (Table 1). Since soil depth had little impact on the relationship between oil palm-derived SOC and fine root biomass, differences in oil palm-derived SOC stabilization rates between depths were mainly related to differences in fine root biomass. The linear relationship across age showed that the proportion of oil-palm derived SOC stabilizing per amount of fine roots present in the plantation remains constant at least during 9 years (Fig. 3b).

Oil palm derived-SOC and fine root biomass were similar under frond piles and weeded circles in the top 10 cm in the 9-year old plantation (Fig. 3a). However, fine roots appeared at 4.5 m away from palms 4 years later than under weeded circles. To remove the effects of root density and duration of C inputs between zones, oil palm-derived SOC stocks were also divided by the cumulative roots biomass stocks and compared between management zones (Fig. 4).

After this normalization to the amount of fine root biomass, the stabilization of oil palmderived SOC was similar between frond piles and interzones, two zones receiving the same amount of mineral fertilizers and sampled at the same distance to the tree. Consequently, only little C from the large amount of C present in dead fronds is eventually stabilized in SOC. The stabilization of oil palm-derived SOC was lower under weeded circles, the zone receiving the largest amount of mineral fertilizer, than under harvest paths, the only zone experiencing no direct application of mineral fertilizers.

Soil microorganisms and fertility

Microbial activity (basal respiration and metabolic quotient) and biomass (microbial biomass C and microbial C:N ratio) parameters were similar under frond piles and weeded circles (Fig. 5). In these two zones, the same amount of microbial biomass respired more C compared to harvest paths and interzones. Consequently, the metabolic quotient was highest under frond piles and weeded circles.

Available P (Bray II) was an order of magnitude higher in weeded circles compared to the other zones, as expected by the higher fertilization application in that zone. Weeded circles also had the lowest amount of K2SO4-extractable C (DOC) despite high SOC content. (Table S1). Between the three zones (HP, IZ and FP) located at the same distance to the palm trees but varying in their management, P availability was two times lower in harvest paths (11.1 ± 0.3 µg g -1 ) as compared to frond piles, but well above the P availability in native savanna sites

(2.2 ± 0.3 µg g -1 ), despite the absence of direct fertilization. The DOC amount was similar between the three zones and the C:N ratio was 1 unit higher under frond piles, indicating that only small change in SOC quality occurred between management zones.

Drivers of the accumulation of oil palm-derived SOC

In the three zones located at the same distance to the palm trees (HP, IZ and FP), oil palmderived SOC stocks (C3) were highly correlated with fine root biomass (r = 0.82), as well as to most microbial parameters, the amount of soil organic matter (C and N contents) and available P (Fig. 6). The basal respiration, the metabolic quotient and the microbial biomass N increased with higher oil palm-derived SOC but not with total C content, underlying the role of fresh organic C to maintain microbial activity. By contrast, higher C:N ratio in microbial biomass and, to a lesser extent, higher microbial biomass were associated to less oil palm-derived SOC and less fine roots biomass but with more savanna-derived SOC remaining after conversion to oil palm.

A path analysis was performed to disentangle the direct and indirect effects of fine roots, soil microorganisms and nutrient availability on the accumulation of oil palm-derived SOC at 4.5 m away from palm trees in the 9-year old plantation. Because of the high association among all microbial parameters (Fig. 6), the scores of each sample on the first axis (PC1 = 68 %) of a principal component analysis (PCA) including all microbial related parameters were used as a proxy for soil microorganisms in the path analysis (Fig. S4). The fitting parameters of the model were good (Pχ 2 = 0.33, RMSEA = 0.09, TLI = 0.97, CFI = 0.99), and the model explained 71% of the variance in oil palm-derived SOC (Fig. 7). Fine root biomass was strongly influenced by nutrient availability (available P) and was an important driver of microbial properties (i.e. increase in microbial biomass N, basal respiration, specific respiration and metabolic quotient and, in the opposite, decrease of microbial biomass C and C:N ratiosee also Fig. 6). By contrast, nutrient availability had only a marginal (p = 0.09) and opposite effect on soil microorganisms. While fine root biomass influenced SOM quality (C:N ratio), microorganisms were not affected by SOM quality. The direct effect of fine roots on oil palm-derived SOC stabilization was 2.8 times stronger than the marginal effect (p = 0.06) of soil microorganisms. Accordingly, nutrient availability had an important indirect impact on the accumulation of oil palm-derived SOC by favouring root development and thereby C inputs without enhancing SOC mineralization.

DISCUSSION

Drivers of SOC stabilization

Soil organic C stocks strongly varied depending on the management zones of the plantation.

Specific management impacted the amount of new oil palm C input and its stabilization into SOC but had little effect on the decomposition rate of old savanna-derived SOC (Fig. 2).

Carbon stabilization was mainly driven by C inputs from fine roots rather than by changes in C outputs from microbial mineralization (Fig. 7). The relationship between fine roots and oil palm derived SOC remained weakly affected by soil depth, plantation age and management zones (Table 1, Fig. 3 and4). Fine root biomass was a good proxy for soil belowground C inputs. Oil palm fine roots absorb nutrient and water, and therefore have definite growth and short-term self-pruning. In contrast, coarse roots, which have the function of conduction, have indefinite growth and long term self-pruning [START_REF] Jourdan | Root system architecture and gravitropism in the oil palm[END_REF][START_REF] Jourdan | Architecture and development of the oil-palm (Elaeis guineensis Jacq.) root system[END_REF]. The relationship between cumulated fine root biomass and oil-palm derived SOC stocks would not remain constant once stocks reach equilibrium. The relationship, however, was linear (Fig. 3), indicating that SOC stocks were still far from equilibrium after 9 years of cultivation and that oil palm-derived SOC will continue to increase.

Microbial metabolism and biomass were strongly affected by fine roots. Higher organic C availability in terms of oil palm-derived SOC and C inputs under frond piles and in the weeded circle were associated with microbial communities characterized by a high mineralization activity but of low efficiencythe so termed r strategy [START_REF] Loeppmann | Substrate quality affects kinetics and catalytic efficiency of exo-enzymes in rhizosphere and detritusphere[END_REF]. The high metabolic quotient results from either a low C use efficiency or a high microbial biomass turnover. By contrast, zones with low C availability such as the harvest path and the interzones were associated with microbial communities more efficient to maintain their biomass despite lower mineralization rates.

Their lower metabolic quotients and their higher C:N ratios suggest that the scarcity of C increased the proportion of K strategists and fungi within microbial communities [START_REF] Mouginot | Elemental stoichiometry of fungi and bacteria strains from grassland leaf litter[END_REF][START_REF] Six | Bacterial and Fungal Contributions to Carbon Sequestration in Agroecosystems[END_REF]. While the main effect of nutrient availability on microbial communities was indirect by increasing roots C inputs, it tended to have also a minor direct but contrasting effect on microbial communities. Lower C to nutrients ratio favours high C use efficiency of microorganisms, which would explain the larger microbial biomass in interzones despite the low microbial activity and root density [START_REF] Sinsabaugh | Carbon use efficiency of microbial communities: Stoichiometry, methodology and modelling[END_REF]. Nonetheless, management impacts on microbial communities and their resource consumption strategies levelled-off resulting in similar fraction of SOC mineralized in all management zones. Weeded circles were the only area of the plantation showing slightly lower net SOC stabilization and savanna-derived SOC stocks (Fig. 2 and4). The very high amount of fertilizer applied already in the early stage of the plantation might have slightly fastened SOC turnover. Nonetheless, management effects were small in regards to the 10-fold increase of P availability in this zone as compared to the rest of the plantation. The fact that nutrient availability had little impact on the SOC mineralization and stabilization indicates that microorganisms were mostly C-limited. Indeed, soils were depleted in SOC and the whole surface area of the plantations, even areas not directly fertilized, experienced an increase in P availability as compared to native savannas.

In summary, management and its impact on soil belowground C inputs and nutrient availability did not lead to priming of recalcitrant SOC nor to faster turnover of fresh organic SOC, except around palms in the weeded circles. Hence, SOC distribution was not driven by an altered decomposition of SOC pools but by different rates of fine root growth and the resulting soil C inputs depending on management.

Drivers of soil C inputs

Oil palm rooting system follows a relatively rigid and genetically determined development [START_REF] Jourdan | Root system architecture and gravitropism in the oil palm[END_REF][START_REF] Jourdan | Architecture and development of the oil-palm (Elaeis guineensis Jacq.) root system[END_REF]. Primary lateral roots start to grow one year after germination at a rate of 3 mm d -1 , confirming that roots of palms, which are planted a year after germination, reach 4.5 m in the fourth year of a plantation [START_REF] Jourdan | Architecture and development of the oil-palm (Elaeis guineensis Jacq.) root system[END_REF].

Nevertheless, root development showed plasticity, reacting to management as shown by the fast development of fine roots under frond piles. Root development is generally enhanced in nutrients-rich zones [START_REF] Hodge | The plastic plant: Root responses to heterogeneous supplies of nutrients[END_REF]. Phosphorous availability was an important driver of root development, favoring soil C inputs from fine roots, and consequently the accumulation of new SOC (Fig. 7). Its distribution did not exactly reflect the pattern of fertilizer application (Table S1). Despite no direct fertilization applications, harvest paths were enriched in available P as compared to the reference savanna sites and only frond piles exhibited a significantly higher P availability. This discrepancy may result from the initial soil preparation with dolomite and phosphate rocks and seems to be maintained over time by the recycling of organic P inputs from rhizodeposition, root turnover, and frond mineralization.

The positive impact of frond piles on root development is likely not limited to P availability.

Fronds piles increase the availability of major cations [START_REF] Law | Spatial variability of soil organic carbon in oil palm: a comparison between youg and mature stands[END_REF]. Similarly, it was shown that the application of empty fruit bunches on the surface increases soil moisture [START_REF] Tao | Effects of soil management practices on soil fauna feeding activity in an Indonesian oil palm plantation[END_REF] and roots development [START_REF] Kheong | Empty fruit bunch application and oil palm root proliferation[END_REF]. The positive impact of frond piles on SOC stocks has been previously reported [START_REF] Haron | Microbial biomass and soil organic matter dynamics in oil palm (Elaeis guineensis jacq.) plantations, West Malaysia[END_REF][START_REF] Law | Spatial variability of soil organic carbon in oil palm: a comparison between youg and mature stands[END_REF] but authors have already highlighted the small increase of SOC stocks given the huge amount of frond's biomass C concentrated on a small surface area (2-3 kg C m -2 yr - 1 in frond piles). It was suggested that fronds were mineralized mostly aboveground with little contribution to SOC [START_REF] Haron | Microbial biomass and soil organic matter dynamics in oil palm (Elaeis guineensis jacq.) plantations, West Malaysia[END_REF]. Our findings confirm this hypothesis and indicate that the positive impacts of frond piles arise more from the improvement of soil conditions (likely nutrient availability, humidity, protection from erosion) that favours root growth than from their role as a C source.

Increasing SOC stocks

Soil organic C stabilization depends on belowground C inputs and not on the management induced variation of SOC mineralization rates (Fig. 7). This field evidence-based finding is similar to the conclusion of a modelling study on arable cropping system in temperate zones [START_REF] Autret | Alternative arable cropping systems : A key to increase soil organic carbon storage ? Results from a 16 year fi eld experiment[END_REF]. Soil C inputs should be enhanced to increase soil C sequestration in oil palm plantations. The current management practice that consist in piling dead fronds, however, is not efficient to integrate the fronds' organic matter into the soil. This can be explained by several mechanisms. First, organic matter quality in fronds is low due to high C to nutrient ratios [START_REF] Yusuyin | Nutrient Status of Frond Heaps and the Underlying Soils at An 18-Year-Old Oil Palm Field in Central Pahang, Malaysia[END_REF]. This decreases the C use efficiency of microorganisms. Second, bioturbation from soil fauna is limited as indicated by a sharp transition between decomposing fronds and the soil surface observed in the field. Soil fauna abundance, especially of earthworms, and biogenic macroaggregates are lower in oil palm plantations as compared to native savanna and in improved pastures [START_REF] Lavelle | Soil ecosystem services and land use in the rapidly changing orinoco river basin of colombia[END_REF].

Consequently, fronds' organic matter does not benefit from the protection mechanisms that minerals would provide if it was integrated into the soils by fauna activity [START_REF] Schmidt | Persistence of soil organic matter as an ecosystem property[END_REF]. Application of composted mill residues would be a solution to decrease C to nutrient ratios of litter and improve the integration of organic matter into the soil by favouring soil fauna activity, humification and mineral protection.

The availability of organic residues from the palm oil production chain that could be applied in plantations, however, is limited. Leguminous cover crops, as already implemented in the younger plantations, are an alternative to increase aboveground and belowground C and N inputs. Enhancing root growth by mimicking the effects of frond piles is a promising solution to increase SOC if palms have plasticity to allocate more C to their rooting system when soil conditions are favourable or if palm varieties are developed for that purpose. Oil palm fronds could be spread on larger surface area or mixed with other residues, such as empty fruit bunched or fibers, to make a mulch. Even solutions not based on organic matter that would limit soil evaporation, retain nutrients, limit run-offs and erosion might improve root development. Future research should address whether SOC stocks eventually recover to initial SOC levels. Ensuring long-term soil fertility in oil palm plantations is fundamental to avoid a conversion of natural ecosystem constrained by soil degradation in older plantation.

The benefits from increasing SOC in terms of climate change mitigation would cascade far beyond the amount of C sequestrated in the soil because of the gain in soil fertility. Letters indicate significant differences between management zones for each parameter

TABLES

  package), followed by Tukey HSD tests for post hoc pairwise comparisons, to test for the effects of management zones (FP, HP, IZ and WC) on soil parameters (bulk, C3 and C4 stocks, fine root biomass, net C3 stabilization, basal respiration, microbial biomass C and C:N ratio, and metabolic quotient) at each soil depth (0-10, 10-20 and 20-30 cm) in the mature plantation (9-year old) with palm trees as random factor. Normal distribution of residuals and homogeneity of variance were tested by Shapiro and Levene tests, respectively, and data was log-transformed if necessary. Causal relationships between parameters (C stocks vs. age, C3

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .

 12345 Fig. 1. Management zones in a mature oil palm plantation. (a) Four management zones

Fig. 6 .Fig. 7 .Fig

 67 Fig. 6. Pearson correlation matrix among soil variables in the top 10 cm. Correlation

Table 1 .

 1 Oil palm-derived SOC stabilization under weeded circles at each depth depending on time, fine root biomass in the plantation and fine root biomass cumulated since plantation establishment (mean ± SE). All linear regressions were significant at p < 0.001.

	Weeded circles	OP-derived SOC	OP-derived SOC	OP-derived SOC
	Depths	per year		per fine roots		per cumulative fine roots
	cm	g SOC m -2 yr -1	R 2	g SOC g -1 roots	R 2	g SOC g -1 roots	R 2
	0-10 (n=14)	54  6	0.85	0.62  0.10	0.76	0.14  0.02	0.86
	10-20 (n=15)	18  6	0.55	0.53  0.11	0.60	0.11  0.02	0.59
	20-30 (n=15)	12  2	0.75	0.50  0.11	0.57	0.11  0.01	0.80

ACKNOWLEDGEMENTS

This study was financed by the Swiss National Science Foundation (r4d -Ecosystems) "Oil Palm Adaptive Landscape" no. 152019. Contribution of YK was supported by the Russian Science Foundation (project No. 18-14-00362). We thank the owners and agronomists of the Palmeras Sillatava and Samani Italcol de Occidente for granting us access to their land and for their logistical support.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.