A Benoit
email: anne.benoit@ens-lyon.fr

L.-C Canon

L Marchal
email: loris.marchal]@ens-lyon.fr

Non-clairvoyant reduction algorithms for heterogeneous platforms †

Keywords: scheduling, reduction, approximation algorithms, non-clairvoyant algorithms

We revisit the classical problem of the reduction collective operation in a heterogeneous environment. We discuss and evaluate four algorithms that are non-clairvoyant, i.e., they do not know in advance the computation and communication costs. On the one hand, Binomial-stat and Fibonacci-stat are static algorithms that decide in advance which operations will be reduced, without adapting to the environment; they were originally defined for homogeneous settings. On the other hand, Tree-dyn and Non-Commut-Tree-dyn are fully dynamic algorithms, for commutative or non-commutative reductions. We show that these algorithms are approximation algorithms with constant or asymptotic ratios. We assess the relative performance of all four non-clairvoyant algorithms with heterogeneous costs through a set of simulations. Our conclusions hold for a variety of distributions.

INTRODUCTION

Reduction is one of the most common collective operations, together with the broadcast operation. Contrarily to a broadcast, it consists in gathering and summarizing information scattered at different locations. A classical example is when one wants to compute the sum of (integer) values distributed over a network: each node owns a single value and can communicate with other nodes and perform additions to compute partial sums. The goal is to compute the sum of all values. Reductions have been used in distributed programs for years, and standards such as MPI usually include a "reduce" function together with other collective communications (see [START_REF] Pjesivac-Grbovic | Performance analysis of MPI collective operations[END_REF][START_REF] Thakur | Optimization of Collective communication operations in MPICH[END_REF] for experimental comparisons). Many algorithms have been introduced to optimize this operation on various platforms, with homogeneous [START_REF] Rabenseifner | Optimization of Collective Reduction Operations[END_REF] or heterogeneous communication costs [START_REF] Liu | An Approximation Algorithm and Dynamic Programming for Reduction in Heterogeneous Environments[END_REF][START_REF] Legrand | Optimizing the steady-state throughput of scatter and reduce operations on heterogeneous platforms[END_REF]. Recently, this operation has received more attention due to the success of the MapReduce framework [START_REF] Dean | MapReduce: Simplified data processing on large clusters[END_REF][START_REF] Zaharia | Improving MapReduce performance in heterogeneous environments[END_REF], which has been popularized by Google. The idea of MapReduce is to break large workloads into small tasks that run in parallel on multiple machines, and this framework scales easily to very large clusters of inexpensive commodity computers. In this framework, the computation is split into a series of phases consisting in map and reduce operations. We review similar algorithms and use cases in the related work section (Section 2).

Our objective in this paper is to compare the performance of various algorithms for the reduce operations in a non-clairvoyant setting, i.e., when the algorithms are oblivious to the communication and computation costs (the time required to communicate or compute). This models well the fact that communication times cannot usually be perfectly predicted, and may vary significantly over time. We would like to assess how classical static algorithms perform in such settings, and to quantify the advantage of dynamic algorithms (if any). We use both theoretical techniques (worst-case analysis) and simulations on a wide range of random distributions.

The rest of the paper is organized as follows. Section 2 reviews existing reduction algorithms and other related work. Section 3 and 4 describes four algorithms and shows that they are approximation algorithms. Section 5 presents simulated executions of the previous algorithms and compares their respective performance, using several different random distributions for costs. Finally, we conclude and discuss future research directions in Section 6.

RELATED WORK

The literature has first focused on a variation of the reduction problem, the (global) combine problem [START_REF] Bar-Noy | An optimal algorithm for computing census functions in message-passing systems[END_REF][START_REF] Bruck | Efficient global combine operations in multi-port message-passing systems[END_REF][START_REF] Van De Geijn | On global combine operations[END_REF]. Algorithmic contributions have then been proposed to improve MPI implementations and existing methods have been empirically studied in this context [START_REF] Ali | Advanced collective communication in Aspen[END_REF][START_REF] Ritzdorf | Collective operations in NEC's high-performance MPI libraries[END_REF]. Recent works concerning MapReduce either exhibit the reduction problem or highlight the relations with MPI collective functions. We describe below the most significant contributions.

Bar-Noy et al. [START_REF] Bar-Noy | Computing global combine operations in the multiport postal model[END_REF] propose a solution to the global combine problem: similarly to allreduce, all machines must know the final result of the reduction. They consider the postal model with a constraint on the number of concurrent transfers to the same node (multi-port model). However, the postal model does not capture varying degree of overlapping between computations and communications.

Rabenseifner [START_REF] Rabenseifner | Optimization of Collective Reduction Operations[END_REF] introduces the butterfly algorithm for the same problem, with arbitrary array sizes. Several vectors must be combined into a single one by applying an element-wise reduction. Another solution has also been proposed when the number of machines is not a power of two [START_REF] Rabenseifner | More Efficient Reduction Algorithms for Non-Power-of-Two Number of Processors in Message-Passing Parallel Systems[END_REF]. These approaches are specifically adapted for element-wise reduction of arrays. Van de Geijn [START_REF] Chan | Collective communication: theory, practice, and experience[END_REF] also proposes a method with a similar cost. In our case, the reduction is not applied on an array and the computation is assumed to be indivisible.

Sanders et al. [START_REF] Sanders | Two-tree algorithms for full bandwidth broadcast, reduction and scan[END_REF] exploit in and out bandwidths. Although the reduction does not require to be applied on arrays, the operation is split in at least two parts. This improves the approach based on a binary tree by a factor of two.

Legrand et al. [START_REF] Legrand | Optimizing the steady-state throughput of scatter and reduce operations on heterogeneous platforms[END_REF] study steady-state situations where a series of reductions are performed. As in our work, the reduction operation is assumed to be indivisible, transfers and computations can overlap and the full-duplex one-port model is considered. The solution is based on a linear program and produces asymptotically optimal schedules with heterogeneous costs.

Liu et al. [START_REF] Liu | An Approximation Algorithm and Dynamic Programming for Reduction in Heterogeneous Environments[END_REF] propose a 2-approximation algorithm for heterogeneous costs and non-overlapping transfers and computations. Additionally, they solve the problem when there are only two possible speeds or when any communication time is a multiple of any shorter communication time. In the homogeneous case, their solution builds binomial trees, which are covered in Section 3.

In the MPI context, Kielmann et al. [START_REF] Kielmann | MPI's reduction operations in clustered wide area systems[END_REF] design algorithms for collective communications, including MPI Reduce, in hierarchical platforms. They propose three heuristics: flat tree for short messages, binomial tree for long messages, and a specific procedure for associative reductions in which data are first reduced locally on each cluster before the results are sent to the root process. Pjesivac-Grbovic et al. [START_REF] Pjesivac-Grbovic | Performance analysis of MPI collective operations[END_REF] conduct an empirical and analytical comparison of existing heuristics for several collective communications. The analytical costs of these algorithms are first determined using different classical point-to-point communication models, such as Hockney, LogP/LogGP and PLogP. The compared solutions are: flat tree, pipeline, binomial tree, binary tree, and k-ary tree. Thakur et al. [START_REF] Thakur | Optimization of Collective communication operations in MPICH[END_REF] perform a similar study for several MPI collective operations and compare the binomial tree with the butterfly algorithm [START_REF] Rabenseifner | Optimization of Collective Reduction Operations[END_REF] for MPI Reduce. These works, however, do not provide any guarantee on the performance.

Finally, this problem has also been addressed for MapReduce applications. Agarwal et al. [START_REF] Agarwal | A Reliable Effective Terascale Linear Learning System[END_REF] present an implementation of allreduce on top of Hadoop based on spanning trees. Moreover, some MapReduce infrastructures, such as MapReduce-MPI † , are based on MPI implementations and benefit from the improvements done on MPI Reduce.

The design and the analysis of algorithms in dynamic context has already received some attention. The closest related work is probably [START_REF] Canon | Comparative Evaluation of the Robustness of DAG Scheduling Heuristics[END_REF], in which the authors study the robustness of several task-graph scheduling heuristics for building static schedules. The schedules are built with deterministic costs and the performance is measured using random costs. [START_REF] Benoit | Computing the throughput of probabilistic and replicated streaming applications[END_REF] studies the problem of computing the average performance of a given class of applications (streaming applications) in a probabilistic environment. With dynamic environments comes the need for robustness to guarantee that a given schedule will behave well in a disturbed environment. Among others, [START_REF] Canon | Evaluation and optimization of the robustness of DAG schedules in heterogeneous environments[END_REF] studies and compares different robustness metrics for makespan/reliability optimization on task-graph scheduling. Optimizing the performance for task-graph scheduling in dynamic environments is a natural follow-up, and has been tackled notably using the concepts of IC (Internet-based Computing) and area-maximizing schedules [START_REF] Cordasco | On scheduling DAGs for volatile computing platforms: Area-maximizing schedules[END_REF].

MODEL AND ALGORITHMS

We consider a set of n processors (or nodes) P 0 , . . . , P n-1 , and an associative operation ⊕. Each processor P i owns a value v i . The goal is to compute the value

v = v 0 ⊕ v 1 ⊕ • • • ⊕ v n-1
as fast as possible, i.e., to minimize the total execution time to compute the reduction. We do not enforce a particular location for the result: at the end of the reduction, it may be present on any node.

There are two versions of the problem, depending on whether the ⊕ operation is commutative or not. For example, when dealing with numbers, the reduction operation (sum, product, etc.) is usually commutative while some operations on matrices (such as the product) are not. The algorithms proposed and studied below deal with both versions of the problem.

We denote by d i,j the time needed to send one value from processor P i to processor P j . A value may be an initial value or a partial result. When a processor P i receives a value from another processor, it immediately computes the reduction with its current value. To model communication congestion, we assume that each processor can receive at most one result at a time. Without this assumption, we may end up with unrealistic communication schemes where all but one processors simultaneously send their value to a single target processor, which would exceed the capacity of the incoming communication port on the target processor. The communication costs are heterogeneous, that is we may well have different communication costs depending on the receiver (d i,j = d i,j ), on the sender (d i,j = d i ,j ) and non-symmetric costs (d i,j = d j,i ). Even though these costs are fixed, we consider non-clairvoyant algorithms that make decisions without any knowledge of these costs.

The computation time of the atomic reduction on processor P i is denoted by c i . In the case of a non-commutative operation, we ensure that a processor sends its value only to a processor that is able to perform a reduction with its own value. Formally, assume that at a given time, a processor owns a value that is the reduction of v i ⊕ • • • ⊕ v j , which we denote by [v i , v j ]; the processor may only send this value to a processor owning a value

[v k , v i-1 ] or [v j+1 , v k ],
which is called a neighbor value in the following.

During a reduction operation, a processor sends its value at most once, but it may receive several values. It computes a partial reduction each time it receives a value. Thus, the communication graph of a reduction is a tree (see Figure 1): the vertices of the tree are the processors and its edges are the communications of values (initial or partially reduced values). In the example, P 0 receives the initial value from P 1 , and then a partially reduced value from P 2 . In the following, we sometimes identify a reduction algorithm with the tree it produces.

We now present the four algorithms that are studied in this paper. The first two algorithms are static algorithms, i.e., the tree is built before the actual reduction. Thus, they may be applied for commutative or non-commutative reductions. The last two algorithms are dynamic: the tree is built at run-time and depends on the actual durations of the operations.

The first algorithm, called Binomial-stat, is organized with log 2 n rounds. Each round consists in reducing a pair of processors that own a temporary or initial data using a communication and a computation. During round k = 1, . . . , log 2 n , each processor i2 k + 2 k-1 (i = 0, . . . , 2 log 2 n -k -1) sends its value to processor i2 k , which reduces it with its own value: at most 2 log 2 n -k+1 processors are involved in round k. Note that rounds are not synchronized throughout the platform: each communication starts as soon as the involved processors are available and have terminated the previous round. We can notice that the communication graph induced by this strategy is a binomial tree [START_REF] Cormen | Introduction to Algorithms[END_REF]Chapter 19], hence the name of the algorithm. This strategy is illustrated on Figure 2(a).

The second algorithm, called Fibonacci-stat, is constructed in a way similar to Fibonacci numbers. The schedule constructed for order k, denoted by FS k (k > 0), first consists in two smaller order schedules FS k-1 and FS k-2 put in parallel. Then, during the last computation of FS k-1 , the root of FS k-2 (that is, the processor that owns its final value) sends its value to the root of FS k-1 , which then computes the last reduction. A schedule of order -1 or 0 contains a single processor and no operation. This process is illustrated on Figure 2(b). Obviously, the number of processors involved in such a schedule of order k is F k+2 , the (k + 2)th Fibonacci number. When used with another number n of processors, we compute the smallest order k such that F k+2 ≥ n and use only the operations corresponding to the first n processors in the schedule of order k.

The previous two schedules were proposed in [START_REF] Canon | Scheduling associative reductions with homogeneous costs when overlapping communications and computations[END_REF], where their optimality is proved for special homogeneous cases: Binomial-stat is optimal both when the computations are negligible in front of communications (c i = 0 and d i,j = d) and when the communications are negligible in front of computations (c i = c and d i,j = 0). Fibonacci-stat is optimal when computations and communications are equivalent (c i = c = d i,j = d). In the non-commutative case, both algorithms build a tree such that only neighboring partial values are reduced. In the commutative case, any permutation of processors can be chosen. For Fibonacci-stat, the two schedules of order 2 and 3 used in the recursive construction are highlighted in green and red.

Then, we move to the design of dynamic reduction algorithms, i.e., algorithms that take communication decisions at runtime. The first dynamic algorithm, called Tree-dyn, is a simple greedy algorithm. It keeps a slot (initially empty), and when a processor is idle, it looks into the slot. If the slot is empty, the processor adds its index in the slot, otherwise it empties the slot and starts a reduction with the processor that was in the slot (i.e., it sends its value to the processor that was in the slot, and the latter then computes the reduced value). It means that a reduction is started as soon as two processors are available. Since in the obtained reduction tree, any two processors may be paired by a communication, this can only be applied to commutative reductions.

Finally, Non-Commut-Tree-dyn, is an adaptation of the previous dynamic algorithm to noncommutative reductions. In this algorithm, when a processor is idle, it looks for another idle processor with a neighbor value (as described above). Now, we keep an array of idle processors rather than a single slot. If there is an idle neighbor processor, a communication is started between them, otherwise the processor waits for another processor to become idle.

WORST-CASE ANALYSIS FOR COMMUTATIVE REDUCTIONS

We analyze the commutative algorithms in the worst case, and we provide some approximation ratios, focusing on communication times. We let d = min i,j d i,j , D = max i,j d i,j , c = min i c i , and

C = max i c i .
Let us first recall results from [START_REF] Canon | Scheduling associative reductions with homogeneous costs when overlapping communications and computations[END_REF], in the context of identical communication costs (d = D) and identical computation costs (c = C). The duration of the schedule built by Fibonacci-stat at order k is d + (k -1) max(d, c) + c, and the number of reduced elements n is such that F k+1 < n ≤ F k+2 , where F k is the kth Fibonacci number.

In order to ease the following worst-case analysis, we first derive bounds on the minimum duration of any schedule depending on the number of elements to reduce and the minimum costs (n, c and d).

Lemma 1

The order k of a schedule built by Fibonacci-stat can be bounded as follows:

log 2 n log 2 ϕ -1 < k < log 2 n log 2 ϕ + 1,
where ϕ = 1+ √ 5 2
is the golden ratio and n is the number of elements to reduce.

Proof

We know that F k+1 < n ≤ F k+2 , where F k is the kth Fibonacci number. By definition of the Fibonacci numbers, we have

F k = 1 √ 5 (ϕ k -(1 -ϕ) k ). Therefore, 1 √ 5 (ϕ k+1 -(1 -ϕ) k+1 ) < n ≤ 1 √ 5 (ϕ k+2 -(1 -ϕ) k+2 ).
Since -1 < 1 -ϕ < 0, it follows that

1 √ 5 (ϕ k+1 -(1 -ϕ) 2 ) < n < 1 √ 5 (ϕ k+2 -(1 -ϕ) 3 )
as soon as k ≥ 1. We therefore have ϕ k+1 < √ 5n + (1 -ϕ) 2 and √ 5n + (1 -ϕ) 3 < ϕ k+2 . We easily derive the following inequalities:

log 2 n log 2 ϕ + log 2 √ 5 + (1-ϕ) 3 n log 2 ϕ -2 ≥-1 when n≥1 < k < log 2 n log 2 ϕ + log 2 √ 5 + (1-ϕ) 2 n log 2 ϕ -1 ≤1 when n≥1 . Finally, log 2 n log 2 ϕ -1 < k < log 2 n log 2 ϕ + 1.

Proposition 1

The time taken to reduce n elements using any reduction algorithm is larger than max(c, d) log 2 n and than min(c, d) log 2 n log 2 ϕ .

Proof

First, any reduction algorithm has a duration that is lower bounded by the duration of the Binomialstat algorithm on a scenario where each communication costs d and each computation costs c.

Because of the overlap between communication and computation, each step of this reduction has a duration of max(c, d), and the total duration of this reduction is max(c, d) log 2 (n) , hence the first bound.

Second, the optimal schedule when all costs (communication and computation) are equal to min(c, d) is obtained with Fibonacci-stat, and its length is greater than (k + 1) min(c, d). Using Lemma 1, we obtain the second lower bound on the schedule duration.

We are ready to state the approximation ratios for Binomial-stat and Tree-dyn. Recall that a λ-approximation algorithm is an algorithm whose execution time is polynomial in the instance size, and that returns an approximate solution, guaranteed to be, in the worst case, at a factor λ away from the optimal solution.

Theorem 1

Binomial-stat and Tree-dyn are C+D max(c,d) -approximation algorithms, and this ratio can be achieved when min(c, d) = 0.

Additionally, these algorithms are also C+D min(c,d) 1 + 1 log 2 n log 2 ϕ-approximation algorithms, where

ϕ = 1+ √ 5 2
is the golden ratio (log 2 ϕ ≈ 0.69), and this ratio can be achieved when c = d.

Proof

The proof is in two steps: we first prove the approximation ratios by exhibiting an upper bound on the duration of Binomial-stat and Tree-dyn, before showing scenarios that achieve these ratios.

Let us consider an algorithm that would perform reductions through log 2 n synchronized steps. At each step, all processors possessing an element group themselves in pair for sending and reducing their elements (they perform their reductions pairwise). At each step, half the elements are processed. This algorithm takes at most a time (C + D) log 2 n , because, in the worst case, one computation and one communication at each step may be of maximum duration. Binomial-stat is faster than this algorithm because it has the same structure for the reductions (the source and destination of each communication is the same), except that no synchronization is required between each global step. Indeed, communications happen as soon as possible. Tree-dyn is also faster than this algorithm because two available processors may start as soon as possible without delaying the remaining steps. Therefore, the time taken by Binomial-stat and Tree-dyn is not greater than (C + D) log 2 n . The comparison of this upper bound to each of the two lower bounds given by Proposition 1 gives two approximation ratios. The first approximation ratio is C+D max(c,d) . The second ratio writes

(C+D) log 2 n min(c,d) log 2 n/ log 2 ϕ ≤ C+D min(c,d) log 2 ϕ 1 + 1 log 2 n .
We now exhibit a scenario for which these ratios are achieved. Let d i,j = d and d j,i = D, for all i < j, c i = C when i is even and c i = c if i is odd. With both Binomial-stat and Tree-dyn, we consider that any processor P i sends its element to a processor P j such that i > j, which takes a total time (C + D) log 2 n . Indeed, we consider that during the first step, each machine with an odd index sends an element to a machine with an even index, and these latter machines are performing the computations. The optimal solution performs only computations of size c and communications of size d (from P i to P j , with j > i, and only machines with odd index perform computations).

When min(c, d) = 0, the optimal total time is max(c, d) log 2 n and the first ratio is achieved. With c = d, a solution with a Fibonacci tree using only small communications is optimal [START_REF] Canon | Scheduling associative reductions with homogeneous costs when overlapping communications and computations[END_REF]; furthermore, it completes in time (k + 1)d. By Lemma 1, we have k > log 2 n log 2 ϕ -1 and the ratio is Unfortunately, the first approximation ratio of the previous theorem does not hold for Fibonaccistat and Non-Commut-Tree-dyn:

• For Fibonacci-stat, consider the schedule in Figure 2(b). Even without computation costs, the number of communication steps is 4 because P 0 has to receive four messages and these communications must be sequential. If these four communications take a time D (and all other d i,j 's are d), the optimal solution may complete in a time 3d (see for instance Figure 2(a), for Binomial-stat), hence a ratio of 4 3 D d . We prove below (see Theorem 2) that Fibonacci-stat is asymptotically a max(C,D) max(c,d) 1 log 2 ϕ -approximation, where1 log 2 ϕ ≈ 1.44. • For Non-Commut-Tree-dyn, the number of steps may well exceed log 2 (n) as well.

Consider for instance that for a reduction on 8 processors, in the first step neighboring processors communicate together two by two (as depicted on Figure 3). However, P 4 and P 2 terminate slightly before P 6 and P 0 , and initiate a communication. Then P 6 and P 0 must wait until completion of this partial reduction, and then two more steps are required to complete the reduction. Here again, the worst-case ratio is 4 3 D d (with d i,j = D if i > j, and d i,j = d otherwise). We do not prove any approximation ratio for this algorithm, because it does not seem very fair to compare a non-commutative version of the algorithm with the optimal commutative solution. Additionally, this algorithm is also a max(C,D) min(c,d) -approximation.

Proof

The length of the Fibonacci-stat schedule in the worst case is bounded by D + (k -1) max(D, C) + C, with k the order of the Fibonacci schedule [START_REF] Canon | Scheduling associative reductions with homogeneous costs when overlapping communications and computations[END_REF]. By Lemma 1, this length is further bounded by D + log2 n log 2 ϕ max(D, C) + C. Recall, from Proposition 1, that max(c, d) log 2 n is a lower bound on the reduction time, which gives the first approximation ratio.

The length of the Fibonacci-stat schedule in the worst case is also bounded by (k + 1) max(D, C). The second lower bound of Proposition 1 is (k + 1) min(d, c), which gives the second approximation ratio.

Table I summarizes the ratios for the case when c ≤ d (the opposite case is symmetrical). Note that for all algorithms, the second approximation ratio is asymptotically tighter when 1 ≥ c d ≥ log 2 ϕ ≈ 0.69, which corresponds to the situation for which the overlap between computations and communications is significant, whereas the first ratio is tighter for a smaller c d ratio, corresponding to a low overlap between communications and computations.

Using these asymptotic ratios, one may wonder when Fibonacci-stat is expected to have a better performance than Binomial-stat and Tree-dyn. To this end, we consider the ratio between C and D that makes Fibonacci-stat the best strategy. This asymptotically occurs when min C D , D C ≥ Table I. Approximation ratios for commutative algorithms (we assume c ≤ d due to symmetry).

first ratio (low overlap) second ratio (high overlap)

Binomial-stat and Tree-dyn

C+D d C+D c 1 + 1 log 2 n log 2 ϕ Fibonacci-stat max(C,D) d 1 log 2 ϕ + C+D d 1 log 2 n max(C,D) c
communications (i.e., when the maximum computation cost is approximatively between 44% and 228% of the maximum communication cost). Note also that when c = d = C = D, the asymptotic approximation ratio of Binomial-stat is 2 log 2 ϕ ≈ 1.39, which nicely complements the asymptotic approximation ratio of 1 log 2 ϕ ≈ 1.44 for Fibonacci-stat in [START_REF] Canon | Scheduling associative reductions with homogeneous costs when overlapping communications and computations[END_REF].

SIMULATION RESULTS

In the first part of this section, we consider that the d ij and c i costs are distributed according to a gamma distribution, which is a generalization of exponential and Erlang distributions. This distribution has been advocated for modeling job runtimes [START_REF] Lublin | The workload on parallel supercomputers: modeling the characteristics of rigid jobs[END_REF][START_REF] Feitelson | Workload modeling for computer systems performance evaluation[END_REF]. This distribution is positive and it is possible to specify its expected value (µ d or µ c ) and standard deviation (σ d or σ c ) by adjusting its parameters. Further simulations with other distributions concur with the following observations and are presented at the end of this section. This suggests that our conclusions are not strongly sensitive to the distribution choice. All simulations were performed with an ad-hoc simulator on a desktop computer.

Cost dispersion effect. In this first simulation, we are interested in characterizing how the dispersion of the communication costs affects the performance of all methods. In order to simplify this study, no computation cost is considered (µ c = 0). The dispersion is defined through the coefficient of variation (CV), which is defined as the ratio of the standard deviation over the expected value (this latter is set to µ d = 1). The number of processors is n = 64 and the time taken by each method is measured over 1,000,000 Monte Carlo (MC) simulations (i.e., simulations have been performed with 1,000,000 distinct seeds).

On a global level, Figure 4 shows the expected performance of the four studied reduction algorithms with distinct CVs. When the heterogeneity is noticeable (CV greater than 1), the performance decreases significantly. In these cases, schedules are mostly affected by a few extreme costs whose importance depends on the CV. Additionally, the variability in the schedule durations is also impacted by the CV (i.e., two successive executions with the same settings may lead to radically different performance depending on the schedule).

Several observations can be made relatively to each method. As expected, Binomial-stat is similar to Tree-dyn for CVs lower than 0.1. In this case, the improvement offered by Tree-dyn may not outweigh the advantage of following a static plan in terms of synchronization. For CVs greater than 1, both static approaches perform equally with a similar dispersion. For all CVs, Treedyn has the best expected performance while Fibonacci-stat has the worst (which is expected since no computation cost is considered), and Non-Commut-Tree-dyn has the second best expected performance when the CV is greater than 0.3. Finally, when the CV is close to 10, all methods are equivalent as a single communication with a large cost may impact the entire schedule duration. In terms of robustness, we can see that Fibonacci-stat and Non-Commut-Tree-dyn are the two best methods for absorbing variations as their expected durations remain stable longer (until the CV reaches 0.3). This is due to the presence of idleness in their schedules that can be used when required.

Non-negligible computation costs. When considering nonzero computation costs, we reduce the number of parameters by applying the same CV to the computation and to the communication costs (i.e., σc µc = σ d µ d ). As Fibonacci-stat is designed for overlapping computations and communications, we characterize the cases when this approach outperforms Tree-dyn.

Figure 5(a) shows the improvement of Tree-dyn over Fibonacci-stat when varying the CV and the ratio µc µ d , which corresponds to the overlapping degree between computations and communications. The darker the color, the better is Fibonacci-stat compared to Tree-dyn. The contour line with value 1 delimits the dark area for which Fibonacci-stat is better than Tree-dyn on average. This occurs when the computation cost is greater than around half the communication cost and when the variability is limited. When the computation costs are low (proportion of computation costs of 10%, i.e., µc µ d = 0.1), the ratio evolution is consistent with the previous observations. Figure 5(a) is horizontally symmetrical as any case such that µc µ d > 1 (proportion of computation costs greater than 100%) is equivalent to the situation where the communication and the computation costs are swapped (and for which µc µ d < 1). These costs can be exchanged because a communication is always followed by a reduction operation.

Non-Commutative Operation. Finally, we assess the performance of Non-Commut-Tree-dyn by comparing it to all other methods that support a non-commutative operation when varying the dispersion and the overlapping degree as in the previous study.

Figure 5(b) shows the method with the best average performance when varying the CV and the ratio µc µ d . We see that Non-Commut-Tree-dyn has the best performance when the cost dispersion is large. Additionally, the transition from Binomial-stat to Fibonacci-stat is when the computation cost reaches half the communication cost (as on Figure 5(a)). With low computation costs (proportion of computation costs of 10%), the results are also consistent with Figure 4. Table II summarizes the tested distributions along with the corresponding parameters. While the CV varies from 0.01 to 10, the mean is always one. We discarded distributions with negative outcomes such as the normal distribution to have only positive costs. For some distributions, CVs greater than some threshold lead to non-positive support or to invalid parameters (the CV range is thus restrained for those distributions). Finally, we scale some distributions with two additional parameters to control the CV: m is an additive parameter (the minimum value in the case of the exponential distribution) and M is a multiplicative parameter (the maximum value in the case of the Bernoulli distribution).

On Figure 6(a), the same results shown by Figure 5(a) were aggregated for all ten distributions given by Table II. For a given scenario, that is, a given CV and a given µ c /µ d value, we first compute the ratio of the average time for Fibonacci-stat over the average time for Tree-dyn (over 1,000 MC simulations), for all possible distributions. Our objective is to check whether these ten ratios r are close to each other, or if some distributions lead to inconsistent results. To do this, we first compute the range of the ratios, and we normalize this range by the maximum distance between a ratio and one:

dispersion = max r -min r max |r -1| .
For instance, if all ratios are between 1.3 and 1.4, then the range is 0.1 and the maximum distance to 1 is 0.4, which gives a dispersion of 0.25. The rationale behind this normalization is to give a dispersion larger than 1 when results are inconsistent among all distributions, that is, if the Fibonacci-stat is faster than Tree-dyn for some distribution and slower for another distribution. Dark zones on Figure 5(a) correspond to scenarios where there are significant differences between the different distributions. Note that our method of combining a range with a relative measure purposely provides a pessimistic value (it considers the two most extreme distributions). On the contrary, other metrics such as the standard deviation would hide the differences between most when the CV is around 0.25 with 40% of computation (and its symmetric counterpart with 250% of computation); when the CV is around 0.8 with perfect overlapping; and when the CV is greater than 5. The first two areas are related to the proximity between both methods, hence a slight quantitative performance difference is amplified by the relative metric that is used (these transitions between both methods occur for the same set of parameters however). The last area is essentially due to the exponential method for which the improvement of the dynamic method is stronger than with other distributions. On overall, all distributions show significant quantitative consistency (95% of the dispersions are below 1 and the median dispersion is 0.27). Figure 6(b) aggregates results shown by Figure 5(b) for all distributions given by Table II. For a given scenario, we compute the expected best method for each distribution (using 1,000 MC simulations). The "inconsistency ratio" on Figure 6(b) is the proportion of the distributions for which the expected best method differs from the most frequent one. That is, an inconsistency ratio of 0.3 means that the best method for 3 out of the 10 distributions differs from the most frequent best method (which appears in the other 7 distributions). Once again, there is very little difference between the different distributions, and the discrepancies are restricted to the frontiers defining the zones where one strategy dominates the others.

In conclusion, using different distributions yields quantitative and qualitative consistent results, which suggests that our conclusions stand independently of the distribution.

CONCLUSION

In this paper, we have studied the problem of performing a non-clairvoyant reduction on a distributed heterogeneous platform. Specifically, we have compared the performance of traditional static algorithms, which build an optimized reduction tree beforehand, against dynamic algorithms, which organize the reduction at runtime. Our study includes both commutative and non-commutative reductions. We have first proposed approximation ratios for all commutative algorithms using a worst-case analysis. Then, we have evaluated all algorithms through extensive simulations using different random distributions to show when dynamic algorithms become more interesting than static ones. We have outlined that dynamic algorithms generally achieve better makespan, except when the heterogeneity is limited and for specific communication costs (no communication cost for Binomial-stat, communication costs equivalent to computation costs for Fibonacci-stat). The worst-case analysis has also confirmed this last observation.

The same symmetrical situation has been seen for c < d and c > d theoretically and empirically. The worst-case analysis of Binomial-stat and Fibonacci-stat is consistent with the simulation results: the later performs well when computations and communications overlap. However, the worst-case analysis has been unable to highlight the probabilistic advantage of Tree-dyn over Binomial-stat. Therefore, this theoretical study is ineffective for assessing the effect of heterogeneity between methods, contrarily to the empirical study.

As future work, we plan to investigate more complex communication models, such as specific network topologies. It would also be interesting to design a better dynamic algorithm for noncommutative reductions, which avoids the situation when many processors are idle but cannot initiate a communication since no neighboring processors are free.

Figure 2 .

 2 Figure 2. Schedules for Binomial-stat of order 3 and Fibonacci-stat of order 4, both using 8 processors. For Fibonacci-stat, the two schedules of order 2 and 3 used in the recursive construction are highlighted in green and red.

1 P 2 P 3 P 4 P 5 P 6 P 7 Figure 3 .

 12345673 Figure 3. Schedule of a worst-case scenario for Non-Commut-Tree-dyn.

1 log 2 n 2

 122 -approximation algorithm, where ϕ = 1+ √ 5is the golden ratio (1/ log 2 ϕ ≈ 1.44).

Figure 4 .

 4 Figure 4. Average schedule length for each method over 1,000,000 MC simulations with n = 64, µ d = 1, µc = 0, and varying coefficients of variation for the communication costs. The lower part of the ribbons corresponds to the 10% quantile while the upper part corresponds to the 90% quantile for each method.

Figure 5 .

 5 Figure 5. Ratio of the average performance of Fibonacci-stat and Tree-dyn (a) and method with the best average performance (b) over 1,000 MC simulations for each square with n = 64, µ d = 1, CV= σc µc = σ d µ d , varying the coefficient of variation (CV) for the costs and varying µc µ d .

Copyright c 2014

 2014 John Wiley & Sons, Ltd. Concurrency Computat.: Pract. Exper. (2014) Prepared using cpeauth.cls DOI: 10.1002/cpe

† http://www.sandia.gov/ ˜sjplimp/mapreduce.html Copyright c 2014 John Wiley & Sons, Ltd. Concurrency Computat.: Pract. Exper. (2014) Prepared using cpeauth.cls DOI: 10.1002/cpe

Copyright c 2014 John Wiley & Sons, Ltd. Concurrency Computat.: Pract. Exper. (2014) Prepared using cpeauth.cls DOI: 10.1002/cpe

log

ϕ -1 ≈ 0.44. This is the case when there is a significant overlap between computations and

Inconsistency ratio

(b) Non-commutative operation Figure 6. Dispersion of the average relative performance of Fibonacci-stat and Tree-dyn (a) and dispersion of the method with the best average performance (b) over 1,000 MC simulations for each square with n = 64, µ d = 1, CV= σc µc = σ d µ d , varying the coefficient of variation (CV) for the costs, varying µc µ d and with the ten distributions given by Table II. distributions, which are small. Note that most of this figure is almost white, which means that the results are coherent using all distributions. There are mostly three areas where there is divergence:

Table II. Experimented probability distributions with their parameters given a mean fixed to one and a variable CV. Distributions are scaled with the additional parameters m and M (the outcome is incremented by m or multiplied by M ). The CV is 0.5 for each example density.

Name Parameters CV range Density


K Deiters

A R Imre

C Ramboz

U K Deiters

T Kraska

Anomalous fluid properties of carbon dioxide in the supercritical region -Application to geological CO 2 storage and related hazards

Keywords: carbon sequestration, CCS, Widom-line

published or not. The documents may come

Anomalous fluid properties of carbon dioxide in the supercritical region -Application to geological CO 2 storage and related hazards

Introduction

Supercritical carbon dioxide is probably the most common supercritical fluid with a wealth of industrial applications. There are two reasons for that: CO 2 is readily available, and it is a low-cost fluid. Also, it exists as a supercritical fluid over a wide pressure and temperature range, starting just a few degrees above room temperature and above 7.3 MPa. A novel application for which the thermophysical properties of supercritical CO 2 need to be known is the "carbon capture and sequestration" (CCS). CCS has been identified as a realistic way to mitigate present global warming of the Earth. Figure 1 shows a sketch of the cycle of CO 2 capture/production, transportation, then injection stages involved in CCS. Major technological and knowledge gaps, however, have to be filled in order to ensure CCS safety, the main challenge being to prevent CO 2 -leakage over a timelength of ≈ 10 3 years.

Thermophysical properties of supercritical CO 2 mixtures, such as would be encountered in CCS cycle conditions, are identified as one field which particularly requires development of experimental and theoretical research (e.g., [START_REF] Bachu | Sequestration of CO2 in geological media: criteria and approach for site selection in response to climate change[END_REF][START_REF] Wilday | Safety issues for carbon capture and storage[END_REF].

It is known that some properties (e.g., compressibility or density) of supercritical fluids behave anomalously in a narrow temperature-pressure band that starts at the critical point and extends to higher T and P values [START_REF] Leneindre | Thermal conductivity of carbon dioxide and steam in the supercritical region[END_REF][START_REF] Wang | Effect of CO 2 flooding on wave velocities in rocks with hydrocarbons[END_REF][START_REF] Thiéry | A new object-oriented library for calculating analytically high-order multivariable derivatives and thermodynamic properties of fluids with equations of state[END_REF][START_REF] Estrada-Alexanders | Speed of sound in carbon dioxide at temperatures between (220 and 450) K and pressures up to 14 MPa[END_REF][START_REF] Xu | Relation between the Widom line and the dynamic crossover in systems with a liquid-liquid phase transition[END_REF][START_REF] Oschwald | Injection of fluids into supercritical environments[END_REF][START_REF] Mcmillan | Going supercritical[END_REF][START_REF] Simeoni | The Widom line as the crossover between liquid-like and gas-like behaviour in supercritical fluids[END_REF]. These anomalies can be represented in a P-T diagram by curves marking extrema or inflection points of the relevant thermodynamic properties; specific definitions are given in Section 2. These curves are collectively called Widom lines, although this term was originally used for the correlation length anomaly only [START_REF] Xu | Relation between the Widom line and the dynamic crossover in systems with a liquid-liquid phase transition[END_REF]. Different properties have different Widom lines, so this group of curves forms a Widom region, as has been shown for various real and model liquids (Brazhkin and Ryzhov 2011;Brazhkin et al. 2011;[START_REF] May | Riemannian geometry study of vapor-liquid phase equilibria and supercritical behavior of the Lennard-Jones fluid[END_REF]Imre et al. 2012[START_REF] Brazhkin | Universal crossover of liquid dynamics in supercritical region[END_REF][START_REF] Brazhkin | The "true" Widom line for a square-well system[END_REF]. In pressure-temperature space the Widom-lines are located in a wedge-shaped region (called Widom region), pointing to the critical point, marking the region where physico-chemical quantities (e.g., compressibility, heat capacities, density, thermal expansivity, speed of sound) show anomalous behaviour. Other names for this region are pseudocritical region, pseudoboiling region or pseudospinodal region. As was recently shown for water, the description of supercritical fluid flow involving this region requires special attention; some of the established methods for flow and thermal exchange calculations for normal liquids and gases cannot be used because they are not able to handle the abrupt changes related to these anomalies (Imre et al. 2012); for example with some well-established numerical codes used to describe fluid flow in nuclear power plants (like RELAP or WAHA), compressibility, heat capacity or density can be estimated only with a large uncertainty only (around 10% for the first two quantities and around 1% for density) in the vicinity of the Widom-line (Imre and Tiselj 2012). Therefore knowledge of the location of this anomalous region is crucial whenever flow (including leaking) and heat exchange involving supercritical states are concerned, so that process designers can be warned against P-T conditions that should be avoided orif this is not possibleat which advanced calculation methods are necessary. Anticipating similar problems for carbon dioxide in the case of CCS, an accurate mapping of the CO 2 Widom region is required in order to ensure that the supercritical CO 2 will remain controllable under the P-T conditions of deep geological storage (for the P-T range of existing storage places, see the case studies later).

In general, any rapid pressure and/or temperature change in the Widom region might induce an unexpectedly large response. For example, compressibility is the response function for a pressure change, expressed in terms of volume change. When a region of anomalously high compressibility is crossed, even a minute pressure and/or temperature change can lead to a hardly controllable volume change. Depending on the actual P-T path of the leak, either pseudo-boiling or pseudo-condensation (significant density change for a small T or P change, without real phase transition) can prevail. This can have a great influence on the leakage rate, since fluid flow is very sensitive to parameters changing during these pseudo-transitions (like density, compressibility, viscosity, etc.). Crossing the compressibility Widom line during the transportation or injection may cause problems with the flow regulation, because the very rapid compressibility variation can lead to pressure oscillations. Other quantities such as the viscosity, the correlation length and the diffusivity exhibit anomalous behaviour in the Widom region, too. Furthermore, the solubility of CO 2 is strongly density-dependent and changes significantly in the Widom region. This fact is particularly relevant in CCS, since CO 2 injection into deep saline aquifers is one of the potential ways envisaged to mitigate climate changes. Despite the fact that fluid density fluctuations in the Widom region obviously control the fluid capacity to dissolve CO2, we do not discuss their Widom lines in this work. This is because these properties cannot be readily calculated directly from equations of state without some additional assumptions.

Figure 1 shows a sketch of the cycle of CO 2 capture/production, transportation, then injection, including a potential leakage stage. From site A (production or capture site), the gas is transported to the injection site (B). This transport happens under fully controlled conditions, and pressure and temperature can be kept outside the Widom region (during transportation, T is usually below the critical temperature T c ). The injection process is also controlled (from injection site B to geological storage site C) although the well-bottom temperature is fixed by geological constraints. One or more Widom lines might be crossed during injection, but the effect of these anomaly-crossings can be at least partially controlled (orin case of a serious problemthe injection process can be stopped). The third step is the migration of carbon dioxide within the reservoir (from C to D). This is a fully uncontrolled process, but as it takes place inside the storage formation, the potential P-T paths cannot be very long, and therefore it is not likely that any anomaly lines are crossed. The most problematic step is the potential leakage process (D to E) from the reservoir to the surface (or to the sea-floor, in case of off-shore storage). Here the pressure-temperature route of the leakage is fully uncontrolled (although the initial and final conditions are fixed). Therefore, during leakage, crossing this anomalous region might cause some deleterious effect. Knowledge of Widom region-related property changes is necessary for evaluating the safety of anthropogenic CO 2 sequestration. Concerning the fate of CO 2 injected into deep saline aquifers, models initially predicted that, at such sites, the injected CO 2 plume would entirely be dissolved by the brine after ~4000-5000 years (e.g., Torp andGale 1996, Al-Fattah et al. 2011). At the Sleipner pilot site (North Sea) however, CO 2 has been continuously injected since 1996 into shallow, brine-saturated sandstone capped with a shale layer. In 2010 a mapping of the CO 2 plume revealed, that the speed of the dissolution process had been overestimated [START_REF] Ringrose | Sleipner and Snøhvit projects[END_REF]). Thus, the potential problems related to Widom anomalies are worthwhile evaluating.

In this short note we are going to map the anomalous region for pure CO 2 with respect to some thermodynamic properties relevant for flow and heat transfer, namely the density, the isothermal compressibility, the isobaric and the isochoric heat capacities, the isobaric thermal expansivity, and the speed of sound. Furthermore, the effects of some "contaminants" (like water or lower alkanes) will be considered. Finally, the potential problems related to the anomalies in the Widom region will be evaluated for six storage locations: Snøhvit and Sleipner in Norway, Nagaoka in Japan, Ketzin in Germany, and two natural CO 2 -storage analogues, Montmiral in France and Mihályi-Répcelak in Hungary.

Supercritical anomaliestheoretical description

Figure 2 shows the fluid phase diagram of carbon dioxide in P-T projection, calculated with the Wagner-Span carbon dioxide reference equation of state [START_REF] Span | A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1000 K at pressures up to 800 MPa[END_REF], using the ThermoC program [START_REF] Deiters | ThermoC[END_REF]. This equation of state (EoS) is highly accurate even for supercritical states; however, it cannot be used for mixtures. The vapour pressure curve runs from the triple point at about 216 K to the critical point at 304.2 K and 7.38 MPa [START_REF] Suehiro | Critical parameters of {x CO 2 + (1x) CHF 3 } for x = (1.0000, 0.7496, 0.5013, and 0.2522[END_REF]. The border of the supercritical region, conventionally defined by the two conditions P > P c and T > T c , is represented in Figure 2 by the two thick grey lines. An alternative definition of the supercritical region is also depicted in Figure 2.

Here the Widom region is represented by a narrow grey band, with vapour-like and liquid-like regions lying below and above that region, respectively. The latter regions are more or less indistinguishable from the normal pressurized liquid or vapour in the sense that there is no drastic change of their relevant thermodynamic properties (see for example Pártay el al. 2006;Imre et al. 2012).

In the figure also the so-called liquid-vapour stability limit (spinodal) is shown, which represents the thermodynamic limit of overheating. By carefull heating a liquid one can cross the vapour pressure curve without inducing a phase transition (boiling) [START_REF] Debenedetti | Metastable liquids: concepts and principles[END_REF]. In that case, , and as long as the spinodal is not approached, the CO 2 will remain liquid; but this state is metastable, and intenseeven explosiveboiling can happen at any moment. In a similar way, there is a limit of undercooling of the vapour phase, namely the vapour-liquid spinodal. The pressure of metastable liquids can be negative; in this case the metastability is induced by isotropic stretching rather than heating [START_REF] Skripov | Metastable Liquids[END_REF]Trevena 1986;[START_REF] Debenedetti | Metastable liquids: concepts and principles[END_REF]Imre et al. 2002). In this work, the stability lines have been derived from the Wagner-Span EoS [START_REF] Span | A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1000 K at pressures up to 800 MPa[END_REF]. This equation of state was developed, and its parameters optimized, for the stable region of the phase diagram. Therefore the extension of metastable region might not be reliable, as it was already stated in the original Span and Wagner paper (1996) and discussed later in a similar case by [START_REF] Imre | Estimation of the thermodynamic limit of overheating for bulk water from interfacial properties[END_REF]. This statement is valid for other reference equation of state, tooalthough this is often overlooked [START_REF] Imre | Estimation of the thermodynamic limit of overheating for bulk water from interfacial properties[END_REF].

Therefore the real locations of the spinodals might differ from the ones shown in Figure 2.

Metastable states located between the vapour pressure curve and the spinodals [START_REF] Kraska | The relation of interface properties and bulk phase stability: MD simulations of carbon dioxide[END_REF]) are responsible for some hysteresis which can also cause problems during CO 2 processing [START_REF] Benintendi | Non-equilibrium phenomena in carbon dioxide expansion[END_REF]. Although spinodals end at the critical point without extending into the supercritical region, they need to be mentioned here because the Widom lines can be regarded as their "shadows" cast beyond the critical point.

In this work the following Widom lines are considered for supercritical carbon dioxide:

• the compressibility Widom line: The isothermal compressibility  T is infinite along the spinodal curves and at the critical point. The associated Widom line is the locus of  T maxima on the P-T plane. Figure 3a shows such a  T maximum at the slightly supercritical pressure of 7.5 MPa. Figure 3d shows the location of the maxima on the three-dimensional P-T- T space. They evidently become sharper and higher as the critical point is approached.

•

the Widom lines of the isochoric or isobaric heat capacities C P and C V , and of the thermal expansivity  P ,: These properties become infinite at the critical point, and the associated Widom lines are therefore loci of maxima of these functions. Strong heat capacity means that a large amount of heat can be absorbed.

• the Widom line of the speed of sound: The speed of sound c s becomes zero at the critical point or along the spinodals; the associated Widom line is a locus of (nonzero) minima in the supercritical region (Figure 3b). This quantity has been included here as acoustic measurements are very important to map reservoirs.

• the density Widom line: The density has a discontinuity for subcritical temperatures. The associated Widom line is a locus of inflection points of the density function (Figure 3c). The density affects viscosity, diffusivity and solubility of CO 2 and its potential mixtures.

Figure 3 shows that, at constant pressure, the properties  T, c s and undergo drastic changes in a relatively narrow temperature range. Similarly, on isothermal plots, these properties show drastic changes in a narrow pressure range. These changes are not phase transitions, but bear some similarity to them. For example while in phase transition, the fluid exhibits and abrupt density change at a fixed temperature (like in boiling from liquid density to vapour density); in the supercritical region a similar density change might happen, although not abruptly at a temperature, but in a temperature range.

In Figure 4, the Widom-lines for the quantities  T, C P and C V , P , c s , and  are plotted on the P-T plane from the critical point to about 550 K and 50 MPa using the Wagner-Span equation of state [START_REF] Span | A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1000 K at pressures up to 800 MPa[END_REF]. The line corresponding to the isothermal compressibility is the two-dimensional P-T projection of the bold line in Figure 3d. The Widom line is plotted up to 7P c for carbon dioxide (where P c is the critical pressure). Beyond these pressures the extrema fall below the uncertainty of the equation of state and the effect of the Widom anomalies diminishes.

The Widom lines can be approximated by cubic polynomials:

T = A + BP + CP 2 + DP 3 Eq. (1)
The coefficients of the polynomials are given in Table 1. The polynomials include the critical point with an error of 10 mK.

It can be seen in Figure 4, that the anomalous zone (in pressure-temperature space) is bounded by the Widom line of the isothermal compressibility ( T ) from the high temperature side, while on the low temperature side it is bounded mostly by the Widom line of the density.

Processes involving pressure-temperature routes crossing or approaching this region can be influenced by the Widom line related anomalies.

The demonstration of a potential problem related to the crossing of density Widomline can be seen in Figure 5. An isothermal pressure drop (Figure 5a) at T = 306 K from 7.75 MPa to 7.65 MPa generates a density drop from 0.537 to 0.398 g/cm 3 , which means that a 0.1 MPa pressure drop (corresponding to a 10 m depth change in a hydrostatic pressure scale and an even smaller one on a lithostatic scale) can cause a density dropor a volume expansionof more than 25%. This massive volume increase is sometimes called pseudoboiling because of its similarity to the boiling of a liquid [START_REF] Oschwald | Injection of fluids into supercritical environments[END_REF]). In similar manner, an isobaric temperature drop at P = 7.6 MPa from 305.5 to 303.2 K (Figure 5b)

generates a density increase from 0.524 to 0.436 g/cm 3 , which means that a 0. The agreement with respect to the locations and depths of the minima is very satisfactory. The presence of the minimum can be seen even by using the original van der Waals EoS, but obviously, in that case the quantitative agreement is not so good.

Incidentally, it is remarkable that the existence of supercritical anomalies has been known from experiments for several decades, while theoretical studies of the Widom region began merely a decade ago [START_REF] Xu | Relation between the Widom line and the dynamic crossover in systems with a liquid-liquid phase transition[END_REF].

Impure carbon dioxide

The effect of a second component on the location and intensity of supercritical Widom anomalies is worth considering, because the stored carbon dioxide hardly remains pure.

Complex CO 2 -hydrocarbon fluid mixtures are known to be generated during CO 2 injection for enhanced oil recovery, and they are also found in some natural deep CO 2 reservoirs (e.g., Orr and Jensens 1984;NASCENT Report 2005, Novak et al. 2013[START_REF] Han | Non-parametric simulations-based conditional stochastic predictions of geologic heterogeneities and leakage potentials for hypothetical CO 2 sequestration sites[END_REF], Novak et al. 2014). Thus, even if CO 2 is injected as pure, it will mix with other compounds during storage, like nitrogen, oxygen, argon, hydrogen sulphide and water vapour, beside various gases like lower alkanes such as methane to hexane or heptane. Some organic contaminantslike benzene, toluene, xylene, etc.can even be dissolved from reservoir rocks and caprocks [START_REF] Zhong | Mobilization and transport of organic compounds from reservoir rock and caprock in geological carbon sequestration sites[END_REF]).

As a "rule of thumb" concerning the displacement of the Widom-region due to the addition of a small amount of a second component, one might expect that a contaminant with higher critical temperature would shift the Widom-peak to higher temperatures, while one with lower critical temperature would do the opposite. But here caution is advised: the point at which the compressibility diverges is not identical with the critical point of a mixture. A mixture may undergo a phase separation even if its compressibility is positive for all pressures. Therefore, it is possible that portions of Widom lines close to a mixture critical point lie in a two-phase region. Another difficulty is that some contaminants may even cause liquid-liquid phase separation. It is therefore advisable to compute the Widom lines of mixtures at some distance from their critical points only (Kraska andImre 2005, Deiters and[START_REF] Imre | The pseudocritical regions for supercritical water[END_REF].

Critical points of several potential contaminants are shown on Figure 7a (data taken from NIST Chemistry Webbook 2014). Figure 7b illustrates the shifts in the Widom compressibility peaks that are induced by the addition of various amounts of methane or water into CO 2 at a fixed pressure of 10 MPa. Since the Wagner-Span reference equation of state cannot be used for mixtures, here the Xiang-Deiters equation was used as implemented into the ThermoC package (Deiters 2006, Xiang and[START_REF] Xiang | A new generalized corresponding-states equation of state for the extension of the Lee-Kesler equation to fluids consisting of polar and larger nonpolar molecules[END_REF]. The compressibility peaks of methane-CO 2 with various methane contents are shifted to lower temperature compare to that of pure CO 2 , consistent with the fact that methane has lower critical temperature than CO 2 .

Similarly, the compressibility peaks of water-CO 2 mixtures are shifted toward higher temperature, consistent with the fact that the water critical temperature is higher than that of CO 2 .

There is no general rule to predict the displacement of the Widom region on adding a second (third, fourth, …) component and it is clear from Figure 7 that it can be shifted in any direction. Therefore, in CCS case studies, even when the P-T conditions of transport, storage or of a potential leakage path are far from the pure CO 2 Widom region, an unavoidable contamination of the stored carbon dioxide might change the situation in some extent.

Leakage hazards

Leakage is a well-known hazard for CCS. Leakage rates must be kept below 0.01%/year [START_REF] Hepple | Implications of surface seepage on the effectiveness of geological storage of carbon dioxide as a climate change mitigation strategy[END_REF]Benson, 2003, Miles et al. 2005) to ensure that CO 2 is stored over a sufficiently long time-span. The danger is complex; from small-scale leakage causing problems to local vegetation to potentially possible high-energy catastrophic pneumatic eruptions [START_REF] Pruess | On CO 2 fluid flow and heat transfer behavior in the subsurface, following leakage from a geologic storage reservoir[END_REF]. One has to consider that the number of cracks can increase drastically after the introduction of high pressure gas into the reservoir [START_REF] Eshiet | Investigation of geomechanical responses of reservoirs induced by carbon dioxide storage[END_REF]. Here we would like to discuss some leakage problems related to the Widom lines.

In Figure 8, various hydrostatic P-T lines can be seen, calculated with temperature gradients of 25 and 30 K/km, with average surface temperatures of 5 °C and 15 °C at atmospheric pressure of 0.1 MPa using averaged water density [START_REF] Pruess | On CO 2 fluid flow and heat transfer behavior in the subsurface, following leakage from a geologic storage reservoir[END_REF]). One hypothetical off-shore case is also considered, located below a 300 meter deep sea with a seafloor temperature of 2 °C. For potential leakage routes, only the lower and upper P-T conditions are fixed (or rather confined in a range) to be in agreement with the reservoir and surface or sea-bottom P-T data. Given that lithostatic P-T paths plot above hydrostatic ones, one can see that, in the case of pure CO 2 storage, the possibility that Widom-lines should be crossed during leakage paths could only occur in places with a high surface temperature and having relatively high geothermal gradient. Obviously, the situation changes when the stored CO 2 is contaminated, e.g. inside of the storage formation, by compounds having a low T c such as like methane. In this case the Widom-region would be shifted to lower temperatures, compatible with potential leakage routes.

While data for potential leakage routes are either non-existing or confidential, the P-T conditions of some existing CO 2 storage places can be compared with those of the pure CO 2 Widom region (Figure 9). We are not going to give detailed case studies; the selected examples serve only to demonstrate that Widom lines are located within the P-T range relevant in CSS. The studied locations are four existing storage places (Figures 9a-d: Snøhvit and Sleipner offshore Norway, Nagaoka in Japan, Ketzin in Germany) and two natural CO 2 reservoirs as analogues for anthropogenic CO 2 storage (Figures 9e-f: Montmiral in France and Mihályi-Répcelak in Hungary). The P-T conditions are taken from various sources (see hereafter). In the figures, some P-T paths of potential leaks in the selected sites are indicated by arrows; whenever these paths are crossing a Widom line, one should be prepared to handle it. For example in case of leakage in a network of narrow cracks, the temperature decrease will densify the uprising carbon dioxide, while the pressure drop will decrease its density.

Depending on the actual geological conditions such as the pressure regime and temperature gradients, the volume of escaped gas might increase and therefore it might widen the originally narrow cracks. In other leakage routes, the volume might decrease and in the latter case, the leakage should remain of low intensity.

Figures 9a shows the P-T storage conditions in the Utsira formation at Sleipner (North Sea), ranging from 25-42 °C and 8-11 MPa [START_REF] Baklid | Sleipner Vest CO 2 disposal, CO 2 injection into a shallow underground aquifer[END_REF][START_REF] Bickle | Modelling carbon dioxide accumulation at Sleipner: Implications for underground carbon storage[END_REF]Arts et al. 2010;[START_REF] Ringrose | Sleipner and Snøhvit projects[END_REF]. It can be seen that various leakage routes might cross the relevant Widom lines (lines b and c), although with low probability, because leakage probably would start from the upper part of the formation, along a P-T path like line "a", thus avoiding the anomalous region. The reported P-T conditions are close to the Widom region, therefore even a slight inaccuracyor a small amount of contaminantscan change to situation.

Indeed, the seafloor depth changed by as much as 21 cm between 2002 to 2005 at Sleipner, due to deep sediment compaction [START_REF] Alnes | Monitoring gas production and CO2 injection at the Sleipner field using time-lapse gravimetry[END_REF]. Also, currents can sweep away thick piles of sediments at the seabottom. In conclusion, concerning the potential problems caused by Widom line crossing, Sleipner can be considered as a proper location for storing pure supercritical CO 2 . However, this conclusion may need to be refined because of the 2 mole-% CH 4 contaminant (Arts et al. 2010), which might shift the compressibility maximum by approx. 5 K (see Figure 7b).

At the Snøhvit gas field (Barents Sea), the CO 2 that is pumped up with hydrocarbon gas to the sea floor at ~320 m water depth is re-injected into a sandy reservoir lying ~2.6 km below sea level. The initial P-T storage conditions were around 29 MPa and 98-100 °C in 2008, but reservoir pressure has increased to ~39 MPa in 2011 after near-continuous CO 2 injection, on account of limited permeability [START_REF] Hansen | Monitoring CO2 injection into a fluvial brine-filled sandstone formation at the Snøhvit field, Barents Sea[END_REF][START_REF] Hansen | Snøhvit: The history of injecting and storing 1 Mt CO 2 in the fluvial Tubåen Fm[END_REF]. In Figure 9b it can be seen that the reported P-T conditions are just in the middle of the anomalous region; therefore injection and leaking routes should cross them. In reality, contaminants (natural gas compounds) might push the Widom lines to lower temperatures, and the location of the storage state would probably be shifted between the Widom lines of thermal expansivity and compressibility, but into a region where anomalies are already weaker because the state is farther from the critical point. Snøhvit would not be a first pick for pure CO 2 storage, but because of some other considerationsthat the re-injected CO 2 is locally produced, so that transportation is not a problem, and that it can even enhance oil recovery (EOR)it might be a good choice.

The small-scale field of Nagaoka (Japan) was recently closed. Here well-head and well-bottom P-T data are available (7-11 MPa, 32 °C;19 MPa, 48 °C, respectively;[START_REF] Kikuta | Field test of CO 2 injection in Nagaoka, Japan[END_REF], but details about the P-T storage conditions and their evolution are lacking. Figure 9c shows that injection conditions just before or in the well-head or in the well-head might cross several Widom lines because they are close to the problematic critical point region.

Such instabilities can however be managed, as injection occurs under fully controlled conditions. Leaking paths theoretically might also intersect most Widom lines, but only along a very specific P-T path. Therefore from this point of view, Nagaoka is a proper storage place.

The P-T conditions of well Ktzi 200 from the Ketzin German pilot site are represented by a grey band in Figure 9d [START_REF] Kempka | Numerical simulations of CO 2 arrival times and reservoir pressure coincide with observations from the Ketzin pilot site, Germany[END_REF]Kuhn, 2013, Loizzo et al. 2013). Injection is done along the vapour pressure line, through the critical point and then hardly above the Widom lines. A leak might lead to a crossing of Widom lines in this very anomalous region, or to a phase transition (boiling) on crossing the vapour pressure line, or even to explosive boiling of depressurized superheated CO 2 [START_REF] Kraska | The relation of interface properties and bulk phase stability: MD simulations of carbon dioxide[END_REF][START_REF] Benintendi | Non-equilibrium phenomena in carbon dioxide expansion[END_REF]. Concerning these points, the Ketzin pilot is far from optimal, although it might score better in other respects [START_REF] Wiese | Well-based hydraulic and geochemical monitoring of the above zone of the CO 2 reservoir at Ketzin, Germany[END_REF], in particular a detailed geological knowledge of the storage site. It should be mentioned here, that according to some recent measurements, leakageat least in the present stagewas not found in Ketzin [START_REF] Wiese | Well-based hydraulic and geochemical monitoring of the above zone of the CO 2 reservoir at Ketzin, Germany[END_REF].

In case of Nagaoka and Ketzin, the storage conditions might reach subcritical states (pressure and temperature below the critical values); in that case a phase transition might occur. For such conditions, caution is required for flow-related calculations (see for example Schultz et al. 2013).

The study of natural deep CO 2 reservoirs as analogues of long-term storage sites, concerning possible CO 2 leakage, is also important to assess the safety of CCS over a time length of several thousand years. In Figure 9e Reservoir conditions would thus fall in the center of the Widom region and also closer to CO 2 critical conditions. Such an evolution would thus favour potential leakage. The fact that the two natural accumulation places have existed for ages without any known problem indicates the nearly infinitesimal probability of problematic leaking routes for these reservoirs.

Summary and Perspectives

Predicting the flow of supercritical CO 2 in reservoirs is a major challenge in connection with the safety of anthropogenic CO 2 storage. For proper modelling of the sequestration steps -CO 2 production/capture, transport, and injection into a reservoiras well as for assessing the effects of potential leakage it is essential to have reliable and accurate thermodynamic data of CO 2 in the relevant pressure and temperature ranges. conditions of injection and storage, most problems related to the anomalous behaviour of supercritical fluids can be avoided or at least minimized. Obviously, the best solutions would be to find reservoirs with caprocks having very low permeabilities [START_REF] Edlmann | Experimental investigation into the sealing capability of naturally fractured shale caprocks to supercritical carbon dioxide flow[END_REF] or to develop methods where only a small fraction of the fluid would be in the supercritical phase, like the storage of absorbed CO 2 in coal mines (see for example Gonzalez-Nicieza et al. from Table 1) representing the estimated locations of the extrema/inflection points of the marked quantities. [START_REF] Xiang | A new generalized corresponding-states equation of state for the extension of the Lee-Kesler equation to fluids consisting of polar and larger nonpolar molecules[END_REF].

2014).

By adding methane to the carbon dioxide (mass ratios are marked in the figure), the peaks moved to lower temperatures, while by adding water, they moved to higher temperatures.

Parameters: mole fractions of CO 2 and the contaminant, ----: water, -•-•-: methane.

  2 K temperature drop (corresponding to a 10 m depth change at 20 K/km thermal gradient or a 8 m depth change at 25 K/km) can induce a density increase of approximately 17%. In analogy to the pseudo-boiling, this phenomenon should be called pseudo-condensation. During a leakage, density decrease caused by the pressure decrease and the density increase caused by the temperature decrease compete, therefore -depending on the actual P-T path of the leakeither pseudo-boiling or pseudo-condensation can prevail.3. Supercritical anomalies: confronting theory with experimentsIn Figure6, speed of sound values of supercritical CO 2 calculated with the Wagner-Span EoS are compared with experimental results at 340.2 K[START_REF] Wang | Effect of CO 2 flooding on wave velocities in rocks with hydrocarbons[END_REF]), 308.2, 313.2, and 318.2 K (Han et al. 2010), 300 and 350 K (Estrada-Alexanders and Trusler 1998).

  and 9f the P-T conditions of two such natural reservoirs (natural high-purity carbon dioxide accumulation) can be seen, Montmiral in France (Figure9e) and Mihályi-Répcelak in Hungary (Figure9f)(NASCENT Report 2005;[START_REF] Cs | Expected physical and chemical effects of injecting industrial carbon dioxide on pore fluids and reservoir rocks based on the study of Mihályi-Répcelak natural CO 2 occurrence[END_REF]. The CO 2 at the former site has been exploited and it is still exploited in the latter. At Montmiral, CO 2 is trapped in a 2400 m deep reservoir at 100 °C and 36 MPa. The CO 2 is associated with oil, and the exploited CO 2 gas contains ~1.6 mole-% CH 4(Pauwels et al. 2007). One of the production wells at Mihályi-Répcelak is 1.2 km deep, with corresponding CO 2 storage conditions of 73-75 °C and 14.7-15.1 MPa. The CO 2 is quite pure, with only a few mole-% water. Because reservoir pressure is high in both places (2 to 4 times the CO 2 critical pressure), the Widom lines would be crossed far above the critical point region in case of leakage, where the anomalies are less intense. For most of the specific P-T leakage paths considered (like the middle arrow in Mihályi-Répcelak), the anomalies in the crosscut Widom region should not be sharp, but gradual, without any abrupt change. At Montmiral, the concentration in CH 4 is too low to shift the Widom region to lower temperatures, thus potentially placing the reservoir in the middle of the Widom-region (Figure9/e). Finally, because both reservoirs are geopressurized (pressure above hydrostatic), abrupt decompression paths can be envisaged. At Montmiral particularly, a change from geopressurized to hydrostatic regime would correspond to a drop in pressure of ~10 MPa.

For

  supercritical carbon dioxidejust like for other supercritical fluidsthere is a wedge-shaped region called Widom region, pointing to the critical point, where several physico-chemical quantities (e.g., compressibility, heat capacities, density, thermal expansivity, speed of sound) show anomalous behaviour. These properties become extremely sensitive to pressure or temperature variations and can therefore change drastically along supercritical P-T paths. In geological environments the pressure and temperature depends on the actual depts. and the aforementioned properties are also depth-dependent. For example density might change even 25 % by changing the depth with 10 meter. This theoretical evaluation is quite in line with the experimental conclusion following ~15 years of deep CO 2 injection, that buoyancy forces are largely underestimated by current models[START_REF] Ringrose | IEA-GHG Summer School[END_REF].In this paper several Widom lines of supercritical CO 2 have been computed with the Wagner-Span reference EoS. The locations of these lines are compared with the P-T range of four storage sites (Snøhvit and Sleipner offshore Norway, Nagaoka in Japan, Ketzin in Germany, and two natural CO 2 storage analogues, Montmiral in France and Mihályi-Répcelak in Hungary), and the potential consequences of crossing any of the Widom-lines because of leaking CO 2 were discussed. The effect of contamination with natural gas compounds on the location of the Widom region is also discussed. It is shown that, by properly choosing the P-T

Figure 1 .

 1 Figure 1. Schematic representation of carbon dioxide capture (production), transportation, injection, storage, migration and leakage (see text).

Figure 2 .

 2 Figure 2. Phase diagram of fluid carbon dioxide (without solid phases) showing the limits of the conventionally defined supercritical region (P > P c and T > T c domain, bracketed between the two thick grey lines) and the novel description showing the liquid-like, vapour-like and transitional (Widom) region. For subcritical region, various stable and metastable phases marked, separated by the (solid) vapour pressure line, the (dashed) liquid-vapour stability line (also called liquid-vapour spinodal) and by the (dot-dashed) vapour-liquid stability line (vapour-liquid spinodal).

Figure 3 .

 3 Figure 3. (a) Isothermal compressibility  T (b), speed of sound c s and (c) density anomalies  in the vicinity of the critical point, calculated by the Wagner-Span reference equation of state (P = 7.5 MPa; the dashed line indicates the critical temperature). (d) Three-dimensional pressure-temperature-compressibility plot, showing the corresponding Widom line (bold) running on the ridge, connecting the maxima.

Figure 4 .

 4 Figure 4. Some relevant Widom lines of carbon dioxide. Points and error bars: calculated with the Wagner-Span reference equation of state, curves: polynomials (Eq. 1 with parameters

Figure 5 .

 5 Figure5. a: The phenomenon of "pseudo-boiling": an isothermal pressure drop (T = 306 K) of 1 bar generating a density drop of more than 25%. b: The phenomenon of "pseudocondensation": an isobaric temperature drop (P = 7.6 MPa) of 0.2 K causing a 17% density increase.

Figure 6 .

 6 Figure 6.: Comparison of measured (symbols) and calculated (lines) speed of sound values c s as a demonstration of the applicability of the Wagner-Span equation of state to determine the location and extent of Widom anomalies. Sources of the data: 340.2 K[START_REF] Wang | Effect of CO 2 flooding on wave velocities in rocks with hydrocarbons[END_REF]), 308.2, 313.2 and 318.2 K (Han et al. 2010), 300 K and 350 K (Estrada-Alexanders and Trusler, 1998).

Figure 7 .

 7 Figure 7. (a) Critical points of various relevant contaminants (NIST Chemistry Webbook 2011). (b) The effect of one potential "contaminant" in supercritical CO 2 on the location and

Figure 8 .

 8 Figure 8. Comparison of the location of Widom lines of CO 2 with various hydrostatic P-T routes. x: sea-floor conditions, grey circle: critical point.

Figure 9 .

 9 Figure 9. Temperature and pressure relations of CO 2 storage places and natural CO 2 accumulations compared to the location of Widom lines. Potential leakage routes are marked by arrows.

  Figure 2

  Figure 4.

Figure 5 .

 5 Figure 5.

  Figure 6.

  Figure 7.

  Figure 8.

Figure

  Figure 9.

Table 1 :

 1 Polynomial coefficients of Eq. 1 for the Widom lines studied in this work. R 2 : correlation coefficient, range: pressure range of validity. The uncertainties of the coefficients are indicated by the last decimal digits.

	Quantity	A/K	B/K•MPa -1	C/K•MPa -2	D/K•MPa -3 R 2	Range/MPa
	c s	244.5	8.37	-0.062	0.00031	0.99992	7.5-50
	C pm	265.3	6.18	-0.097	0.00004	0.99251	7.5-50
	C Vm	233.2	12.73	-0.464	0.00482	0.99704	7.5-20
		260.6	7.22	-0.169	0.00111	0.99700	7.5-50
	 p	255.7	7.65	-0.149	0.00113	0.99989	7.5-50
	 T	249.2	8.57	-0.221	0.00796	0.99989	7.5-30

Acknowledgments

A. R. I. gratefully acknowledges an "Albert's reunion grant" from the University of Cologne and an invitation as a guest professor at the University of Orléans for two months. This work could not have been completed without the support of these two Universities.