Pegah Kamousi

Sylvain Lazard

Anil Maheshwari

Stefanie Wuhrer

Analysis of Farthest Point Sampling for Approximating Geodesics in a Graph

Keywords: Farthest Point Sampling, Approximate Geodesics, Shortest Paths, Planar Graphs, Approximation Algorithms

A standard way to approximate the distance between any two vertices p and q on a mesh is to compute, in the associated graph, a shortest path from p to q that goes through one of k sources, which are well-chosen vertices. Precomputing the distance between each of the k sources to all vertices of the graph yields an efficient computation of approximate distances between any two vertices. One standard method for choosing k sources, which has been used extensively and successfully for isometryinvariant surface processing, is the so-called Farthest Point Sampling (FPS), which starts with a random vertex as the first source, and iteratively selects the farthest vertex from the already selected sources.

In this paper, we analyze the stretch factor F F P S of approximate geodesics computed using FPS, which is the maximum, over all pairs of distinct vertices, of their approximated distance over their geodesic distance in the graph. We show that F F P S can be bounded in terms of the minimal value F * of the stretch factor obtained using an optimal placement of k sources as F F P S 2r 2 e F * + 2r 2 e + 8r e + 1, where r e is the ratio of the lengths of the longest and the shortest edges of the graph. This provides some evidence explaining why farthest point sampling has been used successfully for isometry-invariant shape processing. Furthermore, we show that it is NP-complete to find k sources that minimize the stretch factor.

Introduction

In this work, we analyze the stretch factor of approximate geodesics computed on triangle meshes or more generally in graphs. In our context, a triangle mesh represents a a Département d'Informatique, Université Libre de Bruxelles, Belgium, E-mail: pegahk@gmail.com b Inria, Loria, Nancy, France, E-mail: sylvain.lazard@inria.fr c School of Computer Science, Carleton University, Canada, E-mail: anil@scs.carleton.ca d Cluster of Excellence MMCI, Saarland University, Germany, E-mail: swuhrer@mmci.uni-saarland.de

G G ′ G ′′ Figure 1:
The graph G = (V, E, F ) is a triangle mesh, while G ′ has a non-triangular face, and G ′′ is not connected.

discretization of a two-dimensional manifold, possibly with boundary, embedded in R d for a constant dimension d using a set of vertices V , edges E, and triangular faces F . Such a triangle mesh can be viewed as a (hyper)graph G = (V, E, F ) such that each edge is adjacent to one or two triangles and the triangles incident to an arbitrary vertex can be ordered cyclically around that vertex. Due to discretization artifacts, the triangle mesh may contain holes, and different parts of the surface may even intersect in the embedding. However, we assume that the graph structure G is planar and connected (see Figure 1 for an illustration). For d = 3, this definition of a triangle mesh is commonly used in geometry processing when analyzing models obtained by scanning real-world objects.

Given a connected planar triangle graph G with n vertices, where every edge has a positive length (or weight), we consider the problem of approximating the geodesic distances between pairs of vertices in G, where distances are measured in the graph theoretic sense. Specially, given an integer k, our goal is to select a set S = {s 1 , . . . , s k } of k vertices of G that minimize the stretch factor, defined as the value max (p,q)∈V, p =q min s i ∈S d(p, s i ) + d(s i , q) d(p, q) , where the function d(., .) measures the shortest geodesic distance between two vertices. Throughout this paper, we use for simplicity the notation max p,q for max (p,q)∈V, p =q .

The problem of approximating geodesic distances on surfaces represented by triangle meshes arises when studying shapes (usually shapes embedded in R 2 or in R 3 ) that are isometric, which means that they can be mapped to each other in a way that preserves geodesics. Isometric shapes have been studied extensively recently [START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF][START_REF] Elad | On bending invariant signatures for surfaces[END_REF][START_REF] Lipman | Möbius voting for surface correspondence[END_REF][START_REF] Sun | A concise and provably informative multiscale signature based on heat diffusion[END_REF][START_REF] Tevs | Intrinsic shape matching by planned landmark sampling[END_REF] because many shapes, such as human bodies, animals, and cloth, deform in a nearisometric way as a large stretching of the surfaces would cause injury or tearing.

In order to analyze near-isometric shapes, it is commonly required to compute, on each shape, geodesics between many pairs of vertices, and to compare the corresponding geodesics in order to compute the amount of stretching of the surface. Consider the problem of computing geodesic distances on a surface represented by a triangle mesh, where distances are measured on the graph induced by the vertices and edges of the triangulation. To solve the single-source shortest-path (SSSP) problem on G, Dijkstra's algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF] takes O(n log n) time. 1 To solve the all-pairs shortest-path (APSP) problem, we can run Dijkstra's algorithm starting from each source point, yielding an O(n 2 log n) time algorithm. While there are more efficient methods, the APSP problem has a trivial Ω(n 2 ) lower bound. In practical applications, a typical triangle mesh may contain from 100,000 to 500,000 vertices and, for such large meshes, it is impractical to use algorithms that take Ω(n 2 ) time.

To allow for a reduced complexity, instead of considering the APSP problem, we consider the problem of pre-computing a data structure that allows to efficiently approximate the distance between any two points. We call this the any-pair approximate shortest path problem in the following.

A commonly used method to solve the any-pair approximate shortest path problem is to select a set of k sources, solve the SSSP problem from each of these k sources, and use this information to approximate pairwise geodesic distances. Given a pair p and q of vertices, their shortest distance is approximated as the minimum, over all the k sources, of the sum of the distances of p and q to a source s i . This method is used to approximate the intrinsic geometry of shapes [START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF][START_REF] Elad | On bending invariant signatures for surfaces[END_REF][START_REF] Mémoli | Comparing point clouds[END_REF]. A natural problem is thus to compute an optimal placement of k sources that minimizes the stretch factor. We refer to this problem as the k-center path-dilation problem and show that this problem is NP-complete (see Theorem 6).

A commonly used heuristic for selecting a set of k sources is to use Farthest Point Sampling (FPS) [START_REF] Gonzales | Clustering to minimize the maximum intercluster distance[END_REF][START_REF] Moenning | Fast marching farthest point sampling[END_REF], which starts from a random vertex and iteratively adds to S a vertex that has the largest geodesic distance to its closest already picked source, until k sources are picked. Given k sources on a graph G, the distance between any two vertices p and q is approximated as the minimum over all k sources, of the distance from p to q through one of the sources.

FPS has been shown to perform well compared to other heuristics for isometry-invariant shape processing in practice [START_REF] Ruggeri | Spectral-driven isometryinvariant matching of 3d shapes[END_REF] [19, Chapter 3], which suggests that the stretch factor obtained by a FPS is small. However, to the best of our knowledge, no theoretical results are known on the quality of the stretch factor, F F P S , obtained by a FPS of k sources, compared to the minimal stretch factor, F * , obtained by an optimal choice of k sources.

In this paper, we prove that

F F P S 2r 2 e (F * + 1) + 8r e + 1,
where r e is the ratio of the lengths of the longest and the shortest edges of G (see Theorem 1). Note that this bound holds for any arbitrary graph.

It should further be observed that if the ratio r e is large, F F P S can be much larger than the optimal stretch factor F * but, on the other hand, F * is likely to be large as well. Indeed, if at least k + 1 edges are arbitrarily small and are not "too close" to each other, F * can be made arbitrarily large; this can be seen by considering the pairs of vertices defined by those small edges.

After discussing related work in Section 2, we prove our two main results, Theorems 1 and 6 in Sections 3 and 4, respectively.

Related Work

Computing geodesics on polyhedral surfaces is a well-studied problem for which we refer to the recent survey by Bose et al. [START_REF] Bose | A survey of geodesic paths on 3d surfaces[END_REF]. In this paper, we restrict geodesics to be shortest paths along edges of the underlying graph.

The FPS algorithm has been used for a variety of isometry-invariant surface processing tasks. The algorithm was first introduced for graph clustering [START_REF] Gonzales | Clustering to minimize the maximum intercluster distance[END_REF], and later independently developed for 2D images [START_REF] Eldar | The farthest point strategy for progressive image sampling[END_REF] and extended to 3D meshes [START_REF] Moenning | Fast marching farthest point sampling[END_REF]. Ben Azouz et al. [START_REF] Ben Azouz | Approximations of geodesic distances for incomplete triangular manifolds[END_REF] and Giard and Macq [START_REF] Giard | From mesh parameterization to geodesic distance estimation[END_REF] used this sampling strategy to efficiently compute approximate geodesic distances, Elad and Kimmel [START_REF] Elad | On bending invariant signatures for surfaces[END_REF] and Mémoli and Sapiro [START_REF] Mémoli | Comparing point clouds[END_REF] used FPS in the context of shape recognition. Bronstein et al. [START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF] and Wuhrer et al. [START_REF] Wuhrer | Posture invariant correspondence of incomplete triangular manifolds[END_REF] used FPS to efficiently compute point-to-point correspondences between surfaces. While it has been shown experimentally that FPS is a good heuristic for isometry-invariant surface processing tasks [START_REF] Ben Azouz | Approximations of geodesic distances for incomplete triangular manifolds[END_REF][START_REF] Giard | From mesh parameterization to geodesic distance estimation[END_REF][START_REF] Elad | On bending invariant signatures for surfaces[END_REF][START_REF] Mémoli | Comparing point clouds[END_REF][START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF][START_REF] Wuhrer | Posture invariant correspondence of incomplete triangular manifolds[END_REF], to the best of our knowledge, the worst-case stretch of the geodesics has not been analyzed theoretically.

The problem we study is closely related to the k-center problem, which aims at finding k centers (or sources) s ′ i , such that the maximum distance of any point to its closest center is minimized. With the notation defined above, the k-center problem aims at finding s ′ i , such that max p (min i d(p, s ′ i )) is minimized. This problem is N P -hard and FPS gives a 2-approximation, which means that the k centers s i found using FPS have the property that max p (min i d(p, s i )) 2 max p (min i d(p, s ′ i )) [START_REF] Gonzales | Clustering to minimize the maximum intercluster distance[END_REF].

In the context of isometry-invariant shape processing, we are interested in bounding the stretch induced by the approximation rather than ensuring that every point has a close-by source. A related problem that has been studied in the context of networks by Könemann et al. [START_REF] Könemann | An approximation algorithm for the edge-dilation k-center problem[END_REF] is the edge-dilation k-center problem, where every point, p, is assigned a source, s p , and the distance between two points p and q is approximated by the length of the path through p, s p , s q , and q. The aim is then to find a set of sources that minimizes the worst stretch, and Könemann et al. show that this problem is N P -hard and propose an approximation algorithm to solve the problem.

Könemann et al. [START_REF] Könemann | An approximation algorithm for the edge-dilation k-center problem[END_REF] also study a modified version of the above problem, which is similar to our problem. In particular, they present an algorithm for computing k sources and claim that it ensures, for our problem, a stretch factor of F K 2F * + 1 [11, Theorem 3]2 ; as before, F * denotes the minimal stretch factor for the k-center path-dilation problem.

However, we believe that their proof has gaps. 3 4 It should nonetheless be stressed that our result is independent of whether this bound on F K holds. Indeed, the relevance of our bound of 2r 2 e (F * + 1) + 8r e + 1 on F F P S is to give some theoretical insight on why FPS has been used successfully in heuristics for isometry-invariant shape processing.

Approximating Geodesics with Farthest Point Sampling

We start this section with some definitions and notation. We consider a connected graph G in which the edges have lengths from a positive and finite interval [ℓ min , ℓ max ], and r e denotes the ratio ℓ max /ℓ min . We require the graph to be connected so that the distance between any two vertices is finite. In this section, we do not require the graph to satisfy any other criteria, but observe that if it is not planar, the running time of FPS will be O(k(m + n log n)), where m is the number of edges.

Given k vertices (sources) s 1 , . . . , s k in the graph, let s p denote the (or a) closest source to a vertex p and let d(p, s i , q) denote the shortest path length from p to q through any source s 1 , . . . , s k , that is min

i (d(p, s i )+d(s i , q)). Let s * 1 , • • • , s * k be a choice of sources that minimizes the stretch factor F * = max p,q d(p, s * i , q)/d(p, q). Furthermore, let s ′ 1 , • • • , s ′ k
be a choice of sources that minimizes max p d(p, s p ). In other words, the set of s * i is an optimal solution to k-center path-dilation problem and the set of s ′ i is an optimal solution to the k-center problem.

In this section, we prove the following theorem.

Theorem 1. Let s 1 , • • • , s k be a set of sources returned by the FPS algorithm on a connected graph G with edge lengths of ratio at most r e . Then max p,q d(p, s i , q) d(p, q) 2r 2 e max p,q d(p, s * i , q) d(p, q) + 2r 2 e + 8r e + 1.

In order to prove this theorem, we first show a somewhat surprising property that, for any set of sources, the stretch factor max p,q d(p,s i ,q) d(p,q) is realized when p and q are adjacent in the graph (Lemma 2). We use this property to bound this stretch factor in terms of max p d(p, s p ) (Lemma 3). On the other hand, we bound the stretch factor of any set of sites in terms of the stretch factor of an optimal set of sources for the k-center problem (Lemma 5). We then combine these results to prove Theorem 1.

Lemma 2. For any sources s 1 , . . . , s k and any given vertex q in G, the maximum ratio max p d(p,s i ,q) d(p,q) is realized for some p that is adjacent to q in G. It follows that the maximum ratio max p,q d(p,s i ,q) d(p,q) is realized for some p and q that are adjacent in G.

Proof. For the sake of contradiction, let q be any fixed vertex and let p be a non-adjacent vertex that realizes the maximum max p d(p,s i ,q) d(p,q) and such that among all the vertices p ′ that realize this maximum, the shortest path from p to q has the smallest number of edges.

Let p be the immediate neighbor of p along the shortest path from p to q. As before, d(p, s j , q) denotes the shortest path length from p to q through any source s 1 , . . . , s k (we use here the notation d(p, s j , q) instead of d(p, s i , q) in order to avoid confusion with d(p, s i , q)). Let ℓ be the length of the edge pp (see Figure 2). We have d(p, s j , q) d(p, s i , q) -ℓ. Dividing by d(p, q) = d(p, q) -ℓ we get

d(p, s j , q) d(p, q) d(p, s i , q) d(p, q) -ℓ - ℓ d(p, q) .
On the other hand, by multiplying d(p, q) = d(p, q) + ℓ by d(p,s i ,q) d(p,q)d(p,q) we have

d(p, s i , q) d(p, q) -ℓ = d(p, s i , q) d(p, q) + ℓd(p, s i , q) d(p, q) • d(p, q) ,
and therefore d(p, s j , q) d(p, q) d(p, s i , q) d(p, q) + ℓ d(p, q)

• d(p, s i , q) d(p, q) -1 d(p, s i , q) d(p, q) , which contradicts our assumption. Indeed, either the inequality is strict and d(p,s i ,q) d(p,q) was not maximum, or the equality holds and the shortest path from p to q has fewer edges than the shortest path from p to q. p q p s i The property of the previous lemma that max p,q d(p,s i ,q) d(p,q) is realized when p and q are neighbors allows us to bound it as follows.

Lemma 3. For any sources s 1 , . . . , s k , we have

2 ℓ max max p d(p, s p ) -1 max p,q d(p, s i , q) d(p, q) 2 ℓ min max p d(p, s p ) + 1.
Proof. For the upper bound, we have d(p, s i , q) d(p, s p ) + d(s p , q) 2d(p, s p ) + d(p, q). Therefore, d(p,s i ,q) d(p,q) 2 d(p,q) d(p, s p ) + 1. This holds for any vertices p and q and thus for those that realize the maximum of d(p,s i ,q) d(p,q) . Furthermore, d(p, q) ℓ min and d(p, s p ) max p d(p, s p ). Hence, max p,q d(p, s i , q) d(p, q)

2 ℓ min max p d(p, s p ) + 1.
For the the lower bound, we have by the triangle inequality that, for any i, d(q, s i ) d(p, s i )-d(p, q). Adding d(p, s i ) on both sides, we get d(p, s i )+d(q, s i ) 2d(p, s i )-d(p, q). By the definition of s p , d(p, s i ) d(p, s p ) for any i, thus d(p, s i ) + d(q, s i ) 2d(p, s p )d(p, q). This holds for any i and thus for the i such that d(p, s i ) + d(q, s i ) is minimum, hence d(p, s i , q) 2d(p, s p ) -d(p, q). Dividing by d(p, q), we get d(p,s i ,q) d(p,q) 2 d(p,q) d(p, s p ) -1. This holds for any p and q and thus for the vertex p that realizes the maximum of d(p, s p ); let p denote such vertex. We then have that d(p,s i ,q) d(p,q) 2 d(p,q) max p d(p, s p ) -1. This holds for any q and in particular for the one that realizes max q d(p,s i ,q) d(p,q) . By Lemma 2, the maximum is realized for a q that is adjacent to p in G, thus, for such a q, 2 d(p,q) 2 ℓmax . It follows that max p,q d(p, s i , q) d(p, q) max q d(p, s i , q) d(p, q)

2 ℓ max max p d(p, s p ) -1.
The following lemma bounds the path length between two vertices u and v passing through s u in terms of the shortest path between u and v through any source.

Lemma 4. For any sources s 1 , . . . , s k , and vertices u, v we have

d(u, s u ) + d(s u , v) d(u, s i , v) + 2d(u, v).
Proof. Denote by s i the source that realizes the minimum

d(u, s i , v) = min i (d(u, s i ) + d(s i , v)).
Since by definition d(u, s u ) d(u, s i ), we only have to show that d(v, s u ) d(s i , v) + 2d(u, v). Using the triangle inequality twice, we have

d(v, s u ) d(v, u) + d(u, s u ) d(v, u) + d(u, s i ) d(v, u) + d(u, v) + d(v, s i ),
which concludes the proof.

These results allow us to bound the stretch factor corresponding to the sources returned by the FPS algorithm with respect to the stretch factor corresponding to an optimal choice of sources for the k-center problem.

d(p, s i , q) d(p, q) 2r e max u,v d(u, s ′ i , v) d(u, v)
+ 6r e + 1.

Proof. Since s 1 , . . . , s k is a set of sources returned by the FPS algorithm, this choice of sources provides a 2-approximation for the k-center problem compared to an optimal solution s ′ 1 , . . . , s ′ k ; in other words, max p d(p, s p ) 2 max p d(p, s ′ p ) [START_REF] Gonzales | Clustering to minimize the maximum intercluster distance[END_REF].

By definition, d(p, s i , q) is the minimum over all (fixed) sources s i of d(p, s i ) + d(s i , q). Thus, d(p, s i , q) d(p, s p ) + d(s p , q). Moreover, by the triangle inequality, d(s p , q) d(s p , p) + d(p, q) thus d(p, s i , q) 2d(p, s p ) + d(p, q). One the other hand, d(p, s p ) max u d(u, s u ), which is less than or equal to 2 max u d(u, s ′ u ) by the 2-approximation property. For clarity, denote by u the vertex that realizes the maximum max u d(u, s ′ u ). We then have d(p, s i , q) 4d(u, s ′ u ) + d(p, q). Now, by the triangle inequality,

d(u, s ′ u ) d(u, v) + d(v, s ′ u ) for any vertex v. Thus 2d(u, s ′ u ) d(u, v) + d(v, s ′ u ) + d(u, s ′ u ) which implies, by Lemma 4, that 2d(u, s ′ u ) 3d(u, v) + d(u, s ′ i , v). Thus, d(p, s i , q) 2d(u, s ′ i , v) + 6d(u, v) + d(p, q) and d(p, s i , q) d(p, q) 2 d(u, v) d(p, q) d(u, s ′ i , v) d(u, v) + 6 d(u, v) d(p, q) + 1.
This inequality holds for any distinct p and q, and any v distinct from u (recall that u is fixed). Thus it holds for the vertices p and q that realize max p,q d(p,s i ,q) d(p,q) and for the v that realizes max v d(u,s ′ i ,v) d(u,v) . Such a v is a neighbor of u by Lemma 2, thus it satisfies d(u, v) ℓ max . Since d(p, q) ℓ min for any distinct p and q, and max

v d(u,s ′ i ,v) d(u,v) max u,v d(u,s ′ i ,v) d(u,v) , we get max p,q d(p, s i , q) d(p, q) 2 ℓ max ℓ min max u,v d(u, s ′ i , v) d(u, v) + 6 ℓ max ℓ min + 1.
This finally allows us to prove the main theorem.

Proof of Theorem 1. By Lemma 5 and using the same notation, we have

max p,q d(p, s i , q) d(p, q) 2r e max u,v d(u, s ′ i , v) d(u, v) + 6r e + 1.
Using the upper bound in Lemma 3 on max u,v d(p, s * i , q) d(p, q) and thus max p,q d(p, s i , q) d(p, q) 2r e ℓ max ℓ min (max p,q d(p, s * i , q) d(p, q) + 1) + 1 + 6r e + 1

d(u,s ′ i ,v) d(u,v) ,
2r 2 e max p,q d(p, s * i , q) d(p, q) + 2r 2 e + 8r e + 1.

The Complexity of k-Center Path-Dilation Problem

In this section we consider the complexity of the k-center path-dilation problem on triangle graphs, i.e., computing an optimal set of sources that minimizes the stretch factor.

The following theorem shows that the decision version of this problem is NP-complete for triangle graphs. Note that this directly yields the NP-completeness for arbitrary graphs (since proving that the problem is in NP is trivial).

Theorem 6. Given a triangle graph G = (V, E, F ), an integer k, and a real value ξ, it is NP-complete to determine whether there exists a set S = {s 1 , . . . , s k } of k sources such that the stretch factor max p,q min i d(p,s i )+d(s i ,q) d(p,q)

is at most ξ.

Note that the problem is in NP since, for any set of k sources, the stretch factor can be computed in polynomial time. To show the hardness, we provide a reduction from the decision problem related to finding a minimum cardinality vertex cover on planar graphs of maximum vertex degree three [START_REF] Garey | The Rectilinear Steiner Tree Problem is NP-Complete[END_REF]. The first step of the reduction uses the following well-known result on embedding planar graphs in integer grids [START_REF] Valiant | Universality Considerations in VLSI Circuits[END_REF].

Lemma 7. A planar graph G = (V, E) with maximum degree 4 can be embedded in the plane using O(|V | 2 ) area in such a way that its nodes have integer coordinates and its edges are drawn as polygonal line segments that lie on the integer grid (i.e, every edge consists of one or more line segments that lie on lines of the form x = i or y = j, where i and j are integers).

Consider a planar graph G with maximum degree 3, and let G r be a planar embedding of G according to Lemma 7, to which we have added, on each edge e ∈ E, an even number of 2k e auxiliary nodes with half-integer coordinates and such that every resulting edge in G r has length 1 or 1/2. (We consider the half-integer grid so that we can ensure that we add an even number of auxiliary nodes on every edge of G so that every resulting edge in G r has length at most 1.) Please refer to Figure 3(a,b) for an illustration. For an edge uv ∈ E, we let path(u, v) denote the path in G r replacing the edge uv. The endpoints of the paths (i.e., the nodes that are not auxiliary), are called regular nodes. Finally, let m = e∈E k e . We have the following lemma. Now let C r be a vertex cover for G r of size m + k, and suppose there exists a path, path(u, v), such that neither u nor v belongs to C r . Then at least k uv + 1 auxiliary nodes from path(u, v) must belong to C r in order to cover all the edges of this path. However, by using only k uv auxiliary nodes from path(u, v) and adding u or v to C r , we still have a vertex cover of the same size, which now contains one of the endpoints of path(u, v).

Continuing this way, we can construct a vertex cover C r of size k + m for G r , which includes at least one endpoint from each path(u, v), for all uv ∈ E. Therefore, C r is a vertex cover for G, when restricted to the nodes of G (regular nodes). Since a minimum of k uv auxiliary nodes are needed to cover any path path(u, v), uv ∈ E (even if both u and v belong to the vertex cover), the number of regular nodes selected is at most k. This concludes the proof.

Finally, we replace each edge in G r with a copy of the gadget ̺ illustrated in Figure 3(c), and denote the resulting graph by G ′ (see figure 3(d)). (We note that each copy is scaled, while maintaining the proportions, to match with the length of the edge it replaces.) Proof of Theorem 6. Consider the graph G ′ , constructed as above from a planar graph G with maximum degree 3 and with a gadget such that |au|+|bu| |ab| = ξ 3. The graph G ′ can be seen as a union of triangles, and it thus a triangle mesh. We prove in the following that G has a vertex cover of size k if and only if G ′ has k + m sources such that its stretch factor is at most ξ. Hence, the vertex cover problem can be reduced in polynomial time to the problem at hand, which concludes the proof.

We first show that if G has a vertex cover of size k, then there is a set of k + m sources in G ′ whose stretch factor is ξ. If G has a vertex cover of size k, then G r has a vertex cover of size k + m by Lemma 8. Recall that this vertex cover of G r can be obtained from the vertex cover of G by adding every other auxiliary node on each edge of G. Let this vertex cover of k + m nodes be the choice of sources in G ′ . Consider a pair p and q of nodes. We consider three cases:

1. p and q belong to the same gadget (the same copy of ̺). Let u, v, a and b denote the nodes of this gadget, as illustrated in Figure 3(c), and suppose, without loss of generality, that u is selected as a source. Then, d(p,u)+d(u,q) d(p,q)

is equal to 1 if p or q coincides with u, it is by definition equal to ξ if (p, q) = (a, b), and it is equal to 3 if p = v and q = a or q = b (see Figure 4 , over all pairs (p, q) in a gadget, is ξ.

2. p and q belong to two adjacent gadgets ̺ 1 and ̺

2 . Let {u 1 , a 1 , b 1 , v 1 , a 2 , b 2 , v 2 }
denote the nodes of the two gadgets. If the node v 1 which belongs to both gadgets is selected as a source then, by symmetry, the analysis of Case 1 yields the same bound of ξ. Otherwise, u 1 and v 2 are sources and, for any two nodes from two different gadgets, the maximum ratio is 2 (the ratio is 1 if p or q is one of the sources and the ratio is 2 in the other cases; see Figure 4(b)). Therefore, again the maximum ratio is at most ξ.

3. p and q belong to neither the same gadget nor to two adjacent gadgets. In this case, at least one of the nodes in a shortest path from p to q is selected as a source, and hence their approximate shortest path equals the geodesic shortest path and min i d(p,s * i )+d(s * i ,q) d(p,q) = 1.

We thus proved that the stretch factor of G ′ , for the selected k + m sources, is ξ.

Conversely, we show that if G ′ has k + m sources such that its stretch factor is at most ξ, then G has a vertex cover of size k. Every gadget must contain at least a source since, otherwise, the vertices a and b of a gadget with no source are such that > ξ. For every gadget in G ′ , if vertex u or v is a source, we select the corresponding vertex in G r , and if vertex a or b is a source, we select any endpoint of the edge of G r corresponding to the gadget. Then, at least one vertex from each edge of the graph G r is selected. Hence, G r has a vertex cover of size k + m, which implies that G has a vertex cover of size k by Lemma 8. This completes the proof.

Conclusions

We analyzed the stretch factor F F P S of approximate geodesics computed as distances through at least one of a set of k sources found using farthest point sampling. We showed that F F P S can be bounded by 2r 2 e F * + 2r 2 e + 8r e + 1, where F * is stretch factor obtained using an optimal placement of the sources and r e is the ratio of the lengths of the longest and the shortest edges in the graph. Furthermore, we showed that it is NP-complete to find such an optimal placement of the sources. Note that in many practical applications r e ≈ 1, which gives some evidence explaining why farthest point sampling has been used successfully for isometry-invariant surface processing.

Figure 2 :

 2 Figure 2: For the proof of Lemma 2.

Lemma 5 .

 5 Let s 1 , • • • , s k be a set of sources returned by the FPS algorithm and s ′ 1 , • • • , s ′ k be an optimal set of sources for the k-center problem. Then max p,q

Figure 3 :

 3 Figure 3: (a) A planar graph G and (b) the grid-embedding of G r . The gadget ̺ replacing the edges of G r and (d) the resulting graph G ′ .

Figure 4 :

 4 Figure 4: For the proof of Theorem 6.

  (a)). Since ξ 3 by definition of the gadget, the maximum of min i d(p,s * i )+d(s * i ,q) d(p,q)

  is argmin s 1 ,...,s k max p d(p, s p ) and thus min s 1 ,...,s k max p d(p, s p ) = max p d(p, s ′

	We now apply the lower bound of Lemma 3 to s * 1 , . . . , s * k which gives
		2 ℓ max	max p	d(p, s * p ) -1 max p,q
							we have
	max p,q	d(p, s i , q) d(p, q)	2r e	2 ℓ min	max p	d(p, s ′ p ) + 1 + 6r e + 1.

By definition, s ′ 1 , . . . , s ′ k is an optimal set of sources for the k-center problem, that p ) max p d(p, s * p ).

Alternatively, we can use the linear time algorithm of[START_REF] Henzinger | Faster shortest-path algorithms for planar graphs[END_REF] for the SSSP problem, as the underlying graph is planar.

Theorem

in[START_REF] Könemann | An approximation algorithm for the edge-dilation k-center problem[END_REF] is stated in a slightly different context but with the notation of that paper, considering πv = Π for every vertex v, the triangle inequality yields the claimed bound.

Note to the reviewers (this will be omitted from the final version of this paper): We have contacted the authors of[START_REF] Könemann | An approximation algorithm for the edge-dilation k-center problem[END_REF], a few weeks ago, and explained them our reasoning as described in the following footnote. We will take into account their response for the final version of this paper.

For a given stretch factor α F * , their algorithm iteratively includes an endpoint of the shortest edge that cannot yet be approximated with a stretch of at most 2α + 1 until no such edges are left. If the solution contains at most k sources, a solution with stretch 2α + 1 2F * + 1 has been found. Their algorithm then essentially does a binary search on the optimal stretch factor F * . However, this search is done in a continuous interval without stopping criteria. Moreover, since it is a priori possible that for any given α < F * , their algorithm returns strictly more than k sources, and that F * may not be exactly reachable by dichotomy, we believe that the stretch factor of FK 2F * + 1 is not ensured.

Acknowledgments

This research was initiated at Bellairs Workshop on Geometry and Graphs, March 10-15, 2013. The authors are grateful to Prosenjit Bose, Vida Dujmovic, Stefan Langerman, and Pat Morin for organizing the workshop and to the other workshop participants for providing a stimulating working environment.


M Wigley

N Kampman

H J Chapman

Benoît Dubacq

M J Bickle

H C Chapman

In situ redeposition of trace metals mobilized by CO2-charged brines

come

Introduction

Geological CO 2 sequestration is proposed as one method for mitigating the effects of fossil fuel burning on global climate [e.g. IPCC ,2 0 0 7 ] . O n eo ft h eo u t s t a n d i n gs c i e n t i fi c questions relates to the potential for injected CO 2 to leak and contaminate overlying freshwater aquifers [Gale,2004;[START_REF] Smyth | AssessingRisk to Fresh Water Resources from Long Term CO 2 DRAFT N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT Injection-Laboratory and Field Studies[END_REF]. CO 2 ,injectedintogeologicalformationsin ag a s e o u so rs u pe r c r i t i c a ls t a t e ,i sl i k e l yt od i s s o l v ei n t oformation fluids at the gas-water contact forming acid brines [Bickle,2 0 0 9 ] . T h e s eb r i n e sc a nm o b i l i z et r a c em e t a l sv i a desorption of metals bound to mineral surfaces, or by dissolution of trace metal-bearing phases [START_REF] Yin | Adsorption of Mercury(II) by Soil: Effects of pH, Chloride, and Organic Matter[END_REF][START_REF] Bradl | Adsorption of Heavy Metal Ions on Soils and Soils Constituents[END_REF]. If these brines were to migrate into overlying freshwater aquifers, there is concern that drinking water supplies could be affected. Models [Wilkin and DiGiulio,2010;Zheng et al.,2009]andexperiments [Little and Jackson,2010;[START_REF] Lu | Potential RiskstoF reshwaterResources as a Result of Leakage from CO 2 Geological Storage: a Batch-Reaction Experiment[END_REF]haveshownthatCO 2 can mobilize trace metals; however field experiments have not shown hazardous levels of contamination [Keating et al.,2 0 1 0 ; [START_REF] Kharaka | Changes in the Chemistry of Shallow Groundwater Related to the 2008 Injection of CO 2 at the ZERT Field Site[END_REF]. Of critical importance for environmental contamination is whether the majority of mobilized metals remain in solution or are re-deposited locally.

CO 2 has been escaping along the Little Wash and Salt Wash grabens for hundreds of thousands of years [Kampman et al.,2 0 0 9 ,2 0 1 2 ; [START_REF] Dockrill | Structural Controls on Leakage from a Natural CO 2 Geologic Storage Site: Central Utah, U.S.A[END_REF]. Near the town of Green River, Utah, red sandstones have been bleached byC O 2 -charged brines which have dissolved hematite grain coatings, with trace metals released from the grain coatings concentrated at the bleaching reaction front [Wigley et al., 2 0 1 2 ] . T h eC O 2promoted reaction fronts are used as an analogue for processes occurring in modern CO 2 injection sites if injected CO 2 were to leak into overlying aquifers. Here we quantify trace

DRAFT

N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT metal enrichments at the bleaching reaction fronts, and use these to calculate corresponding time-integrated fluid fluxes (TIFF; Bickle and Baker ,1 9 9 0 )a n dt h ef r a c t i o no ft r a c e metals which have been re-deposited at the reaction front.

Green River natural analogue site

Within the field area, the Jurassic Entrada sandstone is a well sorted sediment of uniform grain size and mineralogy. Regional burial diagenesis includes the development of Fe-oxide grain coatings [Cullers,1995;[START_REF] Trimble | Geology and Uranium-Vanadium Deposits of the San Rafael River Mining Area[END_REF], giving the rock its red color. Near the town of Green River, Utah, the Entrada sandstone has been bleached a pale yellow by diagenetic fluids which have dissolved the Fe-oxide grain coatings [Wigley et al., 2 0 1 2 ] . B l e a c h i n g occurs at the base of the formation in a broad domal structure,e x p o s e do v e ra1 -2 k m east-west section along the Salt Wash graben, at the crest of the Green River anticline (Fig. 1a).

Several lines of evidence suggest that the bleaching fluid wasaC O 2 -charged brine with minor methane, distinguishing the Green River bleaching from hydrocarbon-related bleached sandstones elsewhere on the Colorado Plateau [START_REF] Wigley | Fluid-Mineral Reactions and Trace Metal Mobilization in an Exhumed Natural CO 2 Reservoir[END_REF][START_REF] Beitler | Bleaching of Jurassic Navajo Sandstone on Colorado Plateau Laramide highs: Evidence of exhumed hydrocarbon supergiants?[END_REF][START_REF] Beitler | Fingerprints of Fluid Flow: Chemical Diagenetic History of the Jurassic Navajo Sandstone, Southern Utah, U.S.A[END_REF]. Carbonates associated with the Green River bleaching have a heavier carbon isotopic composition than documented occurrences of hydrocarbon related bleaching [START_REF] Wigley | Fluid-Mineral Reactions and Trace Metal Mobilization in an Exhumed Natural CO 2 Reservoir[END_REF]. These values are similar to those measured in the actively leaking CO 2 -system near Green River [Assayag et al., petrographically linked to the bleaching suggests that the bleaching fluid was a low temperature (≈27 • C) brine (salinity 2-5-7 wt%) with the vapor phase consisting dominantly of CO 2 , with minor (0-28%) methane present in three inclusions [START_REF] Wigley | Fluid-Mineral Reactions and Trace Metal Mobilization in an Exhumed Natural CO 2 Reservoir[END_REF]. The magnitude of the Fermi-diad separation in the Raman spectra ofC O 2 is calibrated for the pressure of CO 2 in the vapor phase at high temperatures [Rosso and Bodnar ,1 9 9 5 ] .

An extrapolation to calculated temperatures during bleaching suggests that the bleaching fluid had a CO 2 partial pressure of 1.8-25.9±9.7 bar (Figure 2). 2. Methods

Sampling and analytical methods

Profiles from bleached to unbleached sandstones were cored using an electric diamondedged 40mm core plug drill, with samples spaced at 2-15cm intervals. Core plugs near to the bleached-red contact were sub-sampled every 4mm using a 2mm tungsten-carbide drill bit. Between drilling, the bit and core were cleaned with compressed air. Samples were dissolved in HF and HNO 3 using an open beaker method, and analyzed by inductively coupled plasma-atomic emission spectroscopy (major elements), and ICP-mass spectrometry (trace elements) in Cambridge, UK.

Mass of metals in enriched zones

The mass of a metal in the enriched region is calculated by integrating across the enriched zone, assuming a stepped concentration profile. Uncertainties are estimated by propagating the uncertainties on the background and enriched compositions. The background compositions are taken as the average within the bleached zone (bleached sandstone), and the average composition of red sandstone near to the bleached fracture (unaltered sandstone; 1σ uncertainty, background compositions are given in supplementary table 4). The uncertainty in peak height is propagated 1σ instrument uncertainty, 1σ of repeat analyses and recovery relative to USGS standard Sco-1 (Supplementary table 1). The range in the fraction of each metal re-adsorbed is calculated from a Monte-Carlo simulation (n=1000) where each parameter listed above is allowed to vary within its range in uncertainty. The range in the fraction of each metal re-deposited is given by 1σ about the mean of the resulting distribution.

DRAFT N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT

Sequential extraction

The nature of the minerals hosting the re-deposited metals wasd e t e r m i n e db yas equential extraction procedure (modified from [START_REF] Tessier | Sequential Extraction Procedure for the Speciation of Particulate Trace Metals[END_REF] and [START_REF] Förstner | Accumulative Phases for Heavy Metals in Limnic Sediments[END_REF]) ). An average clay composition is used based on the major element chemistry of the <2µm fraction (Al 1.4 Ca 0.2 Fe 0.5 K 0.8 Mg 0.2 Na 0.2 Si 3.8 O 10 (OH) 2 ;r e c a l c u l a t e dt oa n1 2o y y g e n formula). Fe and Ti-oxides, pyrite and quartz are assumed to exist as pure phases.

applied

Results and discussion

Phases dissolved during sequential extraction

The results of the modal decomposition of sequential leaching solutions are shown in

Fraction of mobilized metals re-deposited

The minimum fraction of each trace metal re-deposited can be calculated from bleached halos surrounding fractures where the displacement distance of the reaction front is given by the perpendicular distance of the front from the fracture. The theoretical mass of a metal mobilized along a 1-dimensional path is given by the difference between the background concentrations in the red and bleached sandstones, and the displacement distance. The ratio of the predicted to observed mass of metali nt h ee n r i c h e dz o n eg i v e s am e a s u r eo ft h ef r a c t i o no ft h a tm e t a lr e -d e po s i t e da tt h er e action front.

Some metals analyzed (Cu, V1; Zn, V2 and V3) exhibit re-deposited fractions close to unity, with the adsorbed mass being within error of the total mass dissolved (Table 2). Others (Cu, V2 and V3; Co, V1) show enrichment relative to the mass dissolved by propagation of the reaction front from the fracture. We consider that the true value lies in the upper range of calculated results as the sampling may have missed peaks in concentration, thus underestimating the mass of metal that has been re-deposited. be enriched in some metals and these, along with those mobilized by the reaction front moving away from the fracture, are all re-deposited.

Metals in solution exist either as free metal cations or metal-anion complexes (e.g. with OH -, HCO - 3 or Cl -). Speciation modeling using the calculated range in P CO2 and salinity and assuming that the fluid is buffered to equilibrium with calcite, suggests that Sn, Zn and Co exist as uncomplexed species, with a minor fraction of Co and Zn complexing with chloride and Zn also forming carbonate complexes (Supplementary figure 2); Cu speciation is dominated by chloride species. Realistic variations in methane concentrations do not significantly change calculated metal speciation profiles. Modelling was conducted using PHREEQC [START_REF] Parkhurst | User's Guide to PHREEQC (version 2) -A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations[END_REF] with the llnl database [START_REF] Johnson | Database 'thermo.com[END_REF].

These calculations show that metal complexation does not playam a j o rr o l ei nt h er edeposition process, as metals complexing with hydroxyl, carbonate and chloride anions as well as free metal cations are all re-deposited at the CO 2 -promoted reaction fronts.

Transp ort of mobilized metals is thus controlled by the distance moved by CO 2promoted reaction fronts, rather than by the distance moved by the fluid. The difference between these length scales is controlled by the stoichiometry and kinetics of the fluid-rock reactions, and therefore the chemistry and mineralogy of the host formation. In non-reactive reservoirs mobilized metals will be transported further; however such reservoirs are less likely to contain significant quantities of contaminants.

Vertical transport distances across sub-horizontal reaction fronts

The observation that the metals removed from the bleached sandstone are quantitatively re-deposited can be used to calculate minimum transport distances perpendicular to the horizontal reaction fronts. The length of rock bleached by the fluid interacting DRAFT N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT with horizontal reaction fronts has been calculated using the difference in background concentration in bleached and unbleached samples (mol.m -3 ) and the mass of each metal deposited at the reaction front (mol.m -2 ), assuming 100% re-deposition (l=M/∆C; where l is the length of rock bleached, M is the mass of the metal in the enrichment (mol/m 2 reaction front) and ∆C is the difference in background concentrations between red and bleached sandstone). Some of the metals may have been removed by the sub-horizontal advection of the gravity current although this flux may decrease near the margin of a gravity current. If so the calculation underestimates the transport distance. Four transects across horizontal reaction fronts were analyzed (H1-4; localities in supplementary table

2), and the length of rock bleached for each unit area of reaction front is low (<0.35m; Fig. 5), consistent with the CO 2 -rich fluid migrating laterally as a gravity current with vertical motion is dominated by diffusive processes except intheregionofopenfractures.

Potential impact on water quality

The TIFF (the volume of fluid that has interacted with the reaction front) has been calculated from the retardation factor (R) and the distance the front is displaced (Table 2). R is defined as the volume of fluid required to displace the reaction front through au n i tv o l u m eo fr o c k ,a n di se s t i m a t e df r o mt h er e a c t i o ns t o i chiometry and change in brine-CO 2 concentration across the reaction front [START_REF] Wigley | Fluid-Mineral Reactions and Trace Metal Mobilization in an Exhumed Natural CO 2 Reservoir[END_REF]. TIFF is then given by R×distance, i.e. the volume of fluid required to displace the reaction front through the observed volume of rock.

Metal concentrations in the fluid (C m )a st h er e a c t i o nf r o n ti sd i s p l a c e da r eg i v e nb y :

C m = F m F f (1)
DRAFT N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT where F m is the flux of the metal (mol.m -2 .s -1 ), and F f is the fluid flux (m.s -1 ). The time averaged metal flux is given by the mass of metal in the enrichment (mol/m 2 reaction front); and the corresponding F f is given by the TIFF (Table 2).

Concentrations of

Concluding remarks

The reaction fronts studied here, where fluid velocities are low and diffusive processes are important, have a high Damköhler number (N D ;t h er a t i oo fr e a c t i o nr a t et or a t eo f advection) and a low Peclet number (Pe; the ratio of advectivet od i ff u s i v et r a n s p o r t ) .

The high injection pressures in CO 2 -storage reservoirs will result in high rates of advection and consequently such CO 2 -promoted reaction fronts will have a lower N D and a higher Pe compared to those documented here. However much of the concern over contaminant mobilization relates to migration of CO 2 from storage reservoirs into overlying aquifers where the groundwater flow rates are likely to be considerably lower (Fig. 6), being dominated by regional gradients in hydraulic head, such thatC O 2 -promoted reaction fronts have a transport characteristic comparable to those studied here. In this context, our results show that mobilized metals are likely to be re-deposited over small length scales 0.00 (0) 0.00 (0) 0.00 (0) 0.00 (0) 0.00 (0) a Numbers in parentheses are 1 σ uncertainties for the last digit, calculated from the magnitude of the misfit to bulk rock values.

DRAFT

N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT

  2 0 0 9 ] ,b u ta r ei n c o n s i s t e n tw i t hi s o t o p ically light methane bearing fluids. Mobilization of low pH-soluble trace metals is also inconsistent with methane-rich fluids which have pH ≈7.5 when saturated with hematite (Supplementary figure 1). The occurrence of the bleaching at the base of the stratigraphy also rules out buoyant methane bearing fluids as the cause of the alteration. Raman DRAFT N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT spectroscopy of fluid inclusions hosted in aragonite, quartzo v e r g r o w t h sa n dg y p s u m ,

  The upper bleached-red contact is well exposed and is sharp and sub-horizontal[Wigley et al.,2 0 1 2 ]e x c e p ti nr e g i o n so ff r a c t u r e sa n dg r a n u l a t i o ns e a m s , where differences in porosity and tortuosity enhance or retard vertical transport of the reaction front, respectively (Fig.1b-c; Wigley et al.,I nP r e s s ) . T h el o w e rc o n t a c ta p p e a r sl e s ss h a r pb u ti s rarely exposed. To the East the bleached zone terminates as a series of meter scale fingers into red sandstones, whilst the western extent is not exposed. It is thought that the CO 2 -rich fluid migrated as a gravity current, with the main components of flow to the north along the plunge of the anticline, and to the east, down dip and in the direction of regional groundwater flow[Kampman et al.,2 0 0 9] ,w i t hv e r t i c a lt r a n s p o r td o m i n a t e d by diffusion except in areas containing faults and open fractures (Fig. 1b; Wigley et al.2012).Analyses of profiles across the upper sub-horizontal contact between red and bleached sandstone show enrichments in bulk rock trace metal content within 10 cm of the contact (Fig.3a-f; raw data and sample localities are given in supplementary tables 1 and 2), thought to originate from release of trace metals as grain coating Fe-oxides were dissolved,DRAFTN o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT with subsequent re-deposition at the reaction front [Wigley et al.,2 0 1 2 ] . T h e s et r a c e metals are soluble at low pH, becoming insoluble as pH rises across the reaction front. The scale of enrichments is determined by the volume of fluid that interacted with the reaction front, which had previously reacted with unalteredr oc k . Bleaching also occurs around open fractures (Fig 1c), which are mineralized with carbonate and Fe-oxides, and exhibits the same geochemical, isotopic and petrological signatures as the sub-horizontal bleaching contacts [Wigley et al.,2 0 1 2 ] ,i n d i c a t i n gac o m m o np a rent fluid reservoir and continuity in reactive-transport processes. The fracture-associated bleaching reaction front began within the plane of the fracturea n dm i g r a t e dl a t e r a l l yb y diffusion and advection of the fluid. The asymmetry of the bleached halos is due to the diffusion away from the fractures being imposed on uni-directional horizontal regional flow in the aquifer. Spikes in heavy metal concentration are observed near to the red-bleached transition (Fig. 3a-f), similar to those across horizontal contacts. These observed enrichments of trace metals away from bleached fractures may be related to deposition of metals released by the bleaching reaction front as it moves outwards from the fracture, and/or metals present in the bleaching fluid as it enters the fracture. Key questions include: (1) what was the original distribution of the metals? (2) How, and in which phases, are the metals re-deposited? (3) What fraction of the mobilized metals is subsequently re-deposited? (4) Does the concentration of metals in the fluid reach levels hazardous for human health? DRAFT N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT X-7

  to a selection of unaltered and bleached samples and to metal-enriched samples from sub-horizontal reaction fronts. Leaching was performed at room temperature in 5 steps: (1) Water rinse: 1.5mL de-ionized 18MΩ water with continuous agitation for 2 hours. (2) Exchangeable fraction: 1.5mL 1M sodium acetate solution (NaOAc) at pH 8.2, with continuous agitation for 3 hours. (3) Carbonates: 1.5mL 1M NaOAc adjusted to pH 5 with acetic acid (HOAc), with continuous agitation for 7 hours, repeated 3 times. (4) Oxides: 1.5mL 0.1M oxalate buffer adjusted to pH 3 with oxalic acid. Occasional agitation for 54 hours, repeated 3 times. (5) HCl rinse (remaining carbonate/oxides): 1.5mL 16M HCl with continuous agitation for 2.5 hours. After each leach the sample was centrifuged for 30 minutes at 3000 rpm, and the liquid pipetted carefully into a cleaned centrifuge tube. The sample was then rinsed with MΩ water, re-centrifuged, and the rinse added to the leachate solution. After the final leaching step, the residue was dried down and dissolved in a combination of HF and HNO 3 in an open beaker method. The leachate solutions were evaporated to dryness and re-dissolved in a dilute nitric acid solution. The number of moles of a mineral phase dissolved in each sequential leaching step was determined by a modal decomposition of the compositions of the leachate solutions. This was done by minimizing the difference between theoretical ando b s e r v e dc o m p o s i t i o n s by varying the moles of each mineral dissolved. Mineral compositions of plagioclase, K-feldspar and carbonates are taken from electron microprobea n a l y s e s( s u p p l e m e n t a r y DRAFT N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT table 5

  ates. The oxide leach removes up to 60% of total Fe with a greater fraction released in bleached samples, reflecting the dissolution of secondary coarse oxide phases rather than grain coating oxy-hydroxides. Approximately 30% of the minerals dissolved are Fe-oxides, along with some clays and residual dolomite. A moderate correlation (R 2 =0.34) between Fe and S suggests dissolution of minor pyrite; however a Fe/S ratio of ≈100 precludes pyrite dissolution as the dominant control on the Fe content of the leachate solution. The hydrochloric acid leach dissolves residual carbonate and oxide phases, and minor silicates.AlargerfractionoftotalF eisreleasedinunalteredsamplesinthisstep,reflectingpartial dissolution of the grain coating oxy-hydroxide phases. 3.1.1. Location of re-deposited metals DRAFT N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT The results of the sequential extraction experiments are typified by cobalt (Co; Fig. 4). In both bleached and unbleached samples where no enrichment is observed, Co is released from carbonate, oxides and by dissolution of the residue (Fig.4 a ) ,i n d i c a t i n gt h a tt h e trace metals are equally distributed between carbonate, coarse oxides and grain coating phases. In samples with bulk rock Co enrichment, the Co released in the steps dissolving carbonate and oxide phases is two orders of magnitude greatert h a nb a c k g r o u n dv a l u e s (Fig. 4b). Re-deposited trace metals are therefore incorporated into the secondary oxide and carbonate phases growing at the bleaching reaction front[ Wigley et al.,2 0 1 2 ] .

  Fluids entering the open fractures have previously reacted with unaltered rock and therefore may DRAFT N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT X -11

  Cu and Zn in solution are an order of magnitudeb e l o wt h eW o r l d Health Organization (WHO) drinking water limit (Cu: 2 mg/L, Zn: 5mg/L; WHO,2011) and Sn concentrations are well below acceptable limits (150 mg/L; WHO, 2003). No guidelines exist for Co; however calculated Co concentrations are 10 times the average composition of US freshwater [ATSDR,1 9 9 2 ]( s e ea p pe n d i xAf o raw o r k e de x a m p l e ) .

  in secondary carbonate and oxide phases precipitating at theC O 2 -promoted reaction fronts. Fluxes of CO 2 -charged waters of order of 25 times the volume of an aquifer of DRAFT N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT

Figure 1 .

 1 Figure 1. Illustration and outcrop photos of the exhumed CO 2 reservoir. (a) Bleaching

Figure 3 .

 3 Figure 3.Concentration profiles for selected metals across horizontal reaction fronts (a-c,

Figure 4 .Figure 5 .

 45 Figure 4. Results of sequential leaching for Co. a) Mass of Co released in 3 samples without:

Figure 6 .

 6 Figure 6. Schematic description of a CO 2 injection site. If CO 2 migrates into overlying

table 1 .

 1 The water rinse predominantly dissolves soluble salts, releasing up to 86% of available S with a Ca/S ratio of ≈ 1. Modal analysis on water rinse solutions shows that

75% of all minerals dissolved in this step are gypsum (table

1

).T h ee x c h a n g e a b l el e a c h removes residual soluble salts and minor amounts of calcite and silicates. It also removes loosely adsorbed cations: in some samples up to 60% of total Sri sr e m o v e di nt h i ss t e p .

Up to 80% total Ca and 60% total Mg are removed in the carbonate leach indicating dissolution of calcite and dolomite. Approximately 95% of minerals dissolved are carbon-

Table 1 .

 1 Mol fractions of minerals released in each sequential leaching step a

			Average mol fraction released
		MΩ water rinse Exchangeable leach Carbonate leach Oxide leach HCl leach
	Quartz	0.00 (0)	0.02 (1)	0.00 (0)	0.03 (1)	0.00 (0)
	K-feldspar	0.10 (2)	0.17 (6)	0.01 (1)	0.12 (5)	0.01 (1)
	Plagioclase	0.00 (0)	0.00 (0)	0.00 (0)	0.00 (0)	0.01 (0)
	Clay	0.07 (1)	0.19 (7)	0.01 (0)	0.34 (12)	0.04 (1)
	Calcite	0.04 (1)	0.23 (8)	0.63 (5)	0.00 (0)	0.84 (24)
	Dolomite	0.05 (1)	0.05 (2)	0.18 (1)	0.17 (6)	0.04 (1)
	Fe-dolomite	0.00 (0)	0.01 (0)	0.14 (1)	0.00 (0)	0.00 (0)
	Hematite	0.00 (0)	0.00 (0)	0.01 (0)	0.29 (11)	0.06 (2)
	Pyrite	0.00 (0)	0.00 (0)	0.00 (0)	0.04 (2)	0.00 (0)
	Gypsum	0.73 (12)	0.34 (13)	0.02 (0)	0.00 (0)	0.00 (0)
	TiO 2

Acknowledgments. MW would like to thank Jason Day for advice with sequential leaching and geochemical analyses. Carbon research at Cambridge is supported by Natural Environment Research Council grant NE/F004699/1, part of the UK CRIUS (Carbon Research into Underground Storage) consortium. Niko Kampman is a recipient of financial support from Shell International. The manuscript was significantly improved as a result of comments by two anonymous reviewers.

The metal concentration is given by the metal flux divided by the fluid flux (Eqn. 1).

By integrating over the duration of the bleaching event, this reduces to:

where M is the mass of metal at the reaction front (Kg/m 2 ) and TIFF is the time integrated fluid flux (m 3 /m 2 ), which is given by: R×X, where X is the distance moved by the reaction front.

The mass of copper deposited at reaction front V1 is 0.00152 Kg/m 2 ,a n dt h er e a c t i o n front has moved 0.506 m. The corresponding TIFF is therefore 13m 3 /m 2 and the concnetration of copper in the fluid as the reaction front propagates is 0.000312 kg/m 3 or 0.31 mg/L.

X -25

DRAFT

N o v e m b e r2 3 ,2 0 1 2 ,1 1 : 2 0 a m DRAFT


Xu Wu

Jun Zheng
email: junzheng@seu.edu.cn

Hui Tong

Nathalie Mitton

DAYcast: A Dynamic Transmission Delay Based Broadcast Protocol for Vehicular Ad Hoc Networks

Keywords: broadcast, broadcast storm, DAYcast, effective neighbor, routing, VANET I

Broadcasting is an effective routing paradigm for data dissemination in vehicular ad hoc networks (VANETs). One concern that arises with broadcasting is the broadcast storm problem, which would cause node contentions and data collisions, and thus degrade the transmission efficiency of a network. This paper proposes a Dynamic trAnsmission delaY based broadcast (DAYcast) protocol for a VANET. To alleviate the effect of the broadcast storm and improve the transmission efficiency of the network, DAYcast only allows the effective neighbors of a source vehicle to broadcast a received data packet and the selection of the effective neighbors are based on the position information on the one-hop neighbors of the source vehicle. Meanwhile, it allows each effective neighbor to wait a certain transmission delay before it broadcasts a received packet. The transmission delay of an effective neighbor depends on the distance between the neighbor and the source vehicle, and the number of effective neighbors of the source vehicle. Simulation results show that DAYcast can effectively improve the network performance in terms of network reachability and the successful delivery ratio as compared with existing weighted p-persistence broadcasting (WPB) and slotted 1-persistence broadcasting (SPB).

INTRODUCTION

Vehicular ad hoc networks (VANETs) have widely been considered as a promising networking technology for supporting intelligent transportation systems (ITS) [START_REF] Courtney | A broad view of ITS standards in the US[END_REF]. A VANET consists of a number of vehicles moving on highways or roads in an urban area, which may be equipped with sensing devices, radio transceivers, and positioning systems, and can thus provide a variety of applications, such as traffic control, environment monitoring, and inter-vehicle communication [START_REF] Kong | A collaboration-based hybrid vehicular sensor network architecture[END_REF]. In a VANET, a vehicle or roadside unit (RSU) needs to disseminate relevant information or data to one or multiple destinations via intermediate vehicles or RSUs. For this purpose, a routing protocol is needed to route the information or data to its destinations. Compared with traditional mobile ad hoc networks (MANETs), however, a VANET have some unique characteristics, such as highly dynamic topology, intermittent connectivity, and predictable mobility model, which presents a big challenge for the design of a routing protocol. Therefore, routing becomes an important issue in the design of a VANET.

Broadcasting is an effective routing strategy for disseminating data in VANETs. In broadcasting, a source node broadcasts a data packet to all its neighbor nodes, which further broadcast the packet until the packet reach its destinations. Compared with unicasting, this routing paradigm has the advantage of simplicity in implementation and can achieve a higher packet delivery ratio. However, broadcasting would introduce the broadcast storm problem which is caused by unlimited or excessive dissemination of the same copy of a data packet in a network. This problem would result in the network to be flooded by a large amount of duplicate packets, which would cause severe node contentions and data collisions, and thus degrade the transmission efficiency of the network. For this reason, it is critical to address the broadcast storm problem in the design of a broadcast protocol in order to improve the transmission efficiency of the network.

In this paper, we propose a Dynamic trAnsmission delaY based broadcast (DAYcast) protocol for data dissemination in a VANET. To alleviate the effect of the broadcast storm and improve the transmission efficiency of the network, DAYcast only allows the effective neighbors of a source vehicle to broadcast a received data packet and the selection of the effective neighbors are based on the position information on the one-hop neighbors of the source vehicle. Meanwhile, it allows each effective neighbor to wait a certain transmission delay before it broadcasts a received packet. The transmission delay of an effective neighbor depends on the distance between the neighbor and the source vehicle, and the number of effective neighbors of the source vehicle. Simulation results were conducted to evaluate the performance of the proposed DAYcast protocol.

This article is organized as follows. Section II reviews related work on broadcast protocols for VANETs. Section II presents the proposed DAYcast protocol. Section IV shows simulation results to evaluate the performance of DAYcast. Section V concludes this paper.

II. RELATED WORK

Broadcasting has been widely studied for VANETs and a variety of broadcast protocols have been proposed in the literature [START_REF] Ma | Design and analysis of a robust broadcast scheme for VANET safety-related services[END_REF][START_REF] Yiannos | Speed adaptive probabilistic flooding in cooperative emergency warning[END_REF][START_REF] Korkmaz | Urban multi-hop broadcast protocol for intervehicle communication systems[END_REF][START_REF] Barradi | Highway multihop broadcast protocols for vehicular netwroks[END_REF][START_REF] Wisitpongphan | Broadcast storm mitigation techniques in vehicular ad hoc networks[END_REF][START_REF] Wisitpongphan | On the routing problem in disconnected vehicular ad hoc networks[END_REF][START_REF] Sok-Ian | SCB: store-carry-broadcast scheme for message dissemination in sparse VANET[END_REF][START_REF] Jinyoun | Fast broadcast at the intersection in VANET[END_REF][START_REF] Tonguz | DV-CAST: a distributed vehicular broadcast protocol for vehicular ad hoc networks[END_REF]. In [START_REF] Ma | Design and analysis of a robust broadcast scheme for VANET safety-related services[END_REF], Ma et al. proposed a distributive cross-layer scheme for the design of the control channel in Dedicated Short Range Communication (DSRC) with three levels of broadcast services, including cross-layer message priority setting, dynamic receiver-oriented packet repetitions, and farthest relay with distance-based AD timer for multihop broadcast. In [START_REF] Yiannos | Speed adaptive probabilistic flooding in cooperative emergency warning[END_REF], Yiannos et al. proposed a novel speed adaptive probabilistic flooding (SAPF) algorithm, which makes a decision to rebroadcast a message based on the speed of a vehicle and relies only on local information. In [START_REF] Korkmaz | Urban multi-hop broadcast protocol for intervehicle communication systems[END_REF], Korkmaz et al. proposed an urban multi-hop broadcast (UMB) protocol, which only allows the furthest vehicle from the transmitter to rebroadcast the packet by using a blackburst contention approach. In [START_REF] Barradi | Highway multihop broadcast protocols for vehicular netwroks[END_REF], Barradi et al. proposed a new 802.11 based vehicular multi-hop broadcast protocol, called highway multihop broadcast (HMB), which selects the farthest vehicle with the least speed deviation from the source to forward and acknowledge broadcast frames. Moreover, HMB introduces ICTB/FCTB handshake to resolve the hidden node problem and passive forwarder selection to solve the storm problem. In [START_REF] Wisitpongphan | Broadcast storm mitigation techniques in vehicular ad hoc networks[END_REF], Wisitpongphan et al. quantified the impact of the broadcast storm in terms of message delay and packet loss rate and proposed three probabilistic and timer-based broadcast suppression techniques. However, the solutions proposed in [START_REF] Yiannos | Speed adaptive probabilistic flooding in cooperative emergency warning[END_REF][START_REF] Korkmaz | Urban multi-hop broadcast protocol for intervehicle communication systems[END_REF][START_REF] Barradi | Highway multihop broadcast protocols for vehicular netwroks[END_REF] attempt to alleviate the broadcast storm from the perspective of MAC layer. The solutions proposed in [START_REF] Wisitpongphan | Broadcast storm mitigation techniques in vehicular ad hoc networks[END_REF] are only suitable for dense networks and cannot adapt to the density changing scenarios.

In [START_REF] Wisitpongphan | On the routing problem in disconnected vehicular ad hoc networks[END_REF], Wisitpongphan et al. developed a statistical traffic model on the data collected to indicate the characteristics of sparse VANETs. In [START_REF] Sok-Ian | SCB: store-carry-broadcast scheme for message dissemination in sparse VANET[END_REF], Sou et al. proposed a store-carrybroadcast (SCB) scheme to assist message dissemination by broadcasting over a specific road segment instead of a single vehicle in sparse VANET. An opposite vehicle is used to disseminate the messages to oncoming vehicles traveling on the reverse lane by broadcasting in SCB. In [START_REF] Jinyoun | Fast broadcast at the intersection in VANET[END_REF], Cho et al. proposed an efficient way to broadcast a safety message to all directions at an intersection with a short delay and without a collision. In [START_REF] Tonguz | DV-CAST: a distributed vehicular broadcast protocol for vehicular ad hoc networks[END_REF], Tonguz et al. proposed DV-CAST, which relies only on local topology information for handling broadcast messages in VANETs. DV-CAST can operate in all traffic regimes, including sparse and dense scenarios, but is relatively complicated.

III. DYNAMIC TIMESLOT-BASED BROADCAST PROTOCOL

In this section, we present the proposed DAYcast protocol for VANETs.

A. Network Model

We consider a straightway scenario without intersections, where a number of vehicles are moving on the road, which has one or more lanes in each direction, as shown in Fig. 1. The source node can be a vehicle or an RSU. We only consider the vehicles moving in one direction and do not turn to the vehicles moving on the other direction for data forwarding. Moreover, we do not consider car overtaking behavior. We assume that each vehicle is equipped with a GPS device, which can obtain the geographic information on the vehicle. Each vehicle has the same transmission radius. Only when a vehicle is within the transmission radius of another vehicle can the two vehicles communicate with each other.

B. Routing Strategy

The main idea of the DAYcast protocol is to only allow the effective neighbors of a source vehicle to broadcast a received data packet and allow each intermediate vehicle or node to wait a certain transmission delay before it further broadcasts the received packet in order to reduce node contentions and data collisions, and thus alleviate the effect of the broadcast storm. The selection of the effective neighbors is based on the position information on the onehop neighbors of the source vehicle. The transmission delay of an effective neighbor depends on the distance between the neighbor and the source vehicle, and the number of effective neighbors of the source vehicle. Since the broadcast storm is usually caused by a large number of nodes which simultaneously further broadcast a received packet and contend for a common transmission channel, it is expected to reduce the number of broadcasting nodes and avoid the broadcasting nodes to broadcast data packets simultaneously.

In a VANET, a vehicle is moving at a high speed and its position changes frequently. To support DAYcast, each vehicle needs to obtain its position information (i.e., latitude and longitude) using its GPS device and disseminate the position information to its one-hop neighbors periodically. For this purpose, DAYcast allows a Hello message, which is typically used in a routing protocol, to piggyback the position information in order to reduce the control overhead in the network. According to the received Hello messages, each node that receives a data packet can know the number of its effective one-hop neighbors based on the position information contained in the Hello messages. When a node has a data packet to broadcast, it will first add the number of effective neighbors as well as its own position information (i.e., latitude and longitude) in the data packet, and then broadcast the packet. After a node receives the data packet, it can know the number of its last hop's effective neighbors and determine whether it needs to broadcast and when to broadcast the received packet.

C. Protocol Description

The DAYcast protocol consists of two components: effective neighbor selection and dynamic transmission delay calculation.

a) Effective neighbor selection

The purpose of "effective neighbor" selection for a node that has a data packet to broadcast is to select a group of its one-hop neighbors in the direction of packet transmission as the "effective neighbors" for further broadcasting the packet. To alleviate the broadcast storm, DAYcast does not allow all one-hop neighbors of a source vehicle to further broadcast the packet. Instead, it only allows those one-hop neighbors in the direction of packet transmission to further broadcast the packet, which are called "effective neighbors".

To select the effective neighbors of a node, DAYcast employs a Hello message to exchange position information between different one-hop neighbors of the node. The Hello message contains the node ID, and the node's position information (i.e., latitude and longitude), which can be obtained by using the GPS device of the node. After a vehicle receives the Hello messages from all its one-hop neighbors, it extracts the ID, and the latitude and longitude of the source vehicle, compares the source's position with its own position, and determines the number of its effective neighbors. Fig. 2.

Illustration of effective neighbor selection Fig. 2 gives an example for illustrating effective neighbor selection. Assume that a source vehicle A within the dotted circle receives a data packet from its left vehicles and needs to broadcast the packet to its right vehicles. Within its transmission range, vehicle A has seven one-hop neighbors. In this case, the seven neighbor nodes of vehicle A will send Hello messages to vehicle A. But the vehicles to the left of vehicle A have already received the packet and are not in the direction of further transmission relative to vehicle A. Therefore, vehicle A will select the three one-hop vehicles on its right side as its effective neighbors based on the information contained in the Hello messages, i.e., the number of effective neighbors is 3 in this case. In addition to the number of its effective neighbors, vehicle A can also know the maximum distance to its effective neighbors, which is the distance between A and B in this case, based on the Hello messages it has received. When it broadcasts a data packet, it will piggyback its ID, position information, the number of its effective neighbors, and the maximum distance to its effective neighbors.

b) Dynamic transmission delay Calculation

The purpose of dynamic transmission delay calculation is to dynamically calculate a transmission delay for each effective neighbor to further broadcast its received data packet. Since a source vehicle usually has more than one effective neighbor. To avoid node contentions and data collisions, it is expected that all the effective neighbors do not broadcast their received packet simultaneously. For this purpose, we allow each effective neighbor to broadcast the received packet after a certain transmission delay, which depends on the distance between an effective neighbor and the source vehicle. Obviously, this transmission delay is not equal for different effective neighbors.

Let N j denotes the number of effective neighbors of source vehicle j, D ij is the distance between vehicle i and vehicle j, D max is the maximum distance between an effective neighbor vehicle and a source vehicle. Thus, the transmission delay for effective neighbor i, denoted by T i , can be calculated as

max max ij i j D D T N D τ - = ×
(1) where ( ) ( )

2 2 ij i j i j la la lo lo D = + - - ,
and (la i , lo i ) is the latitude and longitude of vehicle i, and τ is a unit time. Note that N j , D max , and (la i , lo i ) can be extracted from the received data packet.

c) Protocol procedure

The major procedure of the proposed DAYcast protocol is given in Fig. 3 The procedure of DAYcast protocol 1) When a source vehicle has a data packet to broadcast, it first determines the number of its effective neighbors as well as the maximum distance between an effective neighbor vehicle and the source vehicle based on the information contained in the Hello messages from its one-hop neighbors. 2) Then it adds the number of its effective neighbors, the latitude and longitude of its own, and the maximum distance in the data packet, and broadcast the data to its one-hop neighbors. 3) If a vehicle has received the data packet and broadcast the packet before, it will drop the packet. 4) If a vehicle receives a data packet for the first time, it means that the vehicle is an effective neighbor of the source vehicle. In this case, it will calculate its A Transmission direction B transmission delay for further broadcasting the packet based on Eq. ( 1). During the transmission delay, if the vehicle receives the data packet again, it will compare the transmission delay with the latest reception time. If the time difference is smaller than a predefined threshold, the vehicle will broadcast the packet at the end of the transmission delay to ensure the network transmission ratio. Otherwise, the vehicle will give up waiting and broadcasting. 5) After a vehicle receives a data packet, it will not explicitly send an acknowledgment message back to the source vehicle. Instead, the source vehicle can know the reception status of its effective neighbors by overhearing the data packet further broadcast by its effective neighbors. 6) If the source vehicle knows that the data packet has already been further broadcasted by its effective neighbors, it will remove the packet from its buffer. Otherwise, it will broadcast the data packet again.

τ ×         - = j max ij max i N / D D D T in transmission delay T i { V i receive P at t 2 if T i -( t 2 -t

IV. PERFORMANCE EVALUATION

In this section, we evaluate the performance of the DAYcast protocol through simulation results. The simulation experiments are conducted using NS2.35 [START_REF]NS-2[END_REF]. In the simulation experiments, we consider a 2km straightway scenario with a number of vehicles ranging from 40 to 250. The source node is a RSU located at the starting end of the straightway. The vehicle driving track is generated using VanetMobiSim [START_REF] Harri | VanetMobiSim: generating realistic mobility patterns for VANETs[END_REF]. IEEE 802.11 with a transmission rate of 2Mbps and a transmission range of 250m is used as the underlying MAC protocol. The time to live (TTL) of a data packet is set to 100 hops. The running time of each simulation experiment is 500 seconds. All simulation results are averaged over 20 runs. The parameters used in the simulations are summarized in Table I. For performance evaluation, we put a virtual node at the 2 km end of the straightway. We compare the proposed DAYcast protocol with weighted p-persistence broadcasting (WPB) and slotted 1-persistence broadcasting (SPB) proposed in [START_REF] Wisitpongphan | Broadcast storm mitigation techniques in vehicular ad hoc networks[END_REF]. The reforwarding probability in WPB is supposed to be 0.5.The performance metrics used in the evaluation are defined as follows:

 Network reachability: The number of data packets received by the virtual node to the total number of data packets broadcast by the source node;  Successful delivery ratio: The number of vehicles successfully receiving the data packets sent by the source node to the total number of vehicles in the network;  Average delivery delay: The average time taken for a data packet to be delivered from the source node to the virtual node. Fig. 4 shows the network reachability with DAYcast, WPB, and SPB under different node density, respectively. It is seen that when the node density is 20 vehicles/km, the network reachability is very low. This is because in this case the network is not fully connected. With the node density increasing, this performance can reach 80% when the density is between 30 vehicles/km to 75 vehicles/km. When the density continues to increase, the network will become fully connected and the broadcast storm problem will become more serious. In this case, the network reachability with all the three protocols begins to decline. On the other hand, the network reachability with DAYcast is very close to that with SPB, and both DAYcast and SPB outperform WPB. 5 shows the successful delivery ratio with DAYcast, SPB, and WPB under different node density, respectively. It is seen that the delivery ratio with DAYcast is 8% and 15% larger than that with SPB and WPB, respectively. This is because DAYcast only allows the effective neighbors of a source vehicle to further broadcast a received data packet and meanwhile allows each effective neighbor to wait a different transmission delay before it broadcasts a received data packet, which effectively alleviates the effect of the broadcast storm. On the other hand, when the node density is between 30 nodes/km and 75 nodes/km, the delivery ratios with DAYcast and SPB change slightly. This is because in these cases the number of vehicles in the network is not large and the node contentions and data collisions are not serious. Both DAYcast and SPB can handle the contentions and collisions effectively. For WPB, however, the ratio decreases with the increase of the node density. With the node density further increasing, the delivery ratios with all the three solutions decrease. The reason is that in these cases the broadcast storm problem becomes more serious, which make all the three solutions unable to handle it satisfactorily. Fig. 6 shows the average delivery delay with DAYcast, SPB, and WPB under different density, respectively. It is seen that the delay with WPB changes slightly, while that with either DAYcast or SPB increases with the increase of the node density. This is because with either DAYcast or SPB a vehicle needs to wait a certain transmission delay before it broadcasts a received data packet. According to Eq. ( 1), the larger the node density is, the longer the transmission delay, causing the results. In this paper, we proposed a dynamic transmission delay based broadcast (DAYcast) protocol for VANETs. To alleviate the effect of the broadcast storm and improve the transmission efficiency of the network, DAYcast only allows the effective neighbors of a source vehicle to broadcast a received data packet and the selection of the effective neighbors are based on the position information on the onehop neighbors of the source vehicle. Meanwhile, it allows each effective neighbor to wait a certain transmission delay before it broadcasts a received packet. The transmission delay of an effective neighbor depends on the distance between the neighbor and the source vehicle, and the number of effective neighbors of the source vehicle. That is to say, DAYcast chooses the effective neighbors to reduce the number of vehicles which broadcast the received packet, and then calculates the transmission delay to avoid the contention. Thereby, DAYcast eliminates the influence of the broadcast storm to the packet transmission in a VANET. Simulation results show that DAYcast can effectively improve the network performance in terms of network reachability and the successful delivery ratio as compared with existing weighted p-persistence broadcasting (WPB) and slotted 1persistence broadcasting (SPB). In future work, we will consider a more practical road scenario with intersections and study relevant issues for the design of efficient broadcast protocols.

Fig. 1 .

 1 Fig. 1. Network model

Fig. 3 .

 3 Fig. 3.The procedure of DAYcast protocol

Fig. 4 .

 4 Fig.[START_REF] Yiannos | Speed adaptive probabilistic flooding in cooperative emergency warning[END_REF]. Network reachability Fig.5shows the successful delivery ratio with DAYcast, SPB, and WPB under different node density, respectively. It is seen that the delivery ratio with DAYcast is 8% and 15% larger than that with SPB and WPB, respectively. This is because DAYcast only allows the effective neighbors of a source vehicle to further broadcast a received data packet and meanwhile allows each effective neighbor to wait a different transmission delay before it broadcasts a received data packet, which effectively alleviates the effect of the broadcast storm. On the other hand, when the node density is between 30 nodes/km and 75 nodes/km, the delivery ratios with DAYcast and SPB change slightly. This is because in these cases the number of vehicles in the network is not large and the node contentions and data collisions are not serious. Both DAYcast and SPB can handle the contentions and collisions effectively. For WPB, however, the ratio decreases with the increase of the node density. With the node density further increasing, the delivery ratios with all the three solutions decrease. The reason is that in these cases the broadcast storm problem becomes more serious, which make all the three solutions unable to handle it satisfactorily.

Fig. 6 .

 6 Fig. 5. Successful delivery ratio

  and described as follows:

	DAYcast Protocol
	Let τ← unit time ; V i ← vehicle i; t i ← the time when a vehicle receives th e same
	packet for i times; T i ← transmission delay of vehicle i; T d ← transmission delay T a ←implicit acknowledge time; N i ← the efficient n eighbor number of vehicle i; la i ← the latitude of vehicle i; lo i ← the longitud e of vehicle i; D ij ← the distance between V i and V j; D max ← the maximum transmission range;
	Initialization: V i receives or generates a packet P from V j
	if V i receiv es P first then {
	extract N j , la j , lo j ,t 1 ;

TABLE I

 I

	. SIMULATION PARAMETERS
	Parameter	Value
	Network simulator	NS2
	Mobility simulator	VanetMobiSim
	Simulation area	2000m x 10m
	Broadcast interval	0.05s
	Transmission range	250m
	Simulation runs	10
	Average vehicle speed	50km/hr
	Simulation time	500 sec
	Number of vehicles	40, 60, 100, 150, 200, 250

This work was supported by the National Natural Science Foundation of China under Grant No. 61372105, the National Science and Technology Major Project under Grant No. 2011ZX03005-004-03, and the Research Fund of National Mobile Communications Research Laboratory, Southeast University, China, under Grant No. 2014A02.


Special Issue on Massive MIMO

Thomas L. Marzetta, Giuseppe Caire, Mérouane Debbah, Chih-Lin I, and Saif K. Mohammed Demand for wireless communications is projected to grow by more than a factor of forty or more over the next five years. A potential technology for meeting this demand is Massive MIMO (also called Large-Scale Antenna Systems, Large-Scale MIMO, ARGOS, Full-Dimension MIMO, or Hyper-MIMO), a form of multi-user multipleantenna wireless which promises orders-of-magnitude improvements in spectral-efficiency over 4G technology, and accompanying improvements in radiated energy-efficiency. The distinguishing feature of Massive MIMO is that a large number of service-antennas -possibly hundreds or even thousands -work for a significantly smaller number of active autonomous terminals. Upsetting the traditional parity between service antennas and terminals in this manner is a game-changer: The simplest multiplexing pre-coding and de-coding algorithms can be nearly optimal, expensive ultra-linear forty-Watt power amplifiers are replaced by many low-power units, and the favorable action of the law of large numbers can greatly facilitate power-control and resource-allocation.

Massive MIMO is still an emerging field. There are many unanswered theoretical questions and much remains to be done to obtain a reduction to practice. The six papers in this Special Issue are a sampling of the types of problems that are topics of active research. The papers logically fall into three categories: a) Acquisition of Channel State Information, b) Spatial Multiplexing Algorithms, and c) Massive Array Issues and Architectures.

Acquisition of Channel State Information

The acquisition of channel state information (CSI) is perhaps the central activity of any Massive MIMO system since possession of CSI by the massive array enables selective focusing of power on the down link and selective collection of power on the up link. Much of the Massive MIMO work has adopted the convenient block fading model. "Effects of Channel Aging in Massive MIMO Systems" by Truong et al. makes the more realistic assumption that CSI degrades gradually rather than abruptly: the performance of the Massive MIMO system does degrade with time, but in a graceful manner. TDD operation permits Massive MIMO to be scalable in the number of service-antennas to any desired extent. Nevertheless there remains considerable interest in making Massive MIMO work for FDD systems which predominate in Europe and North America. "A Channel State Information Feedback Method for Massive MIMO-OFDM" by Kudo et al. investigates an FDD scheme that entails down link pilots and digital CSI feedback on the up link.

Spatial Multiplexing Algorithms

Spatial multiplexing/de-multiplexing algorithms enable every active terminal to have full access to all of the time/frequency resources and are ultimately responsible for the impressive spectral efficiency of which Massive MIMO is capable. Shannon theory tells us that optimal operation on the down link is obtained by dirty-paper coding, and on the up link by joint de-coding. Massive MIMO permits the use of vastly simpler multiplexing/demultiplexing algorithms such as zero-forcing and conjugate beam-forming/matched-filtering.

Massive Array Issues and Architectures

To make Massive MIMO a practical solution it is desirable to develop scalable architectures for massive arrays studies the impact of constraints on backhaul throughputs between the digital unit and the multiplicity of radio units that comprise the massive array. This impacts both the architecture of the array as well as the functionality of the radio units.

In summary, Massive MIMO is a vital, exciting research field as evinced by the papers in this Special Issue, however much remains to be done.

  "Optimum Number of Users in Zero-Forcing Based Multiuser MIMO Systems with Large Number of Antennas" by Jung et al. utilizes random matrix theory to obtain formulas for predicting the optimum number of terminals to serve when zeroforcing is employed. Most of the Massive MIMO literature concerns concentrated arrays of service antennas. An alternative deployment for Massive MIMO is a distributed one, in which users over a large area are served simultaneously by a random distribution of access points. "Low Complexity Zero-Forcing Beamforming for Distributed Massive MIMO Systems in Large Public Venues" by Li et al. devises computationally efficient ways of performing zero-forcing for a system of distributed service antennas.

  that are based on inexpensive modular blocks. It is tempting to think that the Law of Large Numbers will compensate for the imperfections of the individual modules. "Widely-Linear Beamforming and RF Impairment Suppression in Massive Antenna Arrays", by Hakkarainen et al. studies the sensitivity of large antenna arrays to RF I/Q imbalance. Additionally, RF aware beamforming methods are proposed to reduce the unwanted effects of RF I/Q imbalance. "Before/After Precoded Massive MIMO in Cloud Radio Access Networks" by Park et al.

  Thomas L. Marzetta was born in Washington, D.C. He received the PhD in Electrical Engineering from the Massachusetts Institute of Technology in 1978. His dissertation extended, to two dimensions, the three-way equivalence of autocorrelation sequences, minimum-phase prediction error filters, and reflection coefficient sequences. He worked for Schlumberger-Doll Research (1978 -1987) to modernize geophysical signal processing for petroleum exploration. He headed a group at Nichols Research Corporation (1987 -1995) which improved automatic target recognition, radar signal processing, and video motion detection. He joined Bell Laboratories in 1995 (formerly part of AT&T, then Lucent Technologies, now Alcatel-Lucent). Within the former Mathematical Sciences Research Center, he was Director of the Communications and Statistical Sciences Department. He was an early proponent of Massive MIMO which can provide huge improvements in wireless spectral-efficiency and energy-efficiency over 4G technologies. He was a Member of the IEEE Signal Processing Society Technical Committee on Multidimensional Signal Processing, a Member of the Sensor Array and Multichannel Technical Committee, an Associate Editor for the IEEE Transactions on Signal Processing, an Associate Editor for the IEEE Transactions on Image Processing, and a Guest Associate Editor for the IEEE Transactions on Information Theory Special Issue on Signal Processing Techniques for Space-Time Coded Transmissions (Oct. 2002), for the IEEE Transactions on Information Theory Special Issue on Space-Time Transmission, Reception, Coding, and Signal Design (Oct. 2003), and for the IEEE JSAC Special Issue on Large-Scale Multiple Antenna Wireless Systems (Feb. 2013). He is currently the Lead Guest Editor for the JCN Special Issue on Massive MIMO (Aug. 2013). He was the Recipient of the 1981 ASSP Paper Award from the IEEE Signal Processing Society. He was elected a Fellow of the IEEE in Jan. 2003. He received the 2013 IEEE Guglielmo Marconi Best Paper Award. Giuseppe Caire was born in Torino, Italy, in 1965. He received the B.Sc. degree in Electrical Engineering from Politecnico di Torino (Italy), in 1990, the M.Sc. degree in Electrical Engineering from Princeton University in 1992, and the Ph.D. degree from Politecnico di Torino in 1994. He was a Recipient of the AEI G.Someda Scholarship in 1991, has been with the European Space Agency (ESTEC, Noordwijk, The Netherlands) from May 1994 to February 1995, was a Recipient of the COTRAO Scholarship in 1996 and of a CNR Scholarship in 1997. He has been visiting Princeton University in Summer 1997 and Sydney University in Summer 2000. He has been Assistant Professor in Telecommunications at the Politecnico di Torino, Associate Professor at the University of Parma, Italy, Professor with the Department of Mobile Communications at the Eurecom Institute, Sophia-Antipolis, France, and he is currently a Professor of Electrical Engineering with the Viterbi School of Engineering, University of Southern California, Los Angeles, CA. He served as Associate Editor for the IEEE Transactions on Communications in 1998-2001 and as Associate Editor for the IEEE Transactions on Information Theory in 2001-2003. He received the Jack Neubauer Best System Paper Award from the IEEE Vehicular Technology Society in 2003, and the IEEE Communications Society & Information Theory Society Joint Paper Award in 2004 and in 2011. Giuseppe Caire is Fellow of IEEE since 2005. He has served in the Board of Governors of the IEEE Information Theory Society from 2004 to 2007, and as President of the IEEE Information Theory Society in 2011. His main research interests are in the field of communications theory, information theory, channel and source coding with particular focus on wireless communications.