Garcia Pinto

Luka Stanisic

Arnaud Legrand

Lucas Mello Schnorr

Samuel Thibault

Vincent Danjean

Analyzing Dynamic Task-Based Applications on Hybrid Platforms: An Agile Scripting Approach

In this paper, we present visual analysis techniques to evaluate the performance of HPC task-based applications on hybrid architectures. Our approach is based on composing modern data analysis tools (pjdump, R, ggplot2, plotly), enabling an agile and flexible scripting framework with minor development cost. We validate our proposal by analyzing traces from the fullfledged implementation of the Cholesky decomposition available in the MORSE library running on a hybrid (CPU/GPU) platform. The analysis compares two different workloads and three different task schedulers from the StarPU runtime system. Our analysis based on composite views allows to identify allocation mistakes, priority problems in scheduling decisions, GPU tasks anomalies causing bad performance, and critical path issues.

I. INTRODUCTION

To fulfill the ever-growing need for computation power of High-Performance Computing (HPC) applications, it has become common to rely on hybrid nodes, composed of multicore processors (CPUs) with multiple accelerators (GPUs). However, due to the heterogeneity and complexity of such machines, achieving portable and scalable performance has become extremely challenging. A possible solution, increasingly used by the community, is to program the application at a high level, independently of the hardware architecture, as a Directed Acyclic Graph (DAG) of tasks. It is then the responsibility of another software layer, called the runtime, to dynamically schedule the resulting tasks on the different computing resources taking into account the possible speed heterogeneity and variability as well as to automatically take care of data movements between resources. This allows removing artificial synchronizations from the application code, implementing complex scheduling and data movement algorithms (such as HEFT [START_REF] Topcuoglu | Performance-effective and lowcomplexity task scheduling for heterogeneous computing[END_REF]) that would be hard to manually incorporate in the application. Thanks to dynamic decisions, the irregular behavior of applications and resources is absorbed and exploited by the scheduler during execution, effectively balancing the load among computing resources.

Task-based executions on hybrid platforms are inherently stochastic. Task mapping, for instance, can drastically change from one execution to another. From the performance analysis perspective, the nature of such opportunistic execution schemes makes classical performance analysis, which expects regularity and highlights heterogeneity, invariably ineffective. At the same time, application and runtime developers seek to understand the attained performance to improve the application design and scheduling decisions. This investigation is challenging because it is hard to compare many traces when parameters affecting task creation and scheduling change.

In this article, we explain how we designed a framework that enables an easy and faithful identification of subtle scheduling problems that would otherwise go unnoticed and misunderstood with classical trace visualization approaches. We built the framework 1 on top of modern data analytics tools, combining the R programming language (and in particular the ggplot2 library) and org-mode [START_REF] Schulte | A multi-language computing environment for literate programming and reproducible research[END_REF]. The tool combination comes at a very low development cost when compared to a traditional and monolithic performance visualization tool. The designed views depict task execution along time for each resource, automatically detecting several interleaving critical paths in traces. We demonstrate the effectiveness of our visualization approach by analyzing traces from the dense linear algebra Cholesky factorization of the Chameleon/MORSE package [START_REF] Agullo | Poster: Matrices over runtime systems at exascale[END_REF], implemented using the StarPU task-based runtime [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF]. Two representative factorization workloads are carried out on a hybrid multi-core/multi-GPU architecture. The large workload brings interesting insights on pinpointing resource usage mistakes and comparing three StarPU schedulers (DMDA, DMDAS and Work Stealing). The smaller workload shows that the dynamic critical path analysis provides hints for optimizations.

Section II provides some background on task-based runtimes for hybrid platforms and on the Cholesky algorithm. Section III presents some related work on trace visualization, motivating our own study. Section IV presents our trace visualization proposal for the performance analysis of taskbased runtimes. In Section V we detail two case studies demonstrating the effectiveness of our visualization strategy. Section VI gives a summary of results and future work.

II. BACKGROUND AND EXPERIMENTAL CONTEXT

Traditional bulk-synchronous parallel (BSP) applications, made of supersteps (computation, communication, barrier), is a very common design when computing resources are homogeneous. The trend towards heterogeneous resources with accelerators, such as GPUs, makes the task-based programming paradigm much more suitable. Depending on the sophistication of the runtime scheduling algorithms, applications can more or less efficiently exploit heterogeneous configurations.

Several runtimes targeting hybrid platforms have been developed in the recent years. MAGMA [START_REF] Agullo | Numerical linear algebra on emerging architectures: The plasma and magma projects[END_REF], for instance, com-bines multi-core with GPUs for linear algebra applications. OmpSs [START_REF] Duran | OmpSs: a proposal for programming heterogeneous multi-core architectures[END_REF] provides an extension to OpenMP tasks through new directives that allow supporting multi-core systems combined with GPUs. PaRSEC [START_REF] Bosilca | DAGuE: A generic distributed DAG engine for high performance computing[END_REF] is a generic framework for architecture aware scheduling of tasks on many-core heterogeneous clusters. StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] is a task parallelism runtime initially designed to exploit hybrid architectures and additionally providing an MPI-based extension [START_REF] Augonnet | StarPU-MPI: Task programming over clusters of machines enhanced with accelerators[END_REF] to exploit several nodes at once. Dense/sparse linear algebra have been among the first applications to exploit such runtimes but other applications have been also implemented, such as FEM applications [START_REF] Ohshima | Implementation of FEM Application on GPU with StarPU[END_REF], seismic wave modeling [START_REF] Martínez | Towards seismic wave modeling on heterogeneous many-core architectures using task-based runtime system[END_REF], and others [START_REF] Agullo | Taskbased Conjugate Gradient: from multi-GPU towards heterogeneous architectures[END_REF], [START_REF] Lacoste | Taking advantage of hybrid systems for sparse direct solvers via task-based runtimes[END_REF].

Experimental Context and Workload Details: We propose to build on a relatively representative use case to address the lack of adequate visualization tools for task-based applications. We use traces of the dense linear algebra Cholesky factorization, more specifically of the full-fledged implementation from the MORSE library [START_REF] Agullo | Poster: Matrices over runtime systems at exascale[END_REF], compiled with the CUBLAS kernels. A simplified version of this application is shown in Figure 1a; and the corresponding DAG for a 5 × 5 matrix size is shown on its right (Figure 1b). For each step k of the outer loop, one dpotrf task releases N -k dtrsm and dsyrk tasks, followed by ≈(N -k) 2 /2 dgemm tasks. The dependencies indicate that simultaneous execution of several iterations is possible and that the iteration size decreases at the same time as k increases. So, the potential parallelism gets reduced as the algorithm advances (see the DAG of the figure). Finally, the task's execution time highly depends on the task (dpotrf, dtrsm, dsyrk, and dgemm) and resource types (CPU or GPU). We have used the CUBLAS without the dpotrf code. Hence, dpotrf tasks can only be run on CPUs.

for (k = 0; k < N; k++) {

DPOTRF(RW,A[k][k]);

for (i = k+1; i < N; i++) Since hybrid heterogeneous nodes motivate the development of task-based runtimes, we execute this Cholesky implementation on IdCin2, a machine with two 14-core Intel(R) Xeon(R) CPU E5-2697v3@2.6GHz and three NVIDIA Titan X. From this set of resources, only 25 CPU cores participate in the computation because it is generally more efficient to let StarPU dedicate one core to manage each GPU. StarPU provides several scheduling algorithms that exploit both the DAG structure (through critical-path based heuristics) and performance models. Here, we focus on three of them.

DTRSM(RW,A[i][k], R,A[k][k]); for (i = k+1; i < N; i++) { DSYRK(RW,A[i][i], R,A[i][k]); for (j = k+1; j < i; j++) DGEMM(RW,A[i][j], R,A[i][k], R,A[j][k]); } } (a)
The DMDA (Deque Model Data Aware) scheduler is a list scheduler, i.e., every time a resource is idle, if a task is ready, it will be scheduled on this particular resource. Such a scheduler therefore never leaves a resource idle on purpose, which ensures the well-known (2 -1/p) competitive ratio for homogeneous machines [START_REF] Graham | Bounds for certain multiprocessing anomalies[END_REF]. Deciding which ready task to select has a major influence in practice and the classical heuristic consists of prioritizing tasks based on the critical path. However, the critical path notion is dynamic and obtaining a proper estimation can be quite challenging. With heterogeneous computing resources, such prioritization is generally done with variants of the HEFT (Heterogeneous Earliest Finish Time) strategy [START_REF] Topcuoglu | Performance-effective and lowcomplexity task scheduling for heterogeneous computing[END_REF]. The DMDA algorithm is a very greedy heuristic that schedules tasks in the order they become available, without taking critical path priorities into account. However, it considers data transfer time between CPUs and GPUs as well as the relative performance of resources on each computation kernel when taking its decision. The DMDAS (Deque Model Data Aware Sorted) scheduler is similar to the first strategy, except that it sorts tasks by priority, which can be expensive when the number of tasks is large. It is therefore rather close to the original HEFT algorithm by respecting priorities and taking past scheduling decisions into account. Finally, the WS (Work Stealing) scheduler uses a queue per worker; new tasks are kept local by default. When a worker is idle, it steals tasks from the most loaded worker.

The next section presents related work on performance analysis for BSP and DAG-based HPC programming models. We detail the issues of current solutions, motivating our work.

III. RELATED WORK AND MOTIVATION

Despite the plethora of runtimes to execute task-based applications on heterogeneous resources (see Section II), there are very few established tools to conduct a proper taskaware analysis. Developers usually rely on BSP-based trace visualization tools, whose objectives are different, seeking unexpected heterogeneity where regular, homogeneous behavior is normal. Such tools are therefore unsuitable for visualizing task execution behavior since heterogeneity is the expected scenario for task-based applications. We briefly detail trace visualization strategies for BSP and DAG-based applications to better understand their differences. We also describe the design challenges of novel trace visualization techniques for DAG applications, as well as the typical questions usually raised during the analysis process.

A. Trace visualization analysis for BSP-based applications

Many tools exist to visualize traces from BSP-based applications. Most are focused on MPI applications. The common technique is based on Gantt charts, depicting each thread behavior along time. Behavior is drawn using colors to represent different thread states e.g., MPI operation. Message-passing is depicted with arrows from source to destination, sometimes annotated with the amount of transferred data. Vite [START_REF] Coulomb | Visual trace explorer ViTE[END_REF] is an OpenGL-based open-source tool that has such a view, capable to visualize large traces with thread states and communication. Since this tool relies on the semantic-free Paje language [START_REF] Schnorr | The Paje trace file format[END_REF], it can depict virtually any kind of traces. Paraver [START_REF] Pillet | Paraver: A tool to visualize and analyze parallel code[END_REF] is another open-source Gantt tool that allows filtering, zooming, and graphical trace aggregation. Its trace format supports many HPC programming models. Vampir [START_REF] Knüpfer | The Vampir performance analysis toolset[END_REF] is a proprietary tool with multiple views for OTF2 traces. It is more scalable than similar tools due to its distributed organization.

B. Visualizing task execution traces from DAG applications

There are a few tools to perform analysis and visualization of task-based executions. Typically, they are built with resources not designed for data analysis and rely on either nonscalable or non-scriptable strategies i.e., with mouse pointer interaction. DAGViz [START_REF] Huynh | DAGViz: A DAG Visualization Tool for Analyzing Task-parallel Program Traces[END_REF] offers a visual representation of task-based executions focusing on the DAG structure, which is retrieved using macros (translated to Cilk, Intel TBB or OpenMP) and presented in a hierarchical way. The resulting DAG can be folded/unfolded on-demand to show details and the node color indicates where they are executed. There is no way to retrieve the time dimension and task duration, which can make performance analysis difficult. Temanejo [START_REF] Keller | Temanejo: Debugging of Thread-Based Task-Parallel Programs in StarSS[END_REF] provides a similar timeless DAG interactive views for many taskbased runtimes. Its main features are the debug capabilities, e.g. the user can put a breakpoint in a task and fix its dependencies. These capabilities are very useful during algorithm design on small scale, being unsuited for performance analysis Haugen et al. proposes an interactive Gantt chart enhanced with dependencies, drawn as edges between tasks [START_REF] Haugen | Visualizing execution traces with task dependencies[END_REF]. We believe this approach suffers from three issues. First, in term of scalability, since (e.g., in Cholesky) tasks typically have many (≈N) dependencies, drawing everything and finding interesting tasks and dependencies only through mouse interaction can be very tedious. In practice, only tasks belonging to the critical path are important. Second, only one-level dependencies are depicted, while several levels are required to understand the history leading to the scheduling problem. Third, this tool does not really exploit the heterogeneity of resources.

C. Challenges of DAG execution traces analysis

The performance analysis of task-based applications raises many challenges. The stochastic behavior of scheduling decisions induced by actual resource availability and by the task execution variability hinders the performance reproducibility. Moreover, executions are apparently unstructured, with no clear phases as in the BSP model. At the same time, task dependencies are part of the application and should be exploited by the analyst to understand bottlenecks.

Traces are generally much larger than the available screen space and naively displaying everything generally leads to biased views [START_REF] Schnorr | Visualizing more performance data than what fits on your screen[END_REF]. One, therefore, has to use selection schemes to show only data that is fully relevant from the analysis point of view. In such complex execution traces, many hypotheses (and thus filters) can be proposed regarding the expected behavior. Developing a monolithic tool, such as the ones tailored for BSP analysis, that anticipates all possible performance problems is thus impractical. Moreover, such tool would be quite difficult to maintain and customize. A more flexible solution with scripting capability is thus needed.

During the analysis of DAG traces, the typical questions that arise are different from those of BSP applications. Instead of inspecting whether and why a certain algorithm iteration or phase was slowed down, one needs to analyze if the tasks were properly scheduled and if sufficient parallelism is available. One also needs to examine whether complex data movement techniques, task submission, and many internal runtime mechanisms can be further improved.

A common approach to better understand the application behavior is to compare several execution traces, possibly each one with a different configuration (e.g., scheduling parameter). However, due to the dynamic scheduling, this is generally difficult for task-based traces. To draw relevant conclusions, one needs to synchronize multiple visualizations and filter the unwanted states. Although some support exists in some tools [START_REF] Pillet | Paraver: A tool to visualize and analyze parallel code[END_REF], [START_REF] Pagano | FrameSoC Workbench: Facilitating Trace Analysis through a Consistent User Interface[END_REF], they do not offer enough customization flexibility for such studies.

The next section details our proposal to visualize task-based application traces. It tackles some of the preceding challenges with a framework that combines data analytics tools to create a flexible environment enabling customized analysis.

IV. VISUALIZATION METHODOLOGY

Visualizing data allows to graphically check many assumptions at once. It helps with assumptions that are difficult to formally state or for which defining a proper statistical test would require building on even more hypothesis that would also have to be verified. This is why it is important to start by listing various hypothesis (or expectations) made on the system under study. From such list, a set of visualizations can be generated. We propose visualizations that are designed for the application and scheduler developers, assisting them to rapidly identify performance issues as well as potential solutions.

The set of hypothesis to check is fairly rich in heterogeneous platforms targeted by task-based runtime systems. It is thus important to build a visualization framework that allows to easily and rapidly combine various views and propose new alternative views in an agile way. Moreover, since dynamic scheduling and machine heterogeneity bring a lot of variabilities, the ideal visualization should exploit any potential regularity coming from the application algorithm. For example, as we have seen in Figure 1, each task can be identified by the loop indices i, j, k. Such kind of information is much more useful than the internal runtime task identification and should thus be provided by the application to the runtime so that it can be traced and further exploited during the visualization.

To meet these different design goals, we decided to use the workflow shown in Figure 2. With few modifications on MORSE, tracing is extended to tag tasks with loop indexes at the task creation. StarPU relies on FXT [START_REF] Danjean | An efficient multi-level trace toolkit for multi-threaded applications[END_REF] traces to produce timestamped events in the Paje language [START_REF] Schnorr | The Paje trace file format[END_REF]. A complete DAG is also created with task identifiers coherent with the Paje trace. Instead of building a complex monolithic tool, we follow the UNIX philosophy and script many small tools. Using pjdump, the paje trace is converted into a Comma-Separated Values (CSV) file that can be loaded into R. Thanks to the expressiveness and to the rich set of statistical libraries of the R language, many cleanups, filtering, statistic computations, and adaptation to another application than Cholesky can be done with few lines of code. The ggplot library provides a grammar of graphics and a very high-level way of building plots, enabling us to easily produce custom visualizations. This environment has enough expressiveness to guarantee different but coherent views (colors, scales, etc.). This approach allows building static views in a fully automatic and very efficient way. Although such visualizations could probably be sped up even further by programming everything in C/C++, the used libraries are already well optimized and benefit from the know-how of data analysts. Furthermore, a combination of small scripts is both easier to maintain and adapt to a new necessity or to a particular situation than a rigid monolithic visualization environment.

The static views (typically basic X11 window or a PDF) of our approach have disadvantages when compared to tools described in Section III. Interaction is often crucial for the analyst to find what he is looking for. This is why we also build on plotly, an online analytics tool, that enables the quick conversion of ggplots into interactive, online graphs usable with a classical web browser (see https://plot.ly/ggplot2/ for more details). Two illustrations in this article are also available in an interactive version. We strongly believe that putting interaction at the very end together with the scripting capabilities in the core of the analysis process is the key to carrying out the analysis of complex execution traces.

Finally, loading and merging several traces enables faithful comparisons and even produces perfectly aligned and coherent views. Since execution traces are stochastic in nature, we believe that side by side representations are essential to decide whether a phenomenon is important/recurrent or it is an unlucky situation and can be considered as the execution noise.

V. EXPERIMENTAL VALIDATION AND RESULTS

When dynamically scheduling task-based applications, the kind of question to answer is totally different whenever the resulting DAG is large or small. Large DAGs are expected to be embarrassingly parallel, almost reaching peak performance. Since such DAGs have many tasks, one need to use macroscopic views and indicators to understand how performance can be improved. Small DAGs, on the other hand, have little parallelism. Idle time will inevitably be incurred by task dependencies. For such executions, microscopic views with fine-grained data on task dependencies should rather be used.

We fully analyze two very different Cholesky workloads in the rest of this paper: large (60×60 tiles of size 960×960) and small matrices (12 × 12 tiles, same size). For each, we detail the expected behavior, then propose composite views allowing to check these expectations. The views are then exploited to compare the three aforementioned StarPU schedulers, enabling us to propose potential improvements.

A. Large Workload (Cholesky of 60×60 tiles of size 960×960) 1) Expectations: Uniformity. Task duration is expected to depend solely on their type (dgemm, dsyrk, dtrsm or dpotrf) and on the type of resource (CPU or GPU) on which it is executed. Such assumption should be visually verified, highlighting all tasks whose duration is abnormally large compared to the others of the same type/resource. We treat these tasks as independent outliers, unrelated to other tasks behavior. If not so, it may mean that the whole platform has been perturbed at particular moments or that some resource differs from the others. A task is anomalous if its duration exceeds the sampled third quartile plus 1.5 times the sampled interquartile range. Although this outlier notion is highly debatable and contextspecific, other definitions could be easily incorporated. Dependency problems. Large input matrices generate many tasks, especially when the application starts. We, therefore, want to monitor the number of ready and submitted tasks. For this Cholesky implementation, all tasks are expected to be submitted when the application starts. On scale, the number of task dependencies is extremely large. Automatically selecting which ones to display is haphazard. If a detailed view becomes necessary, we should switch to views described in Section V-B. Progress. The task graph resulting from dense linear algebra always share a common structure (for instance, see Figure 1). In a classical semi-sequential execution, the DAG would be executed much similarly to a breadth-first search. However, it is also possible to carry out a depth-first traversal, favoring task execution on the critical path. Following the pipelining of the sets of tasks submitted by each outer loop iteration can be sufficient to get an overview of how the scheduler is handling the DAG and if it corresponds to the analyst's intuition or not. Potential improvements. Dependencies are expected to be easily handled with large workloads. The major issue is the load balancing among CPUs and GPUs. Since one knows the average time w i,k needed to perform a task of type k on a resource of type i as well as the total number n k of tasks per type, one can consider that a fraction α i,k of tasks of type k will be done on resource i and that the α i,k should thus verify: ∀i :

k α i,k .n i,k .w i,k ≤ T
Since such constraints are linear it is possible to compute the optimal makespan T and allocation α i,k . The T value is called the Area Bound Estimation (ABE) and is a lower bound for the execution time. Another classical lower bound is the Critical Path Bound Estimation (CPE). It is obtained by assigning each task to its faster processing resource and by summing all durations along the DAG. These execution time bounds, in particular, the ABE when the workload is large, are quite useful to estimate how much further improvement can be expected. More accurate lower bounds [START_REF] Agullo | Bridging the Gap between Performance and Bounds of Cholesky Factorization on Heterogeneous Platforms[END_REF] could be used as well, in particular for intermediate size workloads. Moreover, an ideal task allocation is also computed when ABE is defined. Comparing the ideal with the actual allocation may help understanding how scheduling could be improved.

2) Composite View for Pinpointing Scheduling Mistakes: Building on the expectations, we propose a five-area composite view shown in Figure 3 (each area is detailed in the caption).

We see in Figure 3 (a) that the makespan is 62725ms while the ABE is 59464ms. So, one can hope for a 5% improvement. The scheduling seems indeed inefficient since there are periods (white areas in CPUs) when no useful computation is done. These periods correspond to filtered states (for clarity) where threads try to actively fetch data. The total idle state (yellow areas: a summary in the right of the Gantt) for CPUs is about 1%, while for GPUs it ranges from 2 to 6%. This GPU inactivity is likely the main source of potential improvement. From (c), it is clear that this idle time does not come from a sudden lack of ready tasks. Figure 3 (d) clearly indicates that all tasks have been submitted in the beginning and that task execution started immediately after, without waiting for fully unrolling the DAG. As suggested in (b), DAG traversal is rather depth-first. Many outer loop iterations are parallel (the maximum is 30 around 40s), explaining why there are always a sufficient number of ready tasks.

Such GPUs starvation is more likely explained either by data prefetching problem (some tasks are ready but their input data is not yet transferred to GPUs) or possibly by some priority issue (the priorities, used by the scheduler to choose which task to run first when several of them are ready, might be inadequate). The first reason is likely to be the right one here. Indeed, most large idle periods on GPUs and large periods of times where CPUs are not doing useful computations (in white) also coincide with abnormal dgemm tasks (in dark green) on GPUs. Note that the scripting feature allow us to show only tasks of a given type, resource, within a specific time interval, and whose duration is larger than a threshold and therefore clear any doubt the analyst may have. An investigation (such as the one of Section V-B) reveals that, for an unknown reason, the GPUs seem to freeze during a task execution inside the proprietary CUBLAS dgemm kernel, ultimately blocking tasks eagerly waiting for GPU data. Understanding why GPUs sometimes get stuck would certainly solve the issue but this clearly suggests a weakness of the chosen scheduler which assumes that tasks duration have small variability. Using other schedulers may, therefore, alleviate this.

The four plots depicted in Figure 3 (e) show the ideal allocation when calculating the ABE. They show how the GPUs have been overused with dgemm tasks and under-exploited for dsyrk and dtrsm tasks. It, therefore, suggests constraining the dsyrk and dtrsm tasks to run exclusively on GPUs.

3) Comparing Scheduling Strategies and Task Constraints: The previous analysis leads us to vary the scheduler (DMDA, DMDAS, WS) and to force or not the dsyrk/dtrsm allocation on GPUs. Figure 4 provides the six-scenario comparison.

First of all, it is interesting to see how the three schedulers differ in their DAG traversal of the DAG. While the DMDA algorithm has a breadth-first traversal (very few iterations of the outer loop are active at the same time), the DMDAS has a much more depth-first traversal as it takes the priority of the critical path into account. The traversal of the Work Fig. 4: The execution representation of three schedulers (columns), with unconstrained versus constrained dsyrk/dtrsm tasks on GPUs(rows). Each of the six plots is an instance of the view of Figure 3 (refer to that caption to understand different parts).

Work Stealing Constrained                                            Unconstrained                                          
Stealing (WS) is even more depth-first as almost all outer loop iterations are still in progress at the end of the execution. Such way of progressing through the DAG is typical of WS and somehow favors local data accesses even though the algorithm is more dependency myopic than the two other ones.

Second, when constraining the dsyrk and dtrsm to run solely on the GPUs (the plots on the bottom row of Figure 4), task allocation corresponds to the ideal one. However, while such constraint allows both DMDAS and Work Stealing to obtain near optimal executions (within less than 2% of the lower bound), it helped only moderately the DMDA algorithm. Many synchronized idle phases can be observed and imputed to both dependency issues (not enough parallelism is obtained from such a strict breadth-first traversal) and particularly slow tasks (probably slowed down by simultaneous data transfers). Interestingly, very few outlier tasks appear in the DMDAS and WS executions although the latter still seems a bit sensitive to this, as inactivity periods on CPUs (white areas) still correlate with the occurrence of dgemm outliers (darker green) on GPUs.

Finally, we stress that such observations are no coincidence. We randomly ran similar scenarios ten times and although the numbers always slightly differ, the general behavior and conclusions are the same. We also highlight that the area bound estimations (ABE) can vary significantly between two scenarios (e.g., 60s for constrained DMDA vs. 57s for constrained WS), which can be initially surprising since these estimates only depend on the number of tasks and their per-type average execution time on the different resources. The observations can be explained by the use of sample execution time mean, which may vary a bit. From our investigation, this variation is not explained by outliers occurrence but rather biased toward one or another scheduler. We think this is the consequence of a better locality (cache usage) but more complex measurements would be needed to fully evaluate this hypothesis.

B. Small Workload (Cholesky of 12×12 tiles of size 960×960) 1) Expectations: Potential improvements. The area bound (ABE) is optimistic for small workloads since it ignores task dependencies. The critical path bound (CPE) is much more relevant, especially on very small workloads such as this one. Still, knowing how tight they are is quite difficult [START_REF] Agullo | Bridging the Gap between Performance and Bounds of Cholesky Factorization on Heterogeneous Platforms[END_REF]. Thus, comparing to the ideal CPU/GPU allocation is meaningless and we should focus mainly on filtering task dependencies. Idle time everywhere. Having a lot of idle time is expected because of dependencies. It is thus imperative to identify crucial tasks, highlighting the "dynamic critical path", i.e., the last tasks upon which they depended on. Such important tasks may appear either because of the DAG structure or because they have been particularly delayed.

2) Backtrack Dependencies to Pinpoint Scheduling Mistakes: Since StarPU also exports the DAG, the dependency information can be merged into the trace. For Cholesky, the dpotrf are critical tasks that release many other tasks. It seems, therefore, relevant to track their dependencies. For a given task T i , it is possible to compute what was the task T i-1 on which it depends on and that finished the latest, similarly, for T i-1 its latest predecessor T i-2 , etc. Such backward tracking of dependencies allows rebuilding the observed critical path. Given the DAG, in an ideal execution one would expect a dpotrf to be immediately preceded by a dsyrk, immediately preceded by a dtrsm that would, in turn, be immediately preceded by the dpotrf of the previous iteration. Therefore, we compute such filtered backward dependencies for each dpotrf and merge them together. This is illustrated in Figure 5. The makespan is 730ms while the ABE is 434ms and the CPE is 368ms. The bounds may be loose but it seems that there is room for improvement. If we start from the end of the schedule and go backward in time, we can see a dependency path (in blue) that, until timestamp 400ms, fully respects the alternation dpotrf-dsyrk-dtrsm. At the very end, all tasks execute right one after the other, which is optimal. The first "mistake" appears in time 600ms where the dsyrk could have been executed a little earlier. Slightly before, some dtrsm are not executed right after their dpotrf maybe because of data transfer or more likely because of a wrong priority. This critical path does not merge with the one obtained for the dpotrf of the first iterations. Now, when looking at the other (red) dependency path, we can see many times that the tasks are scheduled as soon as possible as if there was some priority problem, which could possibly be solved with another scheduler. The scripting feature allows plotting only dependencies whose duration is larger than a given threshold, avoiding graphical clutter.

We identify another problem with the blue dependency path. At the end, tasks are executed on the appropriate resource (dpotrf on CPUs, and dsyrk and dtrsm on GPUs). However, slightly before time 600ms, critical dsyrk tasks start running on the CPUs, slowing the progression. Likewise, slightly before time 400ms, critical dtrsm tasks are run on CPUs whereas they are known to be very slow on such resources. It seems that this scheduler makes a bad decision and that constraining dtrsm and dsyrk to be executed on GPUs may reduce the total makespan. Although the potential solutions (fix priority, constraining some tasks to GPUs) suggested by this analysis are the same as in the previous use case, the underlying reasons are fundamentally different.

3) Comparing Scheduling Strategies and Task Constraints: Based on the previous analysis, we vary again the three schedulers and forcing dsyrk/dtrsm tasks on GPUs. Figure 6 compares the six resulting combinations for this workload.

By taking a closer look at the unconstrained top row of Figure 6, we can observe that the behavior demonstrated by the DMDA and DMDAS schedulers are not so different. They both have similar runtime, two unmerged critical paths on which priority and critical task allocation problems can be identified. WS also demonstrates a very bad allocation, which is not surprising because it does not take into account the heterogeneity of the platform. There are three dynamic critical paths in the WS scheduler, with many dtrsm and dsyrk running on CPUs. When constraining these two task types to execute only on GPUs (the bottom row of Figure 6), we observe that such restriction does not really help for the DMDA and DMDAS schedulers. Tasks on the critical path are no longer an issue, but both schedulers still have priority problems. The behavior demonstrated by DMDA seems easier to understand: we see some typical list scheduling behavior with critical dpotrf being delayed because CPUs are used for not so critical dgemms. If one could run these tasks earlier, it appears that the whole makespan would be greatly improved.

Surprisingly, Work Stealing strongly benefits from the imposed restriction and now favorably compares against DMDA and DMDAS. It is also interesting to note that WS manages to keep all CPUs busy from the very beginning, unlike the other two schedulers. However, GPUs are not fully exploited, in particular at the end where they should be used to accelerate the dgemms like the DMDA and DMDAS strategies do. If there was a way to prevent dgemm task execution on CPU after time 350ms, we would probably get the best of the two scheduling strategies and be much closer to the optimal execution time.

VI. CONCLUSION

This article presents how an agile scripting framework allows creating faithful and enlightening trace views for the performance analysis of task-based HPC applications running on heterogeneous platforms. The framework follows the UNIX philosophy and builds upon many small existing tools (pjdump, ggplot2, plotly) glued together with the R programming language. Putting interaction at the very end together with scripting capabilities in the core of the analysis process is the key to a flexible design. We have shown how applying this framework in the analysis of large and small scale scenarios based on the MORSE/StarPU implementation of the Cholesky factorization allows to identify and address several non-trivial scheduling problems. The approach has scalability limitations mostly in the graphical representations (static or interactive) and not the scripting methodology for creating them. Plotly, for instance, is incapable of handling too much data, such as when there are many task types or resources.

As future work, we intend to improve the interactive graphical views by using temporal aggregation to reduce trace size and improve the integration between performance visualization Fig. 6: The execution representation (interactive version at http://perf-ev-runtime.gforge.inria.fr/vpa2016/) of three schedulers (DMDA, DMDAS and Work Stealing), as columns, with unconstrained versus constrained dsyrk and dtrsm tasks on GPUs, as rows, for a 12 × 12 tiled Cholesky decomposition. Each of the six plots is an instance of the view described in Figure 5.

                     Unconstrained                     
and runtime debugging. We also intend to extend this approach to analyze the behavior of irregular task-based workloads [START_REF] Agullo | Multifrontal QR factorization for multicore architectures over runtime systems[END_REF].

  Corresponding DAG for N = 5.

Fig. 1 :

 1 Fig. 1: The Cholesky code and its DAG (for N = 5).

Fig. 2 :

 2 Fig. 2: Combining data analytics tools to create trace views.

Fig. 3 :

 3 Fig. 3: Cholesky factorization of a large (60 × 60 tiles of 960 × 960) matrix with the DMDAS scheduler; five views are referred as indicated on the right. (a) This Gantt chart is automatically augmented with (from left to right), the critical path bound (CPE) and of the area bound estimations (ABE) and the makespan. The time percentage in the Idle state per resource is shown on the right. Finally, tasks whose duration is abnormally large compared to the others are highlighted with a darker color. (b) The plot is composed of horizontal segments, one per iteration k of the outer loop, indicating during which time interval the tasks of the iteration have been processed. (c) depicts how many tasks are finished waiting for dependencies, and thus ready for execution and (d) how many are submitted and still not executed. (e) compares the actual number of tasks of each type executed on CPUs and GPUs (the bar) with the optimal repartition obtained when computing the ABE (the bullet).

Fig. 5 :

 5 Fig. 5: Detailed view (see http://perf-ev-runtime.gforge.inria. fr/vpa2016/) of a 12 × 12 Cholesky execution (DMDAS scheduler) with two (red and blue) critical paths.

  

	CPE 2201	ABE 59748	66965	0.7% 0.9% 1.0% 0.9% 0.9% 1.0% 1.0% 0.8% 1.0% 1.0% 1.0% 0.9% 0.9% 1.0% 1.3% 1.3% 1.2% 1.3% 1.4% 1.4% 1.6% 1.5% 1.6% 1.4% 1.6% 20.6% 20.2% 19.9%	CPE 2149	ABE 59464	62725	0.4% 0.6% 0.6% 0.7% 0.9% 1.0% 1.0% 0.9% 1.0% 1.0% 1.0% 0.9% 0.9% 0.9% 1.0% 1.0% 1.0% 0.9% 1.0% 1.1% 1.0% 1.1% 1.0% 1.0% 1.0% 5.9% 1.9% 2.0%	CPE 2146	ABE 58452	60987	1.1% 1.3% 1.2% 1.3% 1.3% 1.5% 1.4% 1.4% 1.5% 1.5% 1.3% 1.4% 1.2% 1.3% 1.5% 1.5% 1.5% 1.5% 1.4% 1.5% 1.5% 1.5% 1.4% 1.4% 1.5% 4.0% 2.2% 2.2%

Code available at http://perf-ev-runtime.gforge.inria.fr/vpa2016/

Acknowledgements: This work was supported by the projects: CAPES/Cofecub 764-13, FAPERGS/Inria ExaSE, CNPq 447311/2014-0, HPC4E 689772 from EU/H2020 and MCTI/RNP-Brazil, and the CNRS/LICIA Intl. Lab. We also thank Emmanuel Agullo, Lionel Eyraud-Dubois, and Suraj Kumar for their help in better understanding the traces and the visualization needs of the runtime scheduler developers.


Olivier Pourret
email: olivier.pourret@unilasalle.fr

David Houben

Characterization of metal binding sites onto biochar using rare earth elements as a fingerprint

The ability of biochar to immobilize metals relies on the amount of functional groups at its surface but the contribution of each functional groups (e.g. carboxylic, phenolic) to metal bonding is poorly known. Using a new approach based on previous works on rare earth element (REE) interactions with humic substances, we aim at elucidating the relative contribution of these binding sites to metal sorption under various conditions (i.e. pH and ionic strengths, IS). Using batch experiments, REE sorption onto biochar was analyzed from pH 3 to 9 and IS 10 -1 mol/L to 10 -3 mol/L. Rare earth element patterns show a Middle REE (MREE) downward concavity at acidic pH and low ionic strength. These patterns are in good agreement with existing datasets quantifying REE binding with humic substances. Indeed, the MREE downward concavity displayed by REE-biochar complexation pattern compares well with REE patterns with various organic compounds. This similarity in the REE complexation pattern shapes suggests that carboxylic groups are the main binding sites of REE in biochar.

Overall, our results indicate that the strength of the metal bonding with biochar increases when pH and IS increase, suggesting that biochar is more efficient for long-term metal immobilization at near neutral pH and high ionic strength.

Introduction

Biochar, the carbon-rich product of biomass pyrolysis, has sparked a broad interest due to its attractive potential for improving soil fertility, carbon sequestration and soil remediation [START_REF] Beesley | A review of biochars' potential role in the remediation, revegetation and restoration of contaminated soils[END_REF]Houben et al., 2013a;[START_REF] Lehmann | Biochar for Environmental Management[END_REF]. Surface functional groups play an important role in the application of biochar as a functional material, especially for the immobilization of toxic metals [START_REF] Chen | Quantification of Chemical States, Dissociation Constants and Contents of Oxygen-containing Groups on the Surface of Biochars Produced at Different Temperatures[END_REF][START_REF] Liu | Development of Biochar-Based Functional Materials: Toward a Sustainable Platform Carbon Material[END_REF]. For example, using biochars produced at 400 and 700 °C from rice straw and cattle manure, [START_REF] Qian | Dual Role of Biochars as Adsorbents for Aluminum: The Effects of Oxygen-Containing Organic Components and the Scattering of Silicate Particles[END_REF] found that aluminum (Al) binding onto biochar was predominantly due to complexation of Al with the hydroxyl and carboxyl group at the surface of biochar particles.

Similarly, the sorption of copper (Cu) and lead (Pb) was found to be directly correlated with the amount of oxygen-containing carboxyl, hydroxyl, and phenolic functional groups at the surface of biochars produced at various temperatures (200 -800°C) from cottonseed hulls [START_REF] Uchimiya | Screening biochars for heavy metal retention in soil: Role of oxygen functional groups[END_REF]. Understanding biochar surface reactivity is therefore a prerequisite to its application in soils. Spectroscopic approaches including Fourier transform infrared spectroscopy (FTIR) and nuclear magnetic resonance spectroscopy have been frequently employed to characterize the functional groups of biochar [START_REF] Brewer | Criteria to Select Biochars for Field Studies based on Biochar Chemical Properties[END_REF][START_REF] Keiluweit | Dynamic Molecular Structure of Plant Biomass-Derived Black Carbon (Biochar)[END_REF][START_REF] Zhao | Heterogeneity of biochar properties as a function of feedstock sources and production temperatures[END_REF]. However, these techniques do not allow determining to what extent each functional group contributes to metal binding under various conditions. Previous research on humic substances has shown that metal binding to carboxylic and phenolic sites was strongly dependent on pH and ionic strength in the aqueous phase [START_REF] Guo | Adsorption of metal ions on lignin[END_REF][START_REF] Ravat | Metal Ions Binding to Natural Organic Matter Extracted from Wheat Bran: Application of the Surface Complexation Model[END_REF][START_REF] Tipping | Cation binding by humic substances[END_REF]. Characterization of humic substances has also revealed that carboxylic groups generally exhibit low metal affinity sites, while phenolic groups host high affinity sites [START_REF] Benedetti | Metal Ion Binding to Humic Substances: Application of the Non-Ideal Competitive Adsorption Model[END_REF][START_REF] Kinniburgh | Ion binding to natural organic matter: competition, heterogeneity, stoichiometry and thermodynamic consistency[END_REF]. Unraveling the relative contribution of carboxylic and phenolic groups to metal binding by biochar under various pH and ionic strengths is therefore crucial to gain insight into the strength of association between metal and biochar and to predict metal interaction with biochar over a wide range of conditions.

An emergent approach for characterizing the surface reactivity of organic substances consists in analyzing the pattern of rare earth element (REE) sorption by these compounds [START_REF] Pédrot | Dynamic structure of humic substances: Rare earth elements as a fingerprint[END_REF]. Rare earth elements are a group of seventeen elements presenting very similar and coherent chemical properties [START_REF] Connelly | Nomenclature of inorganic chemistry -IUPAC recommendations[END_REF]. In this contribution, we will consider the fifteen lanthanides as REEs. Due to the unique properties of REEs, their partition coefficients (Kd; commonly used in estimating the potential sorption of dissolved metals by a solid phase) between humic substances and aqueous solution are highly sensitive to the heterogeneity of binding sites (i.e. the nature of the functional group involved) [START_REF] Marsac | Metal loading effect on rare earth element binding to humic acid: Experimental and modelling evidence[END_REF][START_REF] Sonke | Lanthanide-humic substances complexation. I. Experimental evidence for a lanthanide contraction effect[END_REF]. As a result, specific Kd patterns arise depending on the nature of the major binding sites at the surface of humic substances. For example, log Kd pattern exhibits a medium REE (MREE) downward concavity when REE binding occurs at carboxylic sites [START_REF] Pourret | Rare earth elements complexation with humic acid[END_REF]. On the other hand, REE binding to phenolic groups results in a log Kd pattern exhibiting a regular increase from light REEs (LREEs) to heavy REEs (HREEs) [START_REF] Davranche | Biogeochemical Factors Affecting Rare Earth Element Distribution in Shallow Wetland Groundwater[END_REF][START_REF] Marsac | Metal loading effect on rare earth element binding to humic acid: Experimental and modelling evidence[END_REF]. More specifically, studies have revealed that LREEs are mainly bound by weak carboxylic sites at pHs ranging between 3 and 6, while HREEs are preferentially bound to strong multidentate sites at pH = 3 and to phenolic groups at pH = 6 [START_REF] Marsac | An improved description of the interactions between rare earth elements and humic acids by modeling: PHREEQC-Model VI coupling[END_REF]. These findings have raised the possibility to use REEs as a tool to determine which of the carboxylic, phenolic or chelate functional groups are implicated in the binding of cations to humic substances under various environmental conditions [START_REF] Davranche | Biogeochemical Factors Affecting Rare Earth Element Distribution in Shallow Wetland Groundwater[END_REF][START_REF] Marsac | An improved description of the interactions between rare earth elements and humic acids by modeling: PHREEQC-Model VI coupling[END_REF].

In the current study, a new approach based on previous works on REE interactions with humic substances, was developed to provide new insight into the nature of the sorption sites responsible for metal binding to biochar under various pH and ionic strengths.

Material and methods

Biochar characteristics

Biochar was obtained from VT Green (France) who uses an industrial pyrolysis reactor (Biogreen®, ETIA, France) and sewage sludge as feedstock. Operating conditions include a residence time in the reactor of 10 min and an end temperature of pyrolysis of 500 °C. The biochar was then rinsed with MilliQ water and ground to 0.05 mm particle size. The pH of the obtained biochar was 8.2 (biochar:water ratio of 1:10; w/v). Its electrical conductivity at 20°C was 0.009 mS/cm and its cation exchange capacity was 4.6 mmol c /kg (cobalthexammine method) and its specific surface area was 19 m 2 /g (BET method). The biochar was further acid-rinsed with HNO 3 at pH 2 to remove ash and avoid masking metal sorption by precipitation phenomenon and by metal complexation with inorganic ligands like carbonate [START_REF] Houben | Impact of biochar and root-induced changes on metal dynamics in the rhizosphere of Agrostis capillaris and Lupinus albus[END_REF]. The PZC value of the studied biochar is 7.81, which is consistent with previous studies that evidenced an average value of 8.1 and a range from 6.2 to 9.6 [START_REF] Lu | Characterization of sewage sludgederived biochars from different feedstocks and pyrolysis temperatures[END_REF]. The biochar surface functional groups were measured with FTIR: 0.7 mg of biochar was gently mixed with 100 mg of KBr and the mixture were pressed into a pellet. FTIR spectrum of the biochar was recorded using a Nicolet iS10 spectrophotometer (Thermo Fisher Scientific), scanning from 4000 to 400 cm -1 at a resolution of 1 cm -1 .

The following FTIR broad-band assignment was used to characterize the biochar [START_REF] Jindo | Physical and chemical characterization of biochars derived from different agricultural residues[END_REF][START_REF] Keiluweit | Dynamic Molecular Structure of Plant Biomass-Derived Black Carbon (Biochar)[END_REF]: 3400 to 3410 cm -1 , H-bonded O-H stretching vibrations of hydroxyl groups from alcohols, phenols, and organic acids; 2850 to 2950 cm -1 , C-H stretching of alkyl structures; 1580 to 1590 cm -1 , COO -asymmetric stretching; 1460 cm -1 , C-H deformation of CH 3 group; 1280-1270 cm -1 , O-H stretching of phenolic compounds; and 1000-1100 cm -1 , bending of Si-O stretching. The sewage sludge-derived biochar from this study have comparable properties to those from previous studies [START_REF] Lu | Relative distribution of Pb2+ sorption mechanisms by sludge-derived biochar[END_REF].

Experimental setup description

All chemicals used in this study were of analytical grade, and all the experimental solutions used were prepared with doubly deionized water (MilliQ system, Millipore). REEs-biochar complexes were prepared in glassware containers previously soaked in 10% HNO 3 for 48 h at 60 °C, then rinsed with MilliQ water for 24 h at 60 °C to remove all possible REE contamination sources. Synthetic REE solutions were prepared from nitrate REE standards (10 mg/L, Accu Trace Reference Standard).

Batch experiments were performed to study the pH and ionic strength effect on REE sorption by biochar as performed for humic acid (HA) by [START_REF] Pourret | Rare earth elements complexation with humic acid[END_REF]. REE (50 µg/L of each REE) and biochar (20 mg/L) were added together in solution, at various ionic strength (0.001 mol/L to 0.1 mol/L, adjusted by NaCl) and pH (from 3 to 9, adjusted by HCl and NaOH).

Prior to addition of biochar, the pH of the solution was adjusted to near 4. The initial hydroxide and carbonate concentrations were considered negligible and concentrations of LnOH 2+ and LnCO 3 2-were therefore also considered negligible. The pH was measured with a Mettler Toledo EL 20 pH-meter equipped with a LE407 electrode having an accuracy of 0.01 pH units. Solutions were agitated at 120 rpm for 48 h (equilibrium time determined from preliminary kinetic experiments) at room temperature (20 +/-2°C) to allow equilibrium and REES fractionation between aqueous solution and biochar suspension. Solution aliquots of 10 mL were sampled at the beginning of the experiment and at equilibrium state (i.e., after 48 h).

REEs sorbed by biochar were separated from the remaining dissolved REE by filtration through 0.2 µm acetate cellulose filters.

To check possible sorption of inorganic REE species onto the membrane of the filter or onto tube walls, inorganic REE solutions of known REE concentration were filtered for each experimental conditions and monitored. More than 99% of the REEs present in the solution were recovered in the filtrates. It thus demonstrates that no REE were retained onto the membranes or on the walls of the tubes used.

Sorption behavior was quantified by using the apparent partition coefficient (i.e., Kd), defined as the ratio between REE concentration associated with the biochar and REE concentration in the aqueous solution when the system is at equilibrium, following the equation below:

Kd = (µg REE sorbed ⁄ g biochar)/(µg/mL REE in solution) (1)
Rare earth element concentrations in the filtered samples were determined by Inductively Coupled Plasma Mass Spectrometry (ICP-MS; Thermo Fisher Scientific XSERIES2).

Quantification was carried out by external calibration (REE multi elemental standard solution from Accu Trace Reference, USA) and using indium (2.5 μg/L) as an internal standard. The instrumental accuracy was assessed by analyzing the SPS-SW2 certified reference material for measurement of elements in surface water (SpectraPur standards, Oslo, Norway). The analyses of real samples were carried out provided that the bias of measured concentration was <5% compared to the certified values.

Modeling

Humic ion-binding Model VI [START_REF] Tipping | Humic ion-binding model VI: An improved description of the interactions of protons and metal ions with humic substances[END_REF] was introduced into PHREEQC version 3.1.5 [START_REF] Parkhurst | Description of input and examples for PHREEQC version 3 -A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations[END_REF] for the whole REEs group [START_REF] Marsac | An improved description of the interactions between rare earth elements and humic acids by modeling: PHREEQC-Model VI coupling[END_REF]. This coupling allows the REE-HA binding description to be improved [START_REF] Marsac | An improved description of the interactions between rare earth elements and humic acids by modeling: PHREEQC-Model VI coupling[END_REF]. Indeed, HAs are described as 80 different sites that are expected to participate in the cation-HA binding. To simplify the REE speciation on the HA surface, these 80 sites can be gathered in several groups based on two main characteristics: the chemical nature of the ligands (i.e. carboxylic groups: CG, phenolic groups: PG and carboxy-phenolic groups: CPG) and their denticity (i.e. mono-, bi-and tridentate sites) [START_REF] Marsac | An improved description of the interactions between rare earth elements and humic acids by modeling: PHREEQC-Model VI coupling[END_REF][START_REF] Tipping | Humic ion-binding model VI: An improved description of the interactions of protons and metal ions with humic substances[END_REF]. Modeling parameters are given in Table 1.

Results and discussion

pH drives the REE sorption on biochar

The effect of pH on REE sorption was investigated in an initial pH range between 3 and 9.

Results are shown for the first REE (i.e., La) in Figure 1. pH appears to be the parameter which most strongly modified surface reactivity of biochars [START_REF] Lee | Characterization of Biochars Produced from Cornstovers for Soil Amendment[END_REF]. The proportion of REEs removed from solution increases with increasing pH. The removal capacity is minimal at acidic pH (i.e. below pH 5) and increases along with the increase of pH.

Irrespective of the ionic strength, Figure 1 showed a regular increase of La concentration sorbed to the biochar surface with pH. This general feature is well expressed for all REEs and follows general trend of metal sorbed to biochars [START_REF] Chen | Adsorption of copper and zinc by biochars produced from pyrolysis of hardwood and corn straw in aqueous solution[END_REF][START_REF] Liu | Removal of lead from water using biochars prepared from hydrothermal liquefaction of biomass[END_REF].

The REE sorption increase with increasing pH likely results from deprotonation of the biochar surface functional group (e.g., carboxylic) [START_REF] Harvey | Metal Interactions at the Biochar-Water Interface: Energetics and Structure-Sorption Relationships Elucidated by Flow Adsorption Microcalorimetry[END_REF][START_REF] Lu | Relative distribution of Pb2+ sorption mechanisms by sludge-derived biochar[END_REF] making them available for REE sorption.

Moreover, as displayed on Figure 2, log Kd values increase with pH and log Kd REE patterns highlighted a HREE enrichment (i.e., La/Sm<1 and Gd/Yb<1), especially for acidic pH and higher ionic strength. These patterns are in good agreement with the few existing datasets quantifying the binding of REEs with humic substances [START_REF] Marsac | Metal loading effect on rare earth element binding to humic acid: Experimental and modelling evidence[END_REF]. In particular, the HREE enrichment displayed by REE-biochar complexation pattern determined in this study compares well with results from REE sorption on carboxylic sites hosted by humic acid (HA), fulvic acid (FA) and acetic acid (e.g., [START_REF] Davranche | Biogeochemical Factors Affecting Rare Earth Element Distribution in Shallow Wetland Groundwater[END_REF].

Ionic strength effect on REE fractionation

Ionic strength also appeared to have a general effect on the sorption of REEs onto biochar.

Log Kd increases with decreasing ionic strength (Figure 2). However, if we look in details, some exception arises. Indeed, experiments performed at ionic strength of 10 -3 mol/L and pH 4 and 5 show smaller Kd than that performed at pH 4.5 and ionic strength of 10 -2 mol/L. In addition, experiment performed at pH 9 and ionic strength of 10 -1 mol/L shows larger Kd than those performed under ionic strength of 10 -2 mol/L. Moreover, fractionation between LREEs and HREEs increases with increasing ionic strength (Figure 3). As shown in Figure 3, log Kd(Gd/Yb) and log Kd(La/Sm) values vary by ca. 0.15 and 0.2 log units, reflecting fractionation of the REE patterns.

Previous studies on organic matter have shown that HREE enrichment patterns is a general feature of simple organic compounds possessing carboxylic functional groups [START_REF] Marsac | Metal loading effect on rare earth element binding to humic acid: Experimental and modelling evidence[END_REF]. Such carboxylic functional groups can be regarded as weak sites. They displayed log Kd(Gd/Yb) <1 and log Kd(La/Sm) <1. A similar comparative approach may provide information about the nature of strong sites that occur at high ionic strength; high affinity sites are generally thought to comprise phenolic, hydroxylic and/or multidentate carboxylic functions [START_REF] Uchimiya | Influence of pH, Ionic Strength, and Multidentate Ligand on the Interaction of CdII with Biochars[END_REF] displaying log Kd(Gd/Yb) <1 and log Kd(La/Sm) <1, illustrating a lanthanide contraction effect. Thus, at low ionic strength (i.e., 10 -3 mol/L; Figures 2a and3), as log Kd(Gd/Yb) is closed to 1, REEs are interpreted to be mainly sorbed by weak carboxylic groups. At high ionic strength (i.e., 10 -1 mol/L; Figure 2c and 3), as log Kd(Gd/Yb) decreased below 1, they are interpreted to be mainly bound to strong multidentate carboxylic and phenolic ones (Figure 5). Moreover, ionic strength dependence can be used to differentiate inner-sphere (ionic strength independent) and outer-sphere (ionic strength dependent) sorption mechanisms [START_REF] Uchimiya | Influence of pH, Ionic Strength, and Multidentate Ligand on the Interaction of CdII with Biochars[END_REF]). An inner sphere complex is a stable molecular entity resulting from the ionic or covalent bonding of a metal ion with a surface functional group. In the outer-sphere mechanism at least one water molecule remains between the surface functional group and sorbed ion [START_REF] Stumm | Aquatic Chemistry[END_REF]. Our results reveal thus that the strength of REE bonding with biochar increases with increasing ionic strength. These findings have implications for the use of biochar as a REE immobilizing agent as they indicate that REE immobilization by biochar is expected to be more efficient at high than at low ionic strength at acidic pH whereas at more alkaline pH opposite behavior happened as already observed for other metals [START_REF] Hadjittofi | Activated biochar derived from cactus fibres -Preparation, characterization and application on Cu(II) removal from aqueous solutions[END_REF][START_REF] Uchimiya | Influence of pH, Ionic Strength, and Multidentate Ligand on the Interaction of CdII with Biochars[END_REF][START_REF] Vithanage | Mechanisms of antimony adsorption onto soybean stover-derived biochar in aqueous solutions[END_REF]. For example, for La at pH 8, ionic strength increase from 0.001 mol/L to 0.1 mol/L results in removal decrease from 60% to 30% whereas at pH 5, the same ionic strength increase results in removal increase from 10% to 20%. It can be thus inferred that REE binding to biochar was not weak and potentially reversible.

Rare earth elements as a fingerprint

As previously discussed, REE patterns highlight a shift from LREEs to HREEs with increasing ionic strength and pH (Figures 2 and3). This suggests that (i) at low ionic strength and acidic pH, REEs are mainly sorbed by carboxylic groups, (ii) at high ionic strength and alkaline pH, REEs are mainly sorbed by carboxy-phenolic groups. Moreover, when cations are bound to biochars, changes may occur in the surface reactivity of the biochars themselves such as (i) the formation of multidentate complexes and (ii) the decrease of the net surface charge [START_REF] Uchimiya | Influence of pH, Ionic Strength, and Multidentate Ligand on the Interaction of CdII with Biochars[END_REF].

Previous studies investigating REE binding to humic acids (i.e. considered in this study as analogous material to biochar) reported two kinds of REE-HA patterns depending on the occurring metal loading. (i) At high metal loading, MREE downward concavity patterns are observed [START_REF] Pourret | Rare earth elements complexation with humic acid[END_REF][START_REF] Yamamoto | Systematics of stability constants of fulvate complexes with rare earth ions[END_REF]. (ii) At much lower metal loading, patterns displaying a lanthanide contraction effect are highlighted [START_REF] Sonke | Lanthanide-humic substances complexation. I. Experimental evidence for a lanthanide contraction effect[END_REF][START_REF] Stern | A capillary electrophoresis-ICP-MS study of rare earth element complexation by humic acids[END_REF]. The combined effect of metal loading and HA surface heterogeneity [START_REF] Davranche | Biogeochemical Factors Affecting Rare Earth Element Distribution in Shallow Wetland Groundwater[END_REF] is expressed by stronger but less abundant surface sites (i.e., carboxyphenolic groups) at lower metal loading, and weaker but more abundant surface sites (i.e., carboxylic groups) at higher metal loading. Comparison of experimental log Kd REE-biochar patterns with the pattern for the stability constants (log K) of 101 REE organic ligands [START_REF] Byrne | Comparative complexation behavior of the rare earths[END_REF] highlight that biochar surface is in the same range as these 101 organic ligands and can thus be considered as a group of discrete sites (Figure 4). The major sites identified on the biochar surface are the carboxylic and carboxy-phenolic groups (as also highlighted by FTIR), which can form monodentate or multidentate complexes [START_REF] Uchimiya | Influence of pH, Ionic Strength, and Multidentate Ligand on the Interaction of CdII with Biochars[END_REF]. Using this simplified representation, biochar binding properties can be compared with the binding properties of model organic ligands as already shown in literature [START_REF] Davranche | Biogeochemical Factors Affecting Rare Earth Element Distribution in Shallow Wetland Groundwater[END_REF][START_REF] Marsac | An improved description of the interactions between rare earth elements and humic acids by modeling: PHREEQC-Model VI coupling[END_REF][START_REF] Marsac | Metal loading effect on rare earth element binding to humic acid: Experimental and modelling evidence[END_REF]. Indeed, The REE patterns corresponding to the binding of REE with acetate (carboxylic group model ligand) exhibit a MREE downward concavity, whereas NTA (chelate group model ligand) exhibit a lanthanide contraction effect corresponding to both REE-HA patterns obtained at high and low metal loading, respectively.

Moreover, as illustrated by the compilation of organic ligands [START_REF] Byrne | Comparative complexation behavior of the rare earths[END_REF] and infrared spectroscopy results [START_REF] Gangloff | Characterization and evolution of dissolved organic matter in acidic forest soil and its impact on the mobility of major and trace elements (case of the Strengbach watershed)[END_REF], when the ligand is stronger (i.e., strong aromatic functional groups such as carboxy-phenolic and phenolic groups), HREEs are more strongly bound to the ligand compared to LREEs, as observed on the REE-biochar binding experiments presented in the present study. Moreover, as already proposed by [START_REF] Pourret | Modeling lanthanide series binding sites on humic acid[END_REF], literature compilation of Log K(La/Sm) organic ligand constants [START_REF] Byrne | Comparative complexation behavior of the rare earths[END_REF][START_REF] Martinez | Modeling of rare earth element sorption to the Gram positive Bacillus subtilis bacteria surface[END_REF][START_REF] Pourret | Rare earth elements complexation with humic acid[END_REF][START_REF] Sonke | Lanthanide-humic substances complexation. I. Experimental evidence for a lanthanide contraction effect[END_REF] as plotted against Log K(Gd/Yb) (Figure 4) evidenced that calculated log Kd values in this study are in the range of natural carboxylic acids, further illustrating the major role of carboxylic groups.

To confirm this hypothesis and to gain access to the REE speciation on each biochar site types, the same modeling approach as already performed for humic substances was used [START_REF] Marsac | An improved description of the interactions between rare earth elements and humic acids by modeling: PHREEQC-Model VI coupling[END_REF]. Modeling can account for changes in pH and ionic strength as already proposed by [START_REF] Alam | Application of surface complexation modeling to trace metals uptake by biochar-amended agricultural soils[END_REF] for cadmium and selenium. Considering biochar as analogous to humic substances, further allows to dissociate binding parameter of carboxylic and carboxy-phenolic sites. Such modeling also allows accessing the REE distribution onto each binding group (Figure 5; Table 2). At low ionic strength (i.e. 10 -3 mol/L) and pH 3, REEs are mostly bound to carboxylic groups (33% for La) whereas at higher ionic strength and pH 3, LREEs are bound to carboxylic groups (21%) and HREEs are bound to carboxylic (13% for Lu) and carboxy-phenolic (10 % for Lu) groups. At pH of 5 and irrespective of the ionic strength, LREEs are mostly bound to carboxylic (from 83% at 10 -3 mol/L to 57% at 10 -1 mol/L for La) whereas HREEs are bound to carboxylic (from 38% at 10 -3 mol/L to 28% at 10 - 1 mol/L for Lu) and carboxy-phenolic (from 50% at 10 -3 mol/L to 38% at 10 -1 mol/L for Lu) groups. At a near neutral pH, LREEs are bound to carboxylic groups (from 85% at 10 -3 mol/L to 71% at 10 -1 mol/L for La) and carboxy-phenolic groups (from 14% at 10 -3 mol/L to 19% at 10 -1 mol/L for La) and HREEs via carboxy-phenolic (from 55% at 10 -3 mol/L to 54% at 10 -1 mol/L for Lu), carboxylic (from 28% at 10 -3 mol/L to 18% at 10 -1 mol/L for Lu) and phenolic

(from 18% at 10 -3 mol/L to 25% at 10 -1 mol/L for Lu) groups (Figure 5; Table 2). As already discussed, modeling further indicates that REE immobilization by biochar is expected to be more efficient at high than at low ionic strength at acidic pH whereas at more alkaline pH opposite behavior happened as already observed for other metals [START_REF] Hadjittofi | Activated biochar derived from cactus fibres -Preparation, characterization and application on Cu(II) removal from aqueous solutions[END_REF][START_REF] Uchimiya | Influence of pH, Ionic Strength, and Multidentate Ligand on the Interaction of CdII with Biochars[END_REF][START_REF] Vithanage | Mechanisms of antimony adsorption onto soybean stover-derived biochar in aqueous solutions[END_REF]. It must be noted that modeling considering HA overestimates sorption of REE at pH 7 but overall is a good representation of experimental results. Eventually, the denticity of the complex is dependent on the metal loading. The number of groups coordinated to REEs increased with decreasing loading [START_REF] Marsac | An improved description of the interactions between rare earth elements and humic acids by modeling: PHREEQC-Model VI coupling[END_REF][START_REF] Marsac | Effect of loading on the nature of the REE-humate complexes as determined by Yb3+ and Sm3+ LIIIedge EXAFS analysis[END_REF]. The binding of REEs to blocked phenolic and carboxylic sites of HA suggests that REE binding to stronger sites can be assumed to be the chelator [START_REF] Kautenburger | Influence of metal loading and humic acid functional groups on the complexation behavior of trivalent lanthanides analyzed by CE-ICP-MS[END_REF]. They consist in a combination of phenolic and carboxylic groups [START_REF] Kautenburger | Influence of metal loading and humic acid functional groups on the complexation behavior of trivalent lanthanides analyzed by CE-ICP-MS[END_REF]. EXAFS results considering ytterbium complexation with HA synthesized at low and high metal loadings, show that at low metal loading, REEs are bound to HA through bi-ligand complexes without any chelation effect. In these conditions, REE act as a cation bridge between two organic molecules. At the opposite, at high metal loading, REEs are bound to HA via multi-carboxylic ligands [START_REF] Marsac | Effect of loading on the nature of the REE-humate complexes as determined by Yb3+ and Sm3+ LIIIedge EXAFS analysis[END_REF]. This modelling approach is a considerable advance in the sense that such modeling considering HA surfaces as analogous of biochar surfaces account for changes in pH, ionic strength, and metal-tosorbent ratio.

Implications for the use of biochar as a sorbent for metals

An increasing number of studies consider biochar as a potential sorbent to remove and immobilize metals in water and soil, respectively [START_REF] Ahmad | Biochar as a sorbent for contaminant management in soil and water: A review[END_REF]. However, little is known concerning the stability, and therefore the longevity, of the immobilizing processes when environmental parameters change (Houben et al., 2013b). Like any other organic matter, biochar can possess functional groups, namely carboxyl and phenolic groups [START_REF] Uchimiya | Screening biochars for heavy metal retention in soil: Role of oxygen functional groups[END_REF], which host low and high metal affinity sites, respectively. Here, our results showed that, despite a decrease in the amount of metal sorbed onto biochar, increasing ionic strength favors metal sorption on high affinity sites, as indicated by the higher relative contribution of carboxy-phenolic and phenolic groups. Moreover, our data demonstrated that rising pH from acidic to near neutral values increases both the amount of sorbed metal and the relative contribution of carboxy-phenolic and phenolic groups to metal sorption. These findings have great implications for the use of biochar as a sorbent or an immobilizing agent for metals in water and soil because they reveal that environmental conditions affect not only the amount of sorbed metal onto biochar but also the metal-biochar binding strengths and stabilities. Overall, the strength of the binding of metals with biochar increases when pH and ionic strength increase, indicating that biochar is more efficient for long-term metal immobilization at near neutral pH and high ionic strength. Eventually, considering surface complexation modeling of biochar comparable to humic substances is a promising direction even if some other alternative may be considered [START_REF] Alam | Application of surface complexation modeling to trace metals uptake by biochar-amended agricultural soils[END_REF].

Conclusions

The efficiency of biochar to sorb metals is critically dependent on the type of metal-binding sites at the surface of biochar particles. Although techniques such as FTIR provide a preliminary biochar characterization by quantifying these sites, they give no information about the change in surface reactivity of biochar once added to water or soil. Using a new approach based on previous works on REE interactions with humic substances, this study is, to our knowledge, the first attempting to elucidate the relative contribution of these binding sites to metal sorption onto biochar under various conditions (i.e. various pH and ionic strengths). Our findings indicate that environmental conditions affect not only the amount of sorbed metal onto biochar but also the metal-biochar binding strengths and stabilities. In particular, the strength of the metal binding with biochar increases when pH and ionic strength increase, which has a great implication for the use of biochar as a metal sorbent in the long run. More generally, our study highlights the great potential of REEs as a fingerprint for the characterization of metal binding sites onto sorbing matrixes, including biochar. Table 1 Model VI parameters for humic acid [START_REF] Tipping | Humic ion-binding model VI: An improved description of the interactions of protons and metal ions with humic substances[END_REF].

Table and figure captions

Table 2 Proportions of REE sorbed by carboxylic groups (CG), phenolic groups (PG), carboxy-phenolic groups (CPG) and total (TOT) for various experimental conditions obtained by modeling using PHREEQC-Model VI.

Figure 1

 1 Figure 1 Proportion of La sorbed to biochar as a function of pH for various ionic strength

Figure 2

 2 Figure 2 REE patterns for log Kd (REE sorbed onto biochar) for increasing pH values as a function of ionic strength (a) 10 -3 mol/L, (b) 10 -2 mol/L and (c) 10 -1 mol/L.

Figure 3

 3 Figure 3 Log Kd(Gd/Yb) versus Log Kd(La/Sm) highlighting change among REE

Figure 4

 4 Figure 4 Literature compilation of REE constants with (a) low molecular weight organic

Figure 5

 5 Figure 5 REE patterns of modelled log Kd (REE-humic acid) at (a) pH 5 and IS 10 -3 mol/L and (b) pH pH 7 and IS 10 -3 mol/L, (c) pH 5 and IS 10 -2 mol/L, (d) pH 7 and IS 10 -2 mol/L, (e) pH 5 and IS 10 -1 mol/L and (f) pH 7 and IS 10 -1 mol/L, illustrating evolution of REE

Acknowledgments

We are grateful to Le Vinh Hoang, Bich Ngoc Tran, Quoc Viet Pham, Petru Jitaru and Lyes Lakhal for help during analytical procedures, to VT Green and ETIA for providing the biochar samples.