Emmanuelle Anceaume
email: emmanuelle.anceaume@irisa.fr

Romaric Ludinard
email: romaric.ludinard@inria.fr

Bruno Sericola
email: bruno.sericola@inria.fr

Performance Evaluation of Large-scale Dynamic Systems

In this paper we present an in-depth study of the dynamicity and robustness properties of large-scale distributed systems, and in particular of peer-to-peer systems. When designing such systems, two major issues need to be faced. First, population of these systems evolves continuously (nodes can join and leave the system as often as they wish without any central authority in charge of their control), and second, these systems being open, one needs to defend against the presence of malicious nodes that try to subvert the system. Given robust operations and adversarial strategies, we propose an analytical model of the local behavior of clusters, based on Markov chains. This local model provides an evaluation of the impact of malicious behaviors on the correctness of the system. Moreover, this local model is used to evaluate analytically the performance of the global system, allowing to characterize the global behavior of the system with respect to its dynamics and to the presence of malicious nodes and then to validate our approach.

INTRODUCTION

The adoption of peer-to-peer overlay networks as a building block for architecting Internet scale systems has firstly raised the attention of making these systems resilient to nodes-driven dynamics. This dynamics -that represents the propensity of thousands or millions of nodes to continuously join and leave the system, and which is usually called churn-if not efficiently managed, quickly gives rise to dropped messages and data inconsistency, and thus to an increasing latency and bandwidth due to repairing mechanisms. The other fundamental issue faced by any practical open system is the inevitable presence of malicious nodes [START_REF] Castro | Secure routing for structured peer-to-peer overlay networks[END_REF][START_REF] Singh | Eclipse attacks on overlay networks: Threats and defenses[END_REF][START_REF] Srivatsa | Vulnerabilities and security threats in structured peer-to-peer systems: A quantitiative analysis[END_REF]. Guaranteeing the liveness and safety of these systems requires their ability to self-heal or at least to self-protect against this adversity. Malicious nodes can devise complex strategies to subvert the system. In particular these attacks aimed at exhausting key resources of hosts (e.g., bandwidth, CPU processing, TCP connection resources) to diminish their capacity to provide or receive services [START_REF] Naoumov | Exploiting p2p systems for DDoS attacks[END_REF]. Different approaches have been proposed to face adversarial setting, each one focusing on a particular adversary strategy. Regarding eclipse attacks, a very common technique, called constrained routing table, relies on the uniqueness and impossibility of forging nodes identifiers. It consists in selecting neighbors based on their identifiers so that all of them are close to some particular points in the identifier space [START_REF] Castro | Secure routing for structured peer-to-peer overlay networks[END_REF]. Such an approach has been implemented into several overlays (e.g., [START_REF] Ratnasamy | A scalable content-addressable network[END_REF][START_REF] Stoica | Chord: A scalable peer-to-peer lookup service for internet applications[END_REF][START_REF] Druschel | Past: A large-scale, persistent peer-to-peer storage utility[END_REF]). To prevent messages from being mis-routed or dropped, the seminal works on routing security by Castro et al. [START_REF] Castro | Secure routing for structured peer-to-peer overlay networks[END_REF] and Sit and Morris [START_REF] Sit | Security considerations for peer-to-peer distributed hash tables[END_REF] combine routing failure tests and redundant routing as a solution to ensure robust routing. Their approach has then been successfully implemented in different structured-based overlays (e.g., [START_REF] Fiat | Making chord robust to Byzantine attacks[END_REF][START_REF] Anceaume | PeerCube: an hypercube-based p2p overlay robust against collusion and churn[END_REF][START_REF] Hildrum | Asymptotically efficient approaches to fault-tolerance in peer-to-peer networks[END_REF]). However all these solutions assume that malicious nodes are uniformly distributed in the system. Actually, malicious nodes may concentrate their power to attack some region of the overlay [START_REF] Awerbuch | Towards scalable and robust overay networks[END_REF][START_REF] Anceaume | Performance analysis of large scale peer-to-peer overlays using markov chains[END_REF]. These targeted attacks require additional mechanisms to be tolerated or at least to be confined. Actually, it has been shown [START_REF] Awerbuch | Group spreading: A protocol for provably secure distributed name service[END_REF] that structured overlays cannot survive targeted attacks if the adversary may keep sufficiently long its malicious peers at the same position in the overlay. Indeed, once malicious peers have succeeded in sitting in a focused region of the overlay, they can progressively gain the quorum within this region by simply waiting for honest peers to leave their position, leading to the progressive isolation of honest peers. The two fundamental properties that prevent peers isolation are the guarantee that the distribution of peers identifiers is random, and that peers cannot stay forever at the same position in the system [START_REF] Awerbuch | Towards scalable and robust overay networks[END_REF].

In the present work, we present an in-depth study of the dynamic and robustness aspects of large-scale distributed systems. We investigate adversarial strategies that aim at isolating honest nodes in cluster-based overlays. We propose an analytical model of the local behavior of clusters, based on Markov chains. This local model provides an evaluation of the power of malicious behaviors. We use this local model to evaluate analytically the performance of the global system, allowing to characterize the global behavior of the system with respect to its dynamics and to the presence of malicious nodes and then to validate our approach. Our analysis shows that i) by gathering nodes in clusters, ii) by preventing nodes to stay infinitely long at the same position in the overlay and iii) by introducing randomness in the operations of the overlay, the impact of the adversary on the cluster correctness is very limited. In particular we show that the adversary has no incentive to trigger topological operations on the system which confines its impact on particular regions of the system.

The current work extends the study done in [START_REF] Anceaume | Performance analysis of large scale peer-to-peer overlays using markov chains[END_REF], where the focus was placed on the local behavior of a cluster rather than on the global behavior of the system as done in the present work. Specifically, [START_REF] Anceaume | Performance analysis of large scale peer-to-peer overlays using markov chains[END_REF] studies according to different levels of adversity in the cluster, the impact of randomized operations on the time spent by a cluster in safe and corrupted states, and on the proportion of clusters that split or 1 merge in a safe or corrupted state. The main lessons drawn from that study is that a small amount of randomization is sufficient to defend against malicious nodes, even when they collude together to increase their power. This result is interesting because it notably decreases the complexity of the overlay operations.

The remainder of this paper is as follows: In Section 2, we briefly present the main features of structured based overlays. In Section 3 we present an analytical model of the local behavior of a cluster based on Markov chains. Section 4 uses this model to evaluate the performance of the global overlay. This analysis is applied to a large scale dynamic system presented in Section 5. Section 6 concludes.

CLUSTER-BASED SYSTEMS

Self-organization of Nodes

We consider a dynamic system populated by a large collection of nodes in which each node is assigned a unique and permanent random identifier from an m-bit identifier space. Node identifiers (simply denoted ids in the following) are derived by applying some standard strong cryptographic hash function on nodes intrinsic characteristics. The value of m (160 for the standard SHA1 hash function) is chosen to be large enough to make the probability of identifiers collision negligible. Each application-specific object, or data-item, of the system is assigned a unique identifier, called key, selected from the same m-bit identifier space. Each node p owns a fraction of all the data items of the system. The system is subject to churn, which is classically defined as the rate of turnover of nodes in the system [START_REF] Godfrey | Minimizing churn in distributed systems[END_REF]. For scalability reasons, each node locally knows only a small set of nodes existing within the system. This set is typically called the node's local view or the node's neighborhood. The particular algorithm used by nodes to build their local view and to route messages induces the resulting overlay topology. Structured overlays (also called Distributed Hash Tables (DHTs)) build their topology according to structured graphs (e.g., hypercube, torus). Specifically, nodes self-organize within the structured graph according to a distance function D based on their ids, plus possibly other criteria such as geographical distance. Finally, the mapping between nodes and data derives from the same distance function D.

In cluster-based overlays, clusters of nodes substitute for nodes at the vertices of the structured graph. Each vertex of the structured graph is composed of a set or cluster of nodes. Clusters in the overlay are uniquely identified. Cluster size s evolves according to join and leave events. However for robustness and scalability reasons, s is both lower and upper bounded. The lower bound Smin satisfies some constraint based on the assumed failure model, while the upper bound Smax is typically in O(log N ) where N is the current number of nodes in the overlay. Nodes join the clusters according to distance D. For instance in PeerCube [START_REF] Anceaume | PeerCube: an hypercube-based p2p overlay robust against collusion and churn[END_REF], p joins the (unique) cluster whose identifier is a prefix of p's identifier. Once a cluster size exceeds Smax, this cluster splits into two smaller clusters, each one populated with the nodes that are closer to each other according to distance D. Nodes can freely leave their cluster. Once a cluster undershoots its minimal size Smin, this cluster merges with the closest cluster in its neighborhood.

In the following we assume that join and leave events have an equal chance to occur in any cluster.

Malicious Behaviors

We assume the presence of malicious nodes (also called Byzantine nodes in the distributed computing community) that try to subvert the whole system by exhibiting undesirable behaviors [START_REF] Sit | Security considerations for peer-to-peer distributed hash tables[END_REF][START_REF] Douceur | The sybil attack[END_REF]. In our context this amount to devising strategies to take the control of clusters so that all the queries that go through these controlled clusters can be deviated toward malicious nodes (this refers to eclipse attacks [START_REF] Castro | Secure routing for structured peer-to-peer overlay networks[END_REF][START_REF] Sit | Security considerations for peer-to-peer distributed hash tables[END_REF]), and all the data objects stored at these clusters can be manipulated. They can magnify their impact by colluding and coordinating their behavior. We model these strategies through a strong adversary that controls these malicious peers. A strong adversary is an adversary allowed to deviate arbitrarily far from the protocol specification. We assume that the adversary has large but bounded resources in that it cannot control more than a fraction µ (0 < µ < 1) of malicious nodes in the whole network. Note that in the following we make a difference between the whole system and the overlay. The system U encompasses all the nodes that at some point may participate to the overlay (i.e, U contains up to 2 m nodes), while the overlay N contains at any time the subset of participating nodes (i.e, N size is equal to N ≤ 2 m ). Thus, while µ represents the assumed fraction of malicious nodes in U, the goal of the adversary is to populate the overlay with a larger fraction of malicious nodes in order to subvert its correct functioning. Finally, a node which always follows the prescribed protocols is said to be honest. Note that honest nodes cannot a priori distinguish honest nodes from malicious ones.

Operations of the Overlay

From the application point of view, two key operations are provided by the system: the join operation that enables a node to join the closest cluster in the system to itself according to distance D, and the leave operation, indicating that some node leaves its cluster. Note that other operations are also provided to the application (e.g., the lookup(k) operation that enables a node to search for key k, and the put(x) operation, which allows it to insert data x in the system), however they have no impact on the dynamic of the system and a minor one on its robustness. Thus we only concentrate on analyzing and evaluating the dynamics and robustness of the system through the join and leave operations. Now, from the topology structure point of view, two operations may result in a topology modification and thus must be taken into account to evaluate the dynamics of the system, namely the split and the merge operations. When the size of a cluster exceeds Smax, this cluster splits into two new clusters, and when the size of a cluster goes below Smin, this cluster merges with other clusters to guarantee the cluster resiliency. Note that for robustness reasons, a cluster may have to temporarily exceed its maximal size Smax before being able to split into two new clusters. This guarantees that resiliency of both new clusters is met, i.e., both clusters sizes are at least equal to Smin. For the time being, we do not need to go further into the details of these operations, however Section 5 will precise them.

MODELING THE CLUSTER BEHAVIOUR

We model the effect of join and leave events in a cluster using a homogeneous discrete-time Markov chain denoted by X = {Xn, n ≥ 0}. Markov chain X represents the evolution of the number of honest and malicious nodes in the cluster. The impact of malicious nodes on the cluster correctness is application dependent. For instance a necessary and sufficient condition to prevent agreement among a set of nodes is that strictly more than a third of the population set is malicious [START_REF] Lamport | The Byzantine generals problem[END_REF]. In this Section, the conditions under which cluster correctness holds are lumped in predicate correct. This predicate will be explicit in Section 5 where the application is described. Thus we define the state of a cluster as safe if predicate correct holds, while it is defined as attacked otherwise. Markov chain X alternates between the set of all the safe states, denoted by S, and the set of all the attacked states, denoted by A, until entering the absorbing state (cf. Figure 1). This state, denoted by a, represents the logical disappearance of a cluster from the graph. This occurs when the cluster either splits into two smaller clusters (i.e., it has reached its maximal size Smax) or merges with its closest neighbor (i.e., it has reached its minimal size Smin).

The transition matrix of X, denoted by P , is partitioned with respect to the decomposition of the state space Ω = S ∪ A ∪ a P =

P S P SA P Sa P AS P A P Aa 0 0 1
where PUV is the sub-matrix of dimension |U | × |V | containing the transitions from states of U to states of V , with U, V ∈ {S, A, a}. We simply write PU for PUU . In the same way, the initial probability distribution α is partitioned by writing α = (αS αA αa), where αa = 0 and the sub-vector αU contains the initial probabilities of states U ∈ {S, A}.

MODELING THE OVERLAY NETWORK

Notations

We consider an overlay populated with n clusters denoted by D1, . . . , Dn and subject to join and leave events. Each cluster Di implements the same join, leave, split and merge operations. We assume that join and leave events are uniformly distributed throughout the overlay. Specifically, when a join or leave event occurs in the overlay it affects cluster Di with probability pi = 1/n. Thus we consider, for n ≥ 1, n Markov chains X (1) , . . . , X (n) identical to X, i.e., with the same state space Ω, the same transition probability matrix P and the same initial probability distribution α. However these chains are not independent since, at each instant, only one Markov chain, chosen with probability 1/n, is allowed to make a transition.

We denote by N

(n)

S (m) and N (n)

A (m) the respective number of safe and polluted clusters just after the m-th join or leave event, i.e., the respective number of Markov chains that are in set S and in set A at instant m. More formally, these random variables are defined, for m ≥ 0, by

N (n) S (m) = n h=1 1 {X (h) m ∈S} and N (n) A (m) = n h=1 1 {X (h) m ∈A} ,
where notation 1 {A} represents the indicator function, which is equal to 1 if condition A is true and 0 otherwise. We denote by ✶S (resp. ✶A) the column vector of dimension |S∪A| whose ith entry is equal to 1 (resp. 0) if i ∈ S and 0 (resp. 1) if i ∈ A. We denote by T the sub-matrix of P containing the transitions between the states of S ∪ A, i.e., the matrix obtained from P by removing the row and the column corresponding to the absorbing state a. The dimension of T is thus |S ∪ A|. Finally, we define αT = (αS αA). In [START_REF] Anceaume | Performance analysis of large scale peer-to-peer overlays using markov chains[END_REF], we studied the expectations of both random variables

N (n) S (m) and N (n) A (m).
In what follows we focus on their joint distribution.

Focus on the First Topological Change

We study the evolution of the n Markov chains in the set of safe states and in the set of polluted states at time m respectively when none of the n Markov chains are absorbed at time m.

For v = 0, . . . , n, we are interested in the probability p

(n) v (m) defined by p (n) v (m) = P{N (n) S (m) = n -v, N (n) A (m) = v},
which represents the probability that v Markov chains are in polluted states at time m and none of the n Markov chains are absorbed at time m.

For every m ≥ 0 and ℓ ≥ 1, we define the set S m,ℓ as

S m,ℓ = {m = (m1, . . . , m ℓ ) ∈ ◆ ℓ | m1 + • • • + m ℓ = m}.
The probability p

(n)
v (m) is given by the following theorem where, as usual, we take as convention that an empty product is equal to 1. We introduce the notation

qS(k) = P{X k ∈ S} = αT T k ✶S and qA(k) = P{X k ∈ A} = αT T k ✶A.
Theorem 1. For every m ≥ 0 and v = 0, . . . , n, we have

p (n) v (m) (1) = n v 1 n m m∈Sm,n m! m1! • • • mn! v r=1 qA(mr) n r=v+1 qS(mr).
Proof. We introduce the set H = {1, . . . , n} and, for v = 0, . . . , n, the set H(v) of all subsets of H with v elements, i.e.

H(v) = {C ⊆ H | |C| = v}. It is easily checked that |H(v)| = n v
. With this notation, we have

p (n) v (m) = C∈H(v) P{X (r) m ∈ A ∀r ∈ C, X (r) m ∈ S ∀r ∈ H\C}.
The probability mass distribution, used to choose at each instant the Markov chain that has to do a transition, being uniform, all the probabilities in the above sum are equal. We thus have, by taking C = {1, . . . , v} and using Theorem 1 of [START_REF] Anceaume | Markov chains competing for transitions: Application to large scale distributed systems[END_REF], Theorem 2. For every v = 0, . . . , n -1, we have

p (n) v (m) = n v P{X (1) m ∈ A, ..., X (v) m ∈ A, X (v+1) m ∈ S, ..., X (n) m ∈ S} = n v 1 n m m∈Sm,n m! m1! • • • mn! v r=1 qA(mr)
p (n) v (m) = n n -v m ℓ=0 m ℓ 1 n ℓ 1 - 1 n m-ℓ qS(ℓ)p (n-1) v (m -ℓ),
and, for v = n,

p (n) n (m) = m ℓ=0 m ℓ 1 n ℓ 1 - 1 n m-ℓ qA(ℓ)p (n-1)
n-1 (m-ℓ).

Proof. Extracting in the sum (1) the index mn, renaming it ℓ and using Theorem 1 for integers n -1 and mℓ, we get, for v = 0, . . . , n -1,

p (n) v (m) = n v 1 n m m ℓ=0 m ℓ qS(ℓ) × m∈S m-ℓ,n-1 (m -ℓ)! m1! • • • mn-1! v r=1 qA(mr) n-1 r=v+1 qS(mr) = n v 1 n m m ℓ=0 m ℓ qS(ℓ) (n -1) m-ℓ n-1 v p (n-1) v (m -ℓ) = n n -v m ℓ=0 m ℓ 1 n ℓ 1 - 1 n m-ℓ qS(ℓ)p (n-1) v (m -ℓ).
The second part of the proof is obtained similarly.

Let Θn be the first instant at which one of the n Markov chains X (1) , . . . , X (n) gets absorbed. This random variable is defined by Θn = inf{m ≥ 0 | ∃r such that X (r) m = a}. The random variable Θn has been studied in [START_REF] Anceaume | Markov chains competing for transitions: Application to large scale distributed systems[END_REF]. In particular we have shown that its distribution is given for every k ≥ 0 and u = 1, . . . , n -1 by

P{Θn > k} = k ℓ=0 k ℓ u n ℓ 1 - u n k-ℓ ×P{Θu > ℓ}P{Θn-u > k -ℓ}.
By taking n = 2 m and u = 2 m-1 , we have

P{Θ2m > k} = 1 2 k k ℓ=0 k ℓ P{Θ 2 m-1 > ℓ}P{Θ 2 m-1 > k-ℓ}. (2)
Relation ( 2) is particularly interesting for very large values of n. Indeed, the complexity for the computation of the distribution and the expectation of Θn is in O(log 2 n). Note also that this relation can be split into 2 identical sums plus the central term corresponding to ℓ = k/2 when k is even.

We study here the expected number of Markov chains which are in a safe state and in an attacked state at time m respectively when none of the n Markov chains are absorbed at time m.

The quantity E(N

S (m)1 {Θn>m} ) is the expected number of Markov chains that are in a safe state at time m, for Θn > m. In the same way, we define the expected value E(N (n) A (m)1 {Θn>m} ) for attacked states. Remark that, for every v = 0, . . . , n, we have (N

(n) S (m) = v, N (n) A (m) = n-v) ⇔ (N (n) S (m) = v, Θn > m).
Using this remark, we obtain the following result.

Theorem 3. For every m ≥ 0, we have

E(N (n) S (m)1 {Θn>m} ) (3) = n m ℓ=0 m ℓ 1 n ℓ 1 - 1 n m-ℓ qS(ℓ)P{Θn-1 > m -ℓ}. E(N (n) A (m)1 {Θn>m} ) (4) = n m ℓ=0 m ℓ 1 n ℓ 1 - 1 n m-ℓ qA(ℓ)P{Θn-1 > m -ℓ}.
Proof. From the previous remark, we get

E(N (n) S (m)1 {Θn>m} ) = n v=1 vP{N (n) S (m) = v, Θn > m} = n v=1 vp (n) n-v (m) = n-1 v=0 (n -v)p (n) v (m).
Using the first relation of Theorem 2, this leads to

E(N (n) S (m)1 {Θn>m} ) = n m ℓ=0 m ℓ 1 n ℓ 1 - 1 n m-ℓ qS(ℓ) n-1 v=0 p (n-1) v (m -ℓ). Since n-1 v=0 p (n-1) v (m -ℓ) = n-1 v=0 P{Θn-1 > m -ℓ, N (n) S (m -ℓ) = n -1 -v} = P{Θn-1 > m -ℓ},
we get the desired result. Relation (4) is obtained using the same lines.

It has been shown in [START_REF] Anceaume | Markov chains competing for transitions: Application to large scale distributed systems[END_REF] that for all k ≥ 0 lim

n-→∞ P{Θn > k} = (αT T ✶) k
where ✶ is the column vector containing only ones with appropriate dimension determined by the context. Using this result and Relations (3) and ( 4) we get lim

n-→∞ E(N (n) S (m)1 {Θn>m} ) n = qS(0)(αT T ✶) m = αS✶(αT T ✶) m lim n-→∞ E(N (n) A (m)1 {Θn>m} ) n = qA(0)(αT T ✶) m = αA✶(αT T ✶) m .

APPLICATION

Internet growth in recent years has motivated an increasing research in peer-to-peer systems to implement different applications such as file sharing, large scale distribution of data, streaming, and video-on-demand. All these applications share in common large data and the necessity of durable access to them. Implementing these applications over a peer-to-peer system allows to potentially benefit from the very large storage space globally provided by the many unused or idle machines connected to the network. This has led for the last few years to the deployment of a rich number of large scale storage architectures. Typically, reliable storage amounts to replicating data sufficiently enough so that despite adversarial behaviors replicas are still reachable. Given the specificities of the P2P paradigm, the replication schema must guarantee a low message overhead (one cannot afford to store many copies of very large original data), and low bandwidth requirements (upon unpredictable join and leave of nodes it is eventually necessary to recreate new copies of the data to keep a certain number of replicas). Rateless erasure coding [START_REF] Luby | Lt codes[END_REF][START_REF] Maymounkov | Online codes[END_REF][START_REF] Shokrollahi | Raptor codes[END_REF] (also called Fountain coding) meets these constraints. Specifically, an object is divided into k equal size fragments, and recoded into a potentially unlimited number of independent check blocks. Fundamental property of erasure coding is that one may recover an initial object by collecting k ′ distinct check blocks generated by different sources, with k ′ = k(1 + ǫ) with ǫ arbitrarily small.

We propose to take advantage of cluster-based overlays to store these check blocks. Specifically, each of these pieces of data is assigned a unique identifier, called key, selected from the same m-bit identifier space as nodes one (cf. Section 2). Each piece of data o is then disseminated to the cluster Di whose label is closer to o key. To take advantage of the natural redundancy present in each cluster while efficiently handling the churn of the system, the population of each cluster is organized into two sets: the core set and the spare set. Members of the core set are primarily responsible for handling messages routing, and clusters operations. Management of the core set is such that its size is maintained to constant Smin. This constant is defined according to the assumed proportion of malicious nodes in the system. Spare members are the complement number of nodes in the cluster. Size s of the spare set is such that 0 ≤ s ≤ Smax -Smin. For scalability reasons, Smax is in O(log N ) where N is the current number of nodes in the overlay. In contrast to core members, spare members are not involved in any of the overlay operations, however they are in charge of data storage. To handle malicious behaviors, robust join, leave, split and merge operations are designed as detailed in Section 5.1.

In addition to robust operations, we propose to limit the lifetime of node identifiers and to randomize their compu-tation. This is achieved by adding an incarnation number t0 to the fields that appear in the node certificate (certificates are acquired at trustworthy Certification Authorities (CAs)), and by hashing this certificate to generate the initial node identifier id 0 . The incarnation number limits the lifetime of identifiers and thus the position of the nodes in the overlay. The current incarnation k of any node p can be computed as k = ⌈(t -t0)/ℓ⌉, where t is the current time, and ℓ is the length of the lifetime of each node incarnation. Thus, the k th incarnation of node p expires when its local clock reads t0 + kℓ. At this time p must rejoin the system using its (k + 1) th incarnation. By the properties of hash functions, this guarantees that nodes are regularly pushed toward unpredictable and random positions of the overlay. At any time, any node can check the validity of the identifier of any other node q in the system, by simply calculating the current incarnation k of q and generating the corresponding identifier. This leads to the following property.

Property 1 (Limited Sojourn Time). Let Di be some cluster of the system and p some peer in the overlay. Then q belongs to Di at time t if and only if idq matches the label of Di according to distance D (we say that q is valid for Di).

Robust Operations

To protect the system against the presence of malicious nodes in the overlay, we propose to take advantage of nodes role separation at cluster level to design robust operations. Briefly, the join operation is designed so that brute force denial of service attacks are discouraged. The Leave operation impedes the adversary from predicting what is going to be the composition of the core set after a given sequence of join and leave events triggered by its malicious nodes. Finally, as both merge and split operations induce topological changes in the overlay, and more importantly may have an influence on the subset of the identifier space the adversary may gain control over1 , these operations have been designed so that the adversary has, in expectation, no interest to trigger them. Specifically, these four operations make up the overlay protocol and are specified as follows • join(p): when a node p joins the system, it joins the spare set of the closest cluster in the system (according to distance D). Core members of this cluster update their spare view to reflect p's insertion (note that the spare view update does not need to be tightly synchronized at all core members).

• leave(p): When a node p leaves a cluster either p belongs to the spare set or to the core set. In the former case, core members simply update their spare view to reflect p's departure, while in the latter case, the core view maintenance procedure is triggered. This procedure consists in replacing p with one node randomly chosen from the spare set.

• split(Di): when a cluster Di has reached the conditions to split into two smaller clusters Dj and D k , core sets of both Dj and D k are built. Priority is given (s, x, y)

0 < s < ∆ ( s + 1 , x , y + 1 ) p m ( s , x , y ) 1-pm 1 { s > 1 } ( s + 1 , x , y + 1 ) p m ( s + 1 , x , y ) 1-pm 1{ s = 1 } 1 { s < ∆ -1 } ( s , x , y ) 1{ s = ∆ - 1 } 1 { x > c } ( s + 1 , x , y ) 1 - p m ( s + 1 , x , y + 1 ) p m 1 { x ≤ c } p j ( s -1 , x , y ) 1 -p m s ( s , x , y ) d y ( s -1 , x , y -1 ) 1 -d y p m s 1 - p c ( s -1 , x + 1 , y -1 ) 1 {y > 0} ( s -1 , x , y ) 1{y= 0} 1 { x > c } ( s -1 , x + 1 , y -1 ) p m s ( s -1 , x , y ) 1 -pms 1{ x ≤ c } 1 -p m c ( s , x , y ) d x ( s -1 , x + 1 , y -1 )
p m s ( s -1 , x , y )

1pms

1 {x -1≤ c} ( s -1 , x -1 , y ) 1 {y = 0} ( s -1 , x , y -1 ) 1{y>0} 1{x -1> c} 1 -d x p m c p c p ℓ

Probabilities Value

Meaning of the probability µ ratio of Byzantine peers in the universe U Smax maximal size of a cluster Smin size of the core set of a cluster (Smin is a system parameter) ∆ Smax -Smin s current size of the spare set x number of malicious peers in the spare set y number of malicious peers in the core set d probability that the lifetime of a given peer identifier has not expired (per unit of time) pj (resp. p ℓ ) 1/2 join (resp. leave) event probability pc Smin/(Smin + s) probability for a peer to belong to the cluster core set pm µ probability that the joined peer is malicious pmc x/Smin probability for a core member to be malicious pms y/s probability for a spare member to be malicious to core members of Di and completion is done with randomly chosen spares in Di. This random choice is handled through a Byzantine-tolerant consensus run among core members of Di. Spares members of Dj (resp. D k ) are populated with the remaining spares members of Di that are closer to Dj than to D k (resp. closer to D k than to Dj).

• merge(Dj, D k ): when some cluster Dj has reached the conditions to merge (i.e., its spare set is empty), it merges with the closest cluster D k to Dj. The created cluster Di is composed by a core set whose members are the core set members of D k and by a spare set whose members are the union of the spare members of D k and the core set members of Dj.

Specification of the Adversarial Behavior

Based on the operations described here above, we investigate how malicious nodes could proceed to compromise clusters correctness. Correctness of any cluster Di is jeopardized as soon as Di core set is populated by more than a quorum c of malicious nodes where c = ⌊(Smin -1)/3⌋.

As a consequence of assigning an initial unique random node id to nodes and of periodically pushing them to random regions in the overlay, the strategy of the adversary to increase the global representation of malicious identifiers is a combination of the following two actions

• maximizing the number of malicious nodes that sit in the whole overlay, and

• minimizing the likelihood that any attacked cluster Di switches back to a safe state.

Note that there is a trade-off between these two strategies. Indeed, while in the former case, the adversary aims at triggering as many join operations as possible for its malicious nodes to be present into the clusters, in the latter case, the adversary has no interest to let clusters grow in such a way that these clusters will undergo a split operation (cf. Section 5.1). Indeed, the outcome of a split operation cannot increase the subset of identifiers space the adversary has gained control over-at best, it keeps it the same. Thus when an attacked cluster Di is close to split, the adversary does not trigger anymore malicious join operations that would lead malicious nodes to join Di, and prevents honest nodes from joining Di whenever the size s of the spare set satisfies s > 1. This guarantees that Di will not grow because of honest peers, while minimizing the likelihood that Di undergoes a merge operation as well. The reason is that by construction of the merge operation (cf. Section 5.1), when Di triggers a merge operation with its closest neighbor then all the members of Di are pushed to the spare set of the new created cluster which clearly deters the adversary from triggering merge operations. This leads to the following rule Rule 1 (Adversarial Join Strategy). Let Di be a cluster such that at time t its core set contains ℓ > c valid malicious nodes. Any join event issued by node q and received at Di at time t is discarded if (q is honest and s > 1) or (s = Smax -Smin -1).

To summarize, the strategy of the adversary is to maximize the whole subset of the identifiers space it has gained control over. This is achieved by first never asking its malicious peers to leave their cluster except if Property 1 does not hold, and second by having the maximal number of malicious peers join the system except if Rule 1 holds.

Modeling the Evolution of the Cluster

We can now instantiate Markov chain X = {Xn, n ≥ 0} (introduced in Section 3), with the specificities of our application context. Specifically, for n ≥ 1, the event Xn = (s, x, y) represents the state of a cluster after the n-th transition (i.e., the n-th join or leave event), so that the size of the spare set is equal to s, the number of malicious peers in the core set is equal to x and the number of malicious peers in the spare set is equal to y. The state space Ω of X is defined by

Ω = {(s, x, y) | 0 ≤ s ≤ ∆, 0 ≤ x ≤ Smin, 0 ≤ y ≤ s}, where ∆ = Smax -Smin. Predicate correct holds iff x ≤ c. The subset of safe states S is defined by S = {(s, x, y) | 0 < s < ∆, 0 ≤ x ≤ c, 0 ≤ y ≤ s}, while the subset of attacked states A is defined by A = {(s, x, y) | 0 < s < ∆, c + 1 ≤ x ≤ Smin, 0 ≤ y ≤ s}.
Computation of transition matrix P is illustrated in Figure 2. In this tree, each edge is labelled by a probability and each leaf represents the state of the cluster. This figure shows all the states that can be reached from state (s, x, y) and the corresponding transition probabilities. Transition probabilities depend on i) the type of event that occurs (join or leave event from the core or the spare set), ii) the type of nodes involved in this operation (honest or malicious), and iii) the ratio of malicious nodes already present in the core set. The probability associated with each one of these states is obtained by summing the products of the probabilities discovered along each path starting from the root to the leaf corresponding to the target state. For instance, the branch on the very right of the tree corresponds to the situation in which cluster Di is attacked and one of its malicious core member p is no more valid. By Property 1, p leaves Di, however as Di is attacked, the adversary bias the core management procedure by replacing p with another malicious peer from Di spare set. State (s -1, x, y -1) is reached.

Modeling and computation of Property 1 is as follows. Let d be the probability (per unit of time) that the limited lifetime of some peer p has not expired. Hence d is homogeneous to a frequency, and d × ∆t represents the probability that the lifetime of a given peer identifier has not expired during ∆t units of time. Then the probability that for all the peers belonging to a set of size z Property 1 holds is equal to d z .

Initial Distributions

In the experiments conducted for this work, we consider two initial probability distributions. The first one, denoted by α (1) , is such that the cluster starts from a median size and from a failure free state. Namely,

α (1) (s,x,y) = 1 if (s, x, y) = ( ∆ 2 , 0, 0) 0 otherwise.
The second probability distribution, denoted by α (2) , is such that the cluster starts from a median size and the number of malicious nodes in both the core and spare sets follow a binomial distribution. We have n=512, α = α (1) n=1024, α = α (1) n=2048, α = α (1) n=4096, α = α (1) n=256, α = α (2) n=512, α = α (2) n=1024, α = α (2) n=2048, α = α (2) n=4096, α = α (2) (a) Smin = 7, Smax = 14, and d = 30% n=512, α = α (1) n=1024, α = α (1) n=2048, α = α (1) n=4096, α = α (1) n=256, α = α (2) n=512, α = α (2) n=1024, α = α (2) n=2048, α = α (2) n=4096, α = α (2) (b) Smin = 10, Smin = 20, and d = 30% P {Θn > m} d=0%, α = α (1) d=30%, α = α (1) d=50%, α = α (1) d=90%, α = α (1) d=99.9%, α = α (1) d=0%, α = α (2) d=30%, α = α (2) d=50%, α = α (2) d=90%, α = α (2) d=99,9%, α = α (2) (c) Smin = 7, Smax = 14, and n = 4096 P {Θn > m} d=0%, α = α (1) d=30%, α = α (1) d=50%, α = α (1) d=90%, α = α (1) d=99.9%, α = α (1) d=0%, α = α (2) d=30%, α = α (2) d=50%, α = α (2) d=90%, α = α (2) d=99,9%, α = α (2) (d) Smin = 10, Smax = 20, and n = 4096

Experimental Results

In the remaining of this Section, experiments are conducted with a proportion µ = 25% of malicious nodes in the system (cf. Section 2), a core set size Smin ∈ {7, 10}, a maximal cluster size Smax ∈ {14, 20}, and a number n of clusters ∈ {256, . . . , 4096}. Regarding the sojourn time d of nodes, it varies from 0% (nodes have to leave when they are asked to do so) to 99.9% (nodes can stay almost infinitely long at the same position in the system, even if they receive a request to leave).

We first show in Figure 3 the distribution of the first instant at which one of the n Markov chains gets absorbed as a function of m, for different values of n, and starting from different initial states. In Figure 3 (1) and α (2) and the sojourn time d of nodes is set to 30%. As expected, the total number m of join and leave events that need to be triggered before the first cluster splits or merges increases with the size of the system. This is easily explained by the fact that events are uniformly distributed over the clusters and thus the probability for an event to reach a particular cluster decreases as the system size grows. The second observation concerns the very low impact of the initial probability distribution on the distribution of Θn. This result, while intuitively surprising, is justified by the fact that d is relatively small (d = 30%) which prevents malicious nodes from staying infinitely long in the same cluster. Consequently, the natural churn of the system overpasses the influence of malicious nodes. Finally, by comparing both Figures 3(a) and 3(b), we observe that a small increase in ∆ (i.e., from ∆ = 7 to ∆ = 10) significantly augments the first instant at which the first topological change occurs (i.e., a split or merge operation), which is very interesting because it provides a trade-off between the cost implied by the management of core sets and the one implied by topological operations. A deeper investigation into the influence of the initial distribution α on the distribution of Θn is shown in Figures 3(c) and 3(d). The first remark is that for α = α (1) all the curves obtained for the different values of d coincide, meaning that even if malicious nodes may take advantage of longer sojourn times in the same cluster, it takes too much time for them to successfully attack the cluster (due to the design of the robust join and leave operations), and thus they cannot prevent the cluster from splitting and/or merging in response to the natural churn of the system (cf. Section 5.1). Now for α = α (2) , the impact of d on the distribution of Θ(n) is Figure 4 depicts the expected proportion of safe clusters after the m-th transition for a large number of clusters (n = 4096), for a sojourn time d varying from 0% to 99.9%, and for both initial distributions α. Let us observe first that when α = α (1) (cf. Figures 4(a)) and 4(b)), the expected proportion of safe clusters is very high, even for very large values of d. This result combined with the previous observed ones confirms the impact of the join and leave operations on the incapability of the adversary to attack clusters before the first topological change. Indeed, additional graphs (not shown here for space reasons) show that the expected proportion of attacked clusters is very close to 0% which is validated by the relation

E(N (n) S (m)1 {Θn >m} ) n + E(N (n)
A (m)1 {Θn >m} ) n = P{Θn > m}. Now, when clusters start in a failure-prone state, i.e., α = α (2) , Figures 4(c) and 4(d) first show that the initial proportion of safe clusters is equal to 75% (which is in accordance with the initial proportion µ of malicious nodes in the system). Second, if we compare the curves got with α = α (1) and with α = α (2) in presence of a very low induced churn (i.e., d ≥ 90%), then we observe that for a given m, the expected proportion of safe clusters is larger when the system starts from a failureprone environment (i.e., α = α (2) ) than when it starts from a failure-free one. While first intriguing, this result is explained by the fact that the population of malicious nodes is initially larger and as none of the leave events requested for them give rise to the corresponding leave operations, a larger number of events need to be triggered before the first cluster splits or merges. This is a very interesting result as it says that combining robust operations with a very low amount of induced churn is sufficient to keep the clusters in a safe state despite the presence of a large proportion of malicious nodes. Finally, and similarly to Figure 3, increasing ∆ (as shown in Figures 4(b) and 4(d)) augments significantly the time before the first split or merge operation occurs.

CONCLUSION

In this paper we have investigated the power of malicious nodes in their capability to subvert large scale dynamic systems. Given an adversarial objective, here maximizing the whole subset of the identifier space malicious nodes succeed in taking control over, we have analytically evaluated the time (in terms of join and leave triggered events) it takes for the adversary to corrupt part of the system prior to the first topological event. Significance of such an analysis was to determine whether robust operations are capable to prevent pollution from propagating to new clusters, which would ineluctably lead to the progressive subversion of the system. Restricting the analysis to the occurrence of the first topological operation is justified by the fact that topological events occurs rarely (i.e., they require a lot of join and leave operations). It is also sufficient to accurately analyze the behavior of the system as by construction of the proposed robust operations, the adversary has no incentive to provoke the triggering of topological operations. Results of the analysis are extremely positive. First they show that by simply increasing the difference between the core and the spare sets, the time that elapses before the first topological operations is significantly augmented. Second they demonstrate that pushing nodes to random positions in the system, even very infrequently, splits adversarial coalitions. Finally, they show that even when the adversary has succeeded in attacking a cluster, the adversary is incentivized to keep their density low. Thus, pollution cannot propagate to safe clusters in their vicinity.

As a future work, we intend to extend the proposed analysis by considering any kind of graphs, imposing solely that the graph formed by correct nodes is connected. One of the issues that we will have to address is the construction of a uniform node sampling algorithm that guarantees that any correct node in the graph has an equal probability to appear in any routing table of any other correct node. The presence of an omniscient adversary makes the issue challenging essentially because such an algorithm cannot rely anymore on the uniform distribution of node identifiers.

1 Figure 1 :

 11 Figure 1: Aggregated view of Markov chain X.

  ), which completes the proof. Theorem 2 gives a recurrence relation to compute the probabilities p (n) v (m) with a polynomial complexity.

Figure 2 :

 2 Figure 2: Transition diagram for the computation of the transition probability matrix P .

  α

x µ x ( 1 -P

 1 µ) S min -x s y µ y (1µ) s-y for ( ∆ 2 , x, y) ∈ S ∪ A 0 otherwise. {Θn > m} n=256, α = α (1)

P

  {Θn > m} n=256, α = α(1)

Figure 3 :

 3 Figure 3: P{Θn > m} (Relation (2)) as a function of m, n, d and for the two initial distributions α(1) and α(2)

  (a), Smin = 7 and Smax = 14, which gives a number N of nodes in the system varying in [1792, 57344], while in Figure 3(b), Smin = 10 and Smax = 20 leading to N ∈ [2560, 81920]. In both Figures 3(a) and 3(b), the distribution of Θn is plotted for the two initial distributions α

  d) Smin = 10, Smax = 20, and α = α(2)

Figure 4 :S

 4 Figure 4: E(N (n) S (m)1 {Θn >m} ) n as a function of m and d. All these experiments are run with n = 4097 clusters

Indeed, a merge operation doubles the subset of the identifier space a cluster is responsible for, while a split operation divides it per two.


Romain Millot
email: r.millot@brgm.fr

Aimee Hegan

Philippe Négrel

Geothermal waters from the Taupo Volcanic Zone, New Zealand: Li, B and Sr isotopes characterization

Keywords: Geothermal waters, New Zealand, Lithium isotopes, Boron isotopes, Strontium isotopes

In this study, we report chemical and isotope data for 23 geothermal water samples collected in New Zealand within the Taupo Volcanic Zone (TVZ). We analyzed major and trace elements including Li, B and Sr and their isotopic compositions ( 7 Li,  11 B, 87 Sr/ 86 Sr) in high temperature geothermal waters collected from deep boreholes in different geothermal fields (Ohaaki, Wairakei, Mokai, Kawerau and Rotokawa geothermal systems). Lithium concentrations are high (from 4.5 to 19.9 mg/L) and lithium isotopic compositions ( 7 Li) are very homogeneous, being comprised between -0.5 and +1.4‰. In particular, it is noteworthy that, except the samples from the Kawerau geothermal field having slightly higher  7 Li values (+1.4%), the other geothermal waters have a very constant  7 Li signature around a mean value of 0‰ ± 0.6 (2, n=21). Boron concentrations are also high and relatively homogeneous for the geothermal samples, falling between 17.5 and 82.1 mg/L. Boron isotopic compositions ( 11 B) are all negative, and display a range between -6.7 and -1.9‰. These B isotope compositions are in agreement with those of the Ngawha geothermal field in New Zealand. Li and B isotope signatures are in a good agreement with a fluid signature mainly derived from water/rock interactions involving magmatic rocks with no evidence of seawater input. On the other hand, strontium concentrations are lower and more heterogeneous and fall between 2 and 165 µg/L. 87 Sr/ 86 Sr ratios range from 0.70549 to 0.70961. These Sr isotope compositions overlap those of the Rotorua geothermal field in New Zealand, also confirming that some geothermal waters (with more radiogenic strontium) have interacted with bedrocks from the metasedimentary basement. Each of these isotope systems on their own reveals important information about particular aspects of either water source or water/rock interaction processes, but, considered together, provide a more integrated understanding of the geothermal systems from the TVZ in New Zealand.

INTRODUCTION

In the present work, we report chemical and isotope data for 23 geothermal water samples from the Taupo Volcanic Zone (TVZ) in New Zealand. Chemical and isotope data were analyzed for these deep geothermal waters in order to provide further constraints on the characterization of the associated deep geothermal reservoirs. The present study aims therefore to characterize the fluids from the geothermal systems for the TVZ and, more specifically, to constrain the nature and origin of these fluids: two essential parameters for the characterization of a geothermal resource.

Major and trace elements are first investigated in order to determine the chemical signature of the TVZ geothermal samples. A multi-isotopic approach is then used to provide additional information for the characterization of these waters. Sr isotopes ( 87 Sr/ 86 Sr) are investigated in order to better constrain the signature of the reservoir the waters come from, given that Sr isotopes are a good tracer of water origin for groundwaters and geothermal waters [START_REF] Négrel Ph | Occurrence of mineral water springs in the stream channel of the Allier River (Massif Central, France): chemical and Sr isotope constraints[END_REF][START_REF] Ph | Geochemical study in a granitic area, the Margeride, France: chemical element behavior and 87Sr/86Sr constraints[END_REF], Millot et al. 2007). In addition, we also determined Li and B isotopic compositions of the TVZ geothermal waters with the objective of evaluating the utility of these isotopic tools to constrain the water/rocks interactions. Indeed, the isotopic composition of boron ( 11 B) is determined in an attempt to elucidate the source and process controlling B in geothermal waters [START_REF] Aggarwal | Precise and accurate determination of boron isotope ratios by multiple collector ICP-MS: origin of boron in the Ngawha geothermal system, New Zealand[END_REF], and the use of Li isotopic systematics ( 7 Li) is also explored, following recent papers indicating that Li isotopes seem to be an effective tracer of water/rock interactions in geothermal waters (Millot et al. 2007, Millot and Négrel 2007, Millot et al. 2010a).

One of the objectives of the present study is also to compare the different information that can be obtained using a multi-isotopic approach for the characterization of geothermal waters, more precisely, how the combination of these isotopic tracers could help in deciphering the water-rock interaction processes. However, it is important to keep in mind that each isotopic system can only provide information concerning the chemical behavior of the element, but considered together, this multi-isotope approach can provide a more integrated understanding of the TVZ geothermal system in New Zealand.

THE TAUPO VOLCANIC ZONE

Geothermal systems occur in many parts of New Zealand. The conventional geothermal resources of New Zealand are currently utilized to depths between 1 and 3 km, and temperatures up to 330°C. High temperature geothermal fields (T > 250°C) are principally located in the TVZ (Figure 1), with another high temperature field at Ngawha in Northland.

Moderate to low and very low temperature systems (T < 250°C) are more widely scattered. Some are associated with areas of young volcanism: in Northland, Hauraki Plains, and the coastal Bay of Plenty. Many hot springs, particularly in the South Island, are associated with faults and tectonic features. The TVZ is 20-80 km wide and extends from Mt Ruapehu in the south to the Okataina Volcanic Center in the north and continues 200 km offshore. The zone is flanked by thick aprons of welded pyroclastic flows that form shallow dipping plateaus.

Situated in the middle of the North Island of New Zealand, the TVZ is an area of both volcanic and geothermal activity. Rhyolitic volcanic activity in the area is thought to have commenced 1.6 Ma ago. The enhanced activity in the region is a result of an actively extending back-arc rift, due to the subduction of the Pacific plate beneath the North Island of New Zealand. This area is characterized by an extremely high heat flow. The average heat flux from the central zone of the TVZ, which contains most of the geothermal fields, is about 700 mW/m 2 [START_REF] Bibby | Geophysical evidence on the structure of the Taupo volcanic zone and its hydrothermal circulation[END_REF].

Within the TVZ, the different geothermal fields are distributed in two bands approximately 20 km apart. It was estimated by [START_REF] Bibby | Geophysical evidence on the structure of the Taupo volcanic zone and its hydrothermal circulation[END_REF], that 75% of the heat flow occurs in the eastern band and 25% in the west. Samples were collected from five geothermal fields within the TVZ: Kawerau, Rotokawa and Ohaaki from the eastern side, and Mokai and Wairakei from the western side (Figure 1).

Temperatures of over 300°C are found within the fluids of the geothermal systems within the TVZ, generally within the eastern fields e.g. Kawerau (315°C) and Rotokawa (330°C) [START_REF] Kissling | The spatial distribution of the geothermal fields in the Taupo Volcanic Zone, New Zealand[END_REF]. Maximum temperatures at Wairakei reach 265°C [START_REF] Kissling | The spatial distribution of the geothermal fields in the Taupo Volcanic Zone, New Zealand[END_REF], with Mokai (326°C). [START_REF] Giggenbach | Variations in the chemical and isotopic composition of fluids discharged from the Taupo Volcanic Zone, New Zealand[END_REF] used the variations of H 2 O, CO 2 and Cl in discharges from six geothermal systems within the TVZ to identify the existence of two distinct types of deep water supply. Falling along the eastern side of the TVZ: Kawerau, Ohaaki and Rotokawa display higher gas contents than the fields on the western side, including Mokai and Wairakei (illustrated by CO 2 /Cl ratios above 3.9 ±1.5 vs. 0.14 ±0.1). According to [START_REF] Giggenbach | Variations in the chemical and isotopic composition of fluids discharged from the Taupo Volcanic Zone, New Zealand[END_REF], the higher ratios were described as being of andesitic rock origin, with the lower ratios originating from rhyolitic material. The excess of volatiles present in fluid discharges from geothermal and volcanic systems along convergent plate boundaries are likely to be derived preferentially from the marine sedimentary fraction of subducted material [START_REF] Giggenbach | Variations in the chemical and isotopic composition of fluids discharged from the Taupo Volcanic Zone, New Zealand[END_REF]. The Mokai geothermal system lies 25 km northwest of Lake Taupo, and is thought to be 10 km 2 , based on resistivity measurements [START_REF] Kissling | The spatial distribution of the geothermal fields in the Taupo Volcanic Zone, New Zealand[END_REF]. The Rotokawa geothermal system lies 10 km northeast of Lake Taupo and has an area of 25 km 2 , based on resistivity measurements [START_REF] Hunt | Reinjection and gravity changes at Rotokawa geothermal field[END_REF]Bowyer 2007, Heise et al. 2008). A mix of Paleozoic and Mesozoic greywackes form the basement of the system, with a mix of ignimbrites and rhyolite lavas forming overlying layers, which are in turn covered by Holocene tuff (Krupp et al.1986[START_REF] Wilson | Controls on stibnite precipitation at two New Zealand geothermal power stations[END_REF].

Hydrothermally altered greywacke sandstones dominate the basements at both the Kawerau and Ohaaki geothermal fields, however characterization by [START_REF] Wood | Basement structure, lithology and permeability at Kawerau and Ohaaki geothermal fields[END_REF] has shown that, despite similar lithologies, there are very different petrological characteristics between the two. As a result, the Ohaaki greywackes are analogous, for their bulk compositions, to granite, and the Kawerau greywackes to quartz diorite. The basement at Ohaaki also differs from those of the other fields due to the occurrence of (approx. 3%) argillite in fine partings [START_REF] Wood | Basement structure, lithology and permeability at Kawerau and Ohaaki geothermal fields[END_REF].

In Table 1, the geothermal samples are listed with their origin (geothermal system), the borehole depth (m) and deep temperature (°C) estimated by chemical geothermometry.

ANALYTICAL METHODS

Major and trace elements

The samples were collected, after filtration and acidification on site, by MRP and Contact Energy, from the production pipeline using a fluid sampling separator. These geothermal sites are managed by Contact Energy Limited, the Tuaropaki Power Company, the Rotokawa Joint Venture and Mighty River Power Limited. All chemical analyses were performed in the BRGM laboratories using standard water analysis techniques such as Ion Chromatography (Cl), Inductively Coupled Plasma-Atomic Emission Spectroscopy (Li, B and Sr), Inductively Coupled Plasma-Mass Spectrometry (Ca and Mg), and Flame Emission Spectrometry (Na, K, Ca and SiO 2 ). Major species and trace elements were determined on conditioned samples, i.e., after filtration at 0.2 µm for the major anions, and after filtration at 0.2 µm and acidification with Suprapur HNO 3 acid (down to pH = 2) for the major cations and trace elements.

Accuracy and precision for major and trace elements was verified by repeated measurements of standard materials during the course of this study: namely Ion96-3 and LGC6020 for cations and anions and pure Li, B and Sr standard solutions (Merck). The accuracy of the major and trace element data is approx. ± 10%.

Lithium isotopes

Lithium isotopic compositions were measured using a Neptune Multi Collector ICP-MS (Thermo Fischer Scientific). 7 Li/ 6 Li ratios were normalized to the L-SVEC standard solution (NIST SRM 8545, Flesch et al. 1973) following the standard-sample bracketing method (see [START_REF] Millot | Accurate and high precision measurement of lithium isotopes in two reference materials by MC-ICP-MS[END_REF] for more details). The analytical protocol involves the acquisition of 15 ratios with 16 s integration time per ratio, and yields in-run precision better than 0.2‰ (2 m ). Blank values are low, (i.e. 0.2%), and 5 minutes wash time is enough to reach a stable background value.

The samples were prepared beforehand with chemical separation/purification by ion chromatography in order to produce a pure mono-elemental solution. Chemical separation of Li from the matrix was achieved before the mass analysis using a cationic resin (a single column filled with 3 mL of BioRad AG® 50W-X12 resin, 200-400 mesh) and HCl acid media (0.2N) for 30 ng of Li. Blanks for the total chemical extraction were less than 30 pg of Li, which is negligible, since it represents a 10 -3 blank/sample ratio.

Successful quantitative measurement of Li isotopic compositions requires 100% Li recovery.

The column was, therefore, frequently calibrated and repeated analysis of the L-SVEC standard processed through columns shows 100% Li recovery and no induced isotope fractionation due to the purification process.

The accuracy and reproducibility of the entire method (purification procedure + mass analysis) were tested by repeated measurement of a seawater sample (IRMM BCR-403) after separation of Li from the matrix, for which we obtained a mean value of  7 Li = +30.9‰ ± 0.3 (2, n=7) over the analysis period. This mean value is in good agreement with our longterm measurement ( 7 Li = +31.0‰ ± 0.5, 2, n=30, [START_REF] Millot | Accurate and high precision measurement of lithium isotopes in two reference materials by MC-ICP-MS[END_REF]) and with other values reported in the literature (see, for example, [START_REF] Millot | Accurate and high precision measurement of lithium isotopes in two reference materials by MC-ICP-MS[END_REF] for a compilation with  7 Li values for seawater ranging from +28.9 to 33.9‰). Consequently, based on long-term measurements of a seawater standard, we estimate the external reproducibility of our method to be around ± 0.5‰ (2).

Boron isotopes

Boron isotopic compositions were determined on a Finnigan MAT 261 solid source mass spectrometer in a dynamic mode. B isotope compositions were determined on conditioned samples (after filtration at 0.2 µm). For these samples, water volumes corresponding to a mass of 4 µg of B was processed using a two-step chemical purification through Amberlite IRA-743 selective resin. The boron aliquot sample (2 µg) was then loaded onto a Ta single filament with graphite, mannitol and Cs, and the B isotopes were determined by measuring the Cs 2 BO 2 + ion. The total boron blank is less than 10 ng. The values are given using the notation (expressed in ‰) relative to the NBS951 boric acid standard. The 11 B/ 10 B of replicate analyses of the NBS951 boric acid standard after oxygen correction was 4.05122 ± 0.00122 (2, n=27) during this period. The reproducibility of the  11 B determination is ± 0.3‰ (2) and the internal uncertainty is better than 0.3‰ (2 m ).

The accuracy and reproducibility of the whole procedure were verified by the repeated measurements of the IAEA-B1 seawater standard [START_REF] Gonfiantini | Intercomparison of boron isotope and concentration measurements. Part II: Evaluation of results[END_REF] for which the mean  11 B value obtained is +39.22‰ ± 0.32 (2, n=33), in agreement with the accepted value for seawater ( 11 B = +39.5‰, see data compilation reported by [START_REF] Aggarwal | The effect of instrumental mass bias on  11 B measurements: a comparison between thermal ionisation mass spectrometry and multiple-collector ICP-MS[END_REF]).

Strontium isotopes

Chemical purification of Sr (~3 µg) was performed using an ion-exchange column (Sr-Spec) before mass analysis according to a method adapted from [START_REF] Pin | Evaluation of a strontium specific extraction chromatographic method for isotopic analysis in geological materials[END_REF], with total blank <1 ng for the entire chemical procedure. After chemical separation, around 150 ng of Sr was loaded onto a tungsten filament with tantalum activator and analysed with a Finnigan MAT 262 multi-collector mass spectrometer. The 87 Sr/ 86 Sr ratios were normalized to an 86 Sr/ 88 Sr ratio of 0.1194. An average internal precision of ± 10 -5 (2 m ) was obtained and the reproducibility of the 87 Sr/ 86 Sr ratio measurements was tested by repeated analyses of the NBS987 standard, for which we obtained a mean value of 0.710243 ± 0.000010 (2, n=9)

during the period of analysis.

Chemical geothermometry

The concentrations of most dissolved elements in geothermal waters depend on the groundwater temperature and the weathered mineralogical assemblage [START_REF] White | Saline waters of sedimentary rocks. Fluids in Subsurface Environments -A Symposium[END_REF][START_REF] Ellis | Quantitative interpretations of chemical characteristics of hydrothermal systems[END_REF][START_REF] Truesdell | Geochemical techniques in exploration, summary of section III[END_REF][START_REF] Arnórsson | The geochemistry of geothermal waters in Iceland. III. Chemical geothermometry in geothermal investigations[END_REF][START_REF] Fouillac | Chemical geothermometry in CO 2 -rich thermal waters. Example of the French Massif central[END_REF]. Since concentrations can be controlled by temperature-dependent reactions, they could theoretically be used as geothermometers to estimate the deep temperature of the water. In the present study, deep temperature estimates were calculated based on SiO 2 concentrations following silica geothermometer based on quartz, chalcedony,  or  cristobalite and amorphous silica solubility [START_REF] Fournier | Estimation of underground temperatures from the silica content of water from hot springs and wet-steam wells[END_REF][START_REF] Helgeson | Summary and critique of the thermodynamic properties of rock-forming minerals[END_REF], Arnorsson et al. 1983). These temperature estimates agree well with downhole measurements and with previous data reported in the literature [START_REF] Hedenquist | The thermal and geochemical structure of the Broadlands-Ohaaki geothermal system[END_REF][START_REF] Christenson | Ohaaki reservoir chemistry: characteristics of an arc-type hydrothermal system in the Taupo Volcanic Zone, New Zealand[END_REF].

RESULTS AND COMMENTS

Major and trace elements

Geothermal waters are commonly characterized by examining the behavior of the major elements (Table 1). In this context, Figure 2 [START_REF] Mossadik | Les isotopes du bore, traceurs naturels dans les eaux. Mise au point de l'analyse en spectrométrie de masse à source solide et application à différents environnements[END_REF][START_REF] Williams | Boron isotope geochemistry during diagenesis Part II. Application to organic-rich sediments[END_REF][START_REF] Aggarwal | Precise and accurate determination of boron isotope ratios by multiple collector ICP-MS: origin of boron in the Ngawha geothermal system, New Zealand[END_REF], Millot et al. 2007, Millot and Négrel 2007, Millot et al. 2009, Millot et al. 2010a).

Because Na is mainly controlled by water/rock interactions, it is interesting to also investigate the relationships between Na and Li, B and Sr (Figure 3). Firstly, it can be observed (Figure 3a), that all the geothermal waters define a general positive relationship between Na and Li, which suggests that, like Na, Li is mainly controlled by water/rock interactions. Secondly and by contrast, when B concentrations are plotted against Na (Figure 3b), different trends emerge. Indeed, we can observe different positive correlations between B and Na, meaning that B is also controlled by water/rock interactions, but there is not a single general trend at the scale of the whole TVZ as observed for Li. From that graph (Figure 3b) geothermal samples can be divided into different groups of samples: thus, the geothermal waters from Additional information can also be obtained when the concentrations of trace elements are plotted as a function of the deep temperature of the water calculated by chemical geothermometry (Figure 4). Both Li and B concentrations appear to be relatively constant and independent of the temperature of the fluid. On the other hand, Sr concentrations appear to be anti-correlated with the temperature of the fluid (Figure 4b). Such a feature strongly suggests that dilution is occurring in the system by mixing of shallower and colder waters.

However, it is also likely that the lower Sr contents can also be controlled by calcite precipitation. Finally, in Figure 4d, SiO 2 concentrations and deep temperature show a strong correlation (R² = 0.84) but this is an induced correlation, due to the fact that, deep temperature estimates were calculated based on SiO 2 concentrations following silica geothermometers.

Li-B-Sr isotopes

Li isotopes are reported in Table 1 and lithium isotopic compositions ( 7 Li, ‰) are very homogeneous. For all the samples, the range of variation for 7 Li values is small (1.9‰ in total) ranging between -0.52 and +1.42‰, respectively, for sample BR9 (Ohaaki) and KA37

(Kawerau). Omitting the two samples from Kawerau, which have slightly higher  7 Li values (+1.42 and +1.38%), the other geothermal waters have a very constant  7 Li signature around a mean value of 0‰ ± 0.6 (2, n=21). This small range of variation (± 0.6‰) is almost the same as our external reproducibility of our method for Li isotopes analysis (± 0.5‰, 2, section 3.2.). This result means that lithium has either the same origin in these fluids and/or the process(es) that control Li isotope fractionation is (are) the same for all the geothermal water samples under consideration. Compared to scarce literature data for geothermal waters (Millot and Négrel 2007, Millot et al. 2007, 2009, 2010) (1992) for the Rotorua geothermal waters (0.70514-0.70791) and for the Ohaaki geothermal field site (0.70746, [START_REF] Grimes | The Broadlands-Ohaaki geothermal system, New Zealand Part 1. Strontium isotope distribution in well BrO-29[END_REF]) also located in New Zealand.

DISCUSSION

Collectively, the lithium, boron and strontium isotopes can be used to identify the different sources contributing to the Li-B-Sr isotopic signature and to determine the main processes controlling these elements and their isotopes in the geothermal waters of the TVZ.

First, Sr isotopes are investigated in the present work in order to better define the signature of the reservoir the geothermal waters came from, given that Sr isotopes are a good tracer of water origin for groundwaters and geothermal waters [START_REF] Goldstein | The Nd and Sr Isotopic systematics of river-water dissolved material: implications for the sources of Nd and Sr in seawater[END_REF][START_REF] Négrel Ph | Occurrence of mineral water springs in the stream channel of the Allier River (Massif Central, France): chemical and Sr isotope constraints[END_REF][START_REF] Ph | Geochemical study in a granitic area, the Margeride, France: chemical element behavior and 87Sr/86Sr constraints[END_REF][START_REF] Négrel Ph | Rare Earth Elements, neodymium and strontium isotopic systematics in mineral waters: evidence from the Massif Central, France[END_REF]. Second, lithium and boron isotopic compositions ( 7 Li and  11 B) are also considered in order to provide further constraints on the origin of Li and B in these geothermal waters.

As already mentioned, Sr isotopes display a large range in the geothermal waters of the TVZ (from 0.70549 to 0.70961). In addition, we also observed in the previous section that if dilution processes happen in the system (by shallow, colder waters), it could probably explain the distribution of Sr concentrations. On the other hand, the higher 87 Sr/ 86 Sr ratios in geothermal waters are close to the value of modern seawater ( 87 Sr/ 86 Sr = 0.70917, [START_REF] Dia | Seawater Sr isotope variation over the past 300 kyr and influence of global climate cycles[END_REF]). Collectively, these observations could reflect dilution of geothermal waters by seawater (cold and shallow water input with high Sr isotopic ratio). A seawater mixing calculation is displayed in Figure 6, in which 87 Sr/ 86 Sr ratios are plotted vs. Cl/Sr ratios. The geothermal end member is chosen to be sample RK5 (Rotokawa), which has both the lowest 87 Sr/ 86 Sr and the highest Cl/Sr ratios, and is, thus, representative of geothermal waters being least affected by potential seawater contribution. No significant contribution of seawater for strontium is observed in the geothermal systems of the TVZ, it means that Sr and its isotopes are mainly controlled by the signature of the bedrocks themselves.

As already mentioned in section 2., the TVZ is an area of rhyolitic volcanic activity and hydrothermal altered sediments (greywackes) form the basement of the system. On the left side of Figure 6, we show the ranges of 87 Sr/ 86 Sr ratios for the bedrocks of the system: volcanic country rocks (rhyolite, ignimbrite and breccia, 0.7049-0.7057) reported by [START_REF] Graham | Strontium isotope compositions of Rotorua geothermal waters[END_REF], basalts for the TVZ (0.7026-0.7052, [START_REF] Gamble | The geochemistry and petrogenesis of basalts from the Taupo Volcanic Zone and Kermadec Island Arc, S.W. Pacific[END_REF]) and metasedimentary basement rocks that have 87 Sr/ 86 Sr values between 0.705 for greywackes up to 0.725 for argilites [START_REF] Graham | Strontium isotope compositions of Rotorua geothermal waters[END_REF]). Therefore it is likely that the most radiogenic Sr isotope ratios are the result of a significant contribution of waters having interacted with bedrocks having more radiogenic strontium like metasedimentary basement rocks.

Finally, when 87 Sr/ 86 Sr ratios of the TVZ fluids are compared to literature data for other geothermal fluids (Figure 7a), the geothermal waters from the TVZ are similar to those of [START_REF] Graham | Strontium isotope compositions of Rotorua geothermal waters[END_REF] for the Rotorua geothermal waters in New Zealand (0.70514-0.70791), but are significantly different from those of Iceland, which have 87 Sr/ 86 Sr ratios ranging between 0.7032 and 0.7042 [START_REF] Millot | Nd) approach for geothermal fluid characterization in Iceland[END_REF]. However, all of the geothermal waters display a similar range for Cl/Sr ratios, except for those of the Rotokawa geothermal field, which have higher Cl/Sr ratios. It is very likely that the significant difference in the Sr isotopes signature between the geothermal systems from New Zealand and those from Iceland is related to the signature of the volcanic basement rocks. New Zealand volcanic activity commenced 1.6 Ma ago, allowing in-growth of radiogenic Sr from Rb, whereas the geothermal systems in Iceland are located in the central volcanic area, which is only < 0.8 Ma old.

Boron isotopes ( 11 B) range from -6.70‰ to -1.92‰ in the geothermal waters from the TVZ, these values are in a good agreement with a volcanic origin of the waters (< 0‰: [START_REF] Barth | Boron isotope variations in nature: a synthesis[END_REF][START_REF] Barth | Geochemical and boron, oxygen and hydrogen isotopic constraints on the origin of salinity in groundwaters from the crystalline basement of the Alpine Foreland[END_REF]. In addition, according to [START_REF] Aggarwal | Precise and accurate determination of boron isotope ratios by multiple collector ICP-MS: origin of boron in the Ngawha geothermal system, New Zealand[END_REF] for the Ngawha geothermal field, the relatively low  11 B values for the fluids implies no significant marine input into the geothermal reservoir and this is also in accord with other geochemical data, e.g. average Cl/B = 53 compared to seawater Cl/B = 4839. In Figure 7b, it is important to note that the worldwide geothermal waters have a  11 B signature almost entirely comprised between -10 and 0‰, with the exception of the samples from Reykjanes and Svartsengi geothermal fields in Iceland, for which the contribution of seawater is significant [START_REF] Millot | Nd) approach for geothermal fluid characterization in Iceland[END_REF]. And, like the Na concentration (Figure 3), the Cl/B is also variable according to the location of the geothermal field.

Lithium isotopic compositions ( 7 Li, ‰) are very homogeneous for the geothermal waters of the TVZ. In Figure 7c, we can compare Li isotope data for TVZ geothermal waters with geothermal waters from Iceland [START_REF] Millot | Nd) approach for geothermal fluid characterization in Iceland[END_REF]) and with geothermal systems from the Guadeloupe and Martinique islands (volcanic islands belonging to the Lesser Antilles arc, French West Indies, Millot et al. 2010a). The geothermal waters from the TVZ are distinct from those of other geothermal systems. Indeed, the geothermal waters from the TVZ display the lowest  7 Li yet reported for geothermal waters, but they also have the lowest Cl/Li ratios.

In addition, whereas the geothermal fields from Iceland and French West Indies can have a significant contribution of Li from seawater [START_REF] Millot | Nd) approach for geothermal fluid characterization in Iceland[END_REF](Millot et al. , 2010b)), it is obvious that seawater has no influence on the composition of the geothermal waters from New-Zealand.

Several studies of Li-isotope behaviour in near-surface environments have shown that  7 Li values do not directly reflect the signature of the bedrock, but instead are controlled by fractionation during water/rock interactions during the formation of secondary minerals of alteration [START_REF] Huh | Lithium and its isotopes in major world rivers: implications for weathering and the oceanic budget[END_REF], 2001, 2004[START_REF] Pistiner | Lithium isotope fractionation during continental weathering processes[END_REF][START_REF] Rudnick | Extreme lithium isotopic fractionation during continental weathering revealed in saprolites from South Carolina[END_REF][START_REF] Kisakürek | Behaviour of Li isotopes during continental weathering: the Bidar laterite profile, India[END_REF], 2005[START_REF] Pogge Von Strandmann | Riverine behaviour of uranium and lithium isotopes in an actively glaciated basaltic terrain[END_REF][START_REF] Vigier | The relationship between riverine lithium isotope composition and silicate weathering rates in Iceland[END_REF][START_REF] Teng | Contrasting lithium and magnesium isotope fractionation during continental weathering[END_REF][START_REF] Lemarchand | Lithium isotope systematics in a forested granitic catchment (Strengbach, Vosges Mountains, France)[END_REF], Millot et al. 2010b). And it has been suggested that the 7 Li signature in the liquid might be controlled by the preferential retention of light lithium isotope ( 6 Li) into secondary mineral phases during the weathering processes. It has also been shown that the fractionation of lithium isotopes during water/rock interactions also depends on temperature because different secondary minerals might control the uptake or release of lithium in secondary minerals depending on the temperature of interaction and the associated dissolution/precipitation reaction (Chan andEdmond 1988, Chan et al. 1992;[START_REF] Barth | Boron isotope variations in nature: a synthesis[END_REF][START_REF] Chan | Lithium isotope geochemistry of sediments and hydrothermal fluids of the Guaymas Basin, Gulf of California[END_REF][START_REF] Chan | Lithium and lithium isotope profiles through the upper oceanic crust: a study of seawater-basalt exchange at ODP Sites 504B and 896A[END_REF]. All geothermal fluids from the TVZ show homogeneous and low  7 Li values, indicating that temperature is probably not the main factor in controlling the Li isotopic composition of the geothermal fluids. Rather, it is more likely that the consistent low  7 Li values in these geothermal fluids could reflect Li leaching from the same source rock. When plotted on isotope vs. isotope diagrams (Figure 8), the following conclusions can be reached: i) Li isotope signatures are very homogeneous and do not allow the discrimination of any geothermal field, suggesting that the fluids are well-mixed for Li and that it is the process of water/rock interaction at high temperature that is the main factor that determine both Li and its isotopes distribution in geothermal waters of the TVZ; ii) B isotopes are less homogeneous, and some differences could exist between the geothermal fields;

although B and its isotopes are also mainly controlled by water/rock interaction processes, there is a small difference for  11 B values between geothermal waters from Rotokawa and Ohaaki on the one hand and those from Mokai, Wairakei and Kawerau on the other hand; iii)

the Sr isotopes signatures vary widely and confirm distinction between waters from different locations. From a general point of view, the geographical distribution of the samples within the TVZ (eastern vs. western location) shows that it could affect B and Sr but not Li isotopes, meaning that it is likely that the fluids are well-mixed for Li and less for B and Sr. In addition, for Sr isotopes, the isotopic signal also depends on the type of basement rocks.

CONCLUSIONS

The main conclusions of the Li-B-Sr isotopes characterization of the geothermal waters from the Taupo Volcanic Zone are:

• Lithium concentrations are high (ranging from 4.5 to 19.9 mg/L) and lithium isotopic compositions ( 7 Li) are homogeneous ranging between -0.5 and +1.4‰. Li isotope tracing

shows that the input of seawater is negligible in these geothermal waters and that Li and its isotopes are mainly controlled by equilibrium exchange with magmatic rocks at high temperature.

• Boron concentrations are also high (17.5 and 82.1 mg/L) and relatively homogeneous and boron isotopic compositions ( 11 B) are all negative, ranging from -6.7 to -1.9‰. This B isotope signature is in a good agreement with a fluid signature derived mainly from water/rock interactions involving magmatic rocks and no seawater input.

• Strontium concentrations (0.02 to 0.165 mg/L) are lower and more heterogeneous while 87 Sr/ 86 Sr ratios range between 0.70549 and 0.70961. These Sr isotope compositions are similar to those of local magmatic bedrocks and the highest Sr isotope ratios are the result of a significant contribution of waters having interacted with bedrocks having more radiogenic strontium like, metasedimentary basement rocks.

Each of these isotope systems on their own reveals important information about particular aspects of either water source or water/rock interaction processes, but, considered together, provide a more integrated understanding of the geothermal systems from the TVZ in New Zealand. However, the combination of Li, B and Sr isotopic systems highlights the complexity of the study of these geothermal waters, and the use of only one isotopic tool could lead to an incomplete characterization of the geothermal waters.

Figure 6

Sr isotopes plotted as a function of Cl/Sr massic ratio. Fields for volcanic country rocks, TVZ basalts and metasedimentary basement rocks have also been added, see text for comments.

Figure 7

a : Sr isotopes plotted as a function of Cl/Sr massic ratio. Comparison with geothermal waters from Rotorua geothermal system [START_REF] Graham | Strontium isotope compositions of Rotorua geothermal waters[END_REF]) and geothermal systems from Iceland [START_REF] Millot | Nd) approach for geothermal fluid characterization in Iceland[END_REF]).

b : B isotopes plotted as a function of Cl/B massic ratio. Comparison with geothermal waters from Ngawha geothermal system [START_REF] Graham | Strontium isotope compositions of Rotorua geothermal waters[END_REF], Yellowstone [START_REF] Palmer | The boron isotope systematics of the Yellowstone National Park (Wyoming) hydrothermal system: A reconnaissance[END_REF] and geothermal systems from Iceland [START_REF] Millot | Nd) approach for geothermal fluid characterization in Iceland[END_REF]).

c : Li isotopes plotted as a function of Cl/Li massic ratio. Comparison with geothermal waters from Iceland [START_REF] Millot | Nd) approach for geothermal fluid characterization in Iceland[END_REF]) and from the Antilles (French West Indies, Millot et al. 2010a).

Figure 8

Multi-isotopic (Li-B-Sr) characterization of geothermal waters from the Taupo Volcanic Zone.

  Kawerau and Rotokawa seem to define a single trend, those from Ohaaki are different, and finally, those from Mokai and Wairakei seem to also plot on a single trend. These trends correlate with the spatial distribution of geothermal fields within the TVZ. The geothermal systems from Mokai and Wairakei are located in the western side of the TVZ, and those form Kawerau, Rotokawa and Ohaaki are located in the eastern side of the TVZ. Finally, when Sr concentrations are plotted as a function of Na (Figure3c), it seems that there is no global link between these two elements, except that the geothermal waters from Kawerau and Rotokawa show a positive relationship.

Figure captions Figure 1

 1 Figure captions

Figure 2

 2 Figure 2Cl concentrations (mg/L) plotted as a function of Na concentrations (mg/L).

Figure 3

 3 Figure 3Li, B and Sr concentrations (mg/L) plotted as a function of Na concentrations (mg/L). The dashed lines represent the 2 uncertainty on the linear correlations.

Figure 4

 4 Figure 4Li, Sr, B and SiO 2 concentrations (mg/L) plotted as a function of the deep temperature (°C) estimated by geothermometry. Correlations have been added as well as a 2 interval of confidence.

Figure 5 a:

 5 Figure 5 a: Li isotopes plotted as a function of Li concentrations (mg/L). The errors bars correspond to the external reproducibility of our method for Li isotopes analysis ± 0.5‰, 2. b : B isotopes plotted as a function of B concentrations (mg/L). The errors bars correspond to the external reproducibility of our method for B isotopes analysis ± 0.3‰, 2. c : Sr isotopes plotted as a function of Sr concentrations (mg/L). The errors bars are included within the sample symbol and correspond to the external reproducibility of our method for Sr isotopes analysis ± 0.000020, 2.

  Figure 1

  illustrates the relationship between Cl, considered as a conservative element, and Na, which is largely controlled by water/rock interactions. Whereas Cl concentrations range between 666 and 2834 mg/L, Na concentrations range between 388 and 1477 mg/L. From a general point of view, we observe a good relationship between Na and Cl (Figure2), indicating that the geothermal waters display a large range of salinity, with Na and Cl concentrations increasing from Rotokawa, Kawerau, Ohaaki, and Wairakei up to the samples from Mokai.

	Other major elements range from 72 to 338 mg/L for K, from 0.005 to 0.028 mg/L for Mg, and
	from 0.1 to 27.7 mg/L for Ca. Anions concentrations range between 2.6 and 74.8 mg/L for
	SO 4 and between 1.1 and 5.9 mg/L for Br. Silica concentrations are high, ranging from 384 to
	1169 mg/L.

Concerning trace elements, and those of interest for the present work: Li concentrations are high, ranging from 4.5 to 19.9 mg/L. B concentrations are also high, ranging between 17.5 and 82.1 mg/L (BR49, Ohaaki). If we exclude the sample with the highest concentration, B contents are homogeneous, falling between 17.5 and 42.3 mg/L. Finally, Sr concentrations are clearly lower, ranging from 0.002 to 0.165 mg/L. The lithium, boron and strontium concentrations are in agreement with the data reported in the literature for worldwide geothermal waters

  , geothermal waters from the TVZ display low  7 Li values.B isotopes are also reported in Table1. The range of  11 B values is of 4.8‰ in total, from -Rotokawa and Ohaaki geothermal waters have the lowest  11 B values, and, by contrast, geothermal waters from Wairakei, Mokai and Kawerau have the highest  11 B values. Geothermal samples from this study can be compared with literature Holocene basaltic field (Kaikohe volcanic field). Measurements of water samples from the Ngawha geothermal system fall within a limited range of  11 B values between -3.9 and -3.1‰, overlapping the TVZ data.Sr isotopes are also reported in Table1.87 Sr/ 86 Sr ratios range from 0.70549 (WK245,

	6.70‰ (RK14, Rotokawa) to -1.92‰ (WK235, Wairakei). A plot of B isotopes ( 11 B, ‰) as a function of B concentrations (Figure 5b) shows that, with the exception of sample BR49 (Ohaaki) having the highest B concentrations, there is no large variation of both B isotopes concentrations in Figure 5c and show no general trend or any relationship with the spatial distribution of the samples (eastern vs. western location). However, it is noteworthy that geothermal waters from the Wairakei field display the most constant 87 Sr/ 86 Sr ratios, between and concentrations. However, Wairakei) to 0.70961 (BR49, Ohaaki). Sr isotopes are plotted as a function of Sr 0.70549 and 0.70574. Sr isotope data reported here are in agreement with those of Graham

data from the Ngawha geothermal field

[START_REF] Aggarwal | Precise and accurate determination of boron isotope ratios by multiple collector ICP-MS: origin of boron in the Ngawha geothermal system, New Zealand[END_REF]

) for B isotopes. The Ngawha geothermal field is the only high temperature geothermal field in New Zealand that is located outside the TVZ. It is located on the central axis of the Northland peninsula in a Quaternary-

Table caption Table 1

 caption1 Sample list including geothermal site, sample i.d., latitude and longitude coordinates, depth of borehole (m), deep temperature estimates given by geothermometry (°C) and date of sampling. Major and trace elements concentrations (mg/L) in the geothermal water samples are also reported in this table and Li, B and Sr isotopic compositions. d.l. is the detection limit (5 µg/L for Mg). Individual errors (2 m ) are also reported for isotopic data.

Table 1 733

 1

	87 Sr/ 86 Sr 2m		0.707855 0.000006	0.709612 0.000006	0.706058 0.000009	0.706312 0.000007	0.705900 0.000007	0.706215 0.000007	0.706837 0.000007	0.705942 0.000007	0.706079 0.000007	0.706409 0.000008	0.706572 0.000006	0.705490 0.000009	0.705600 0.000008	0.705738 0.000007	0.705587 0.000008	0.705658 0.000008	0.705616 0.000007	0.705636 0.000006	0.706789 0.000007	0.706170 0.000009	0.705919 0.000007	0.709240 0.000007	0.705554 0.000007
	 11 B 2m	‰	-3.43 0.09	-3.50 0.08	-3.93 0.06	-4.71 0.07	-4.84 0.08	-4.23 0.11	-4.35 0.09	-4.30 0.07	-4.39 0.07	-6.34 0.07	-5.45 0.08	-2.84 0.09	-3.23 0.08	-2.41 0.09	-1.92 0.09	-2.48 0.12	-2.15 0.08	-1.98 0.09	-2.36 0.09	-2.84 0.06	-1.99 0.06	-6.70 0.06	-5.25 0.07
	 7 Li 2m	‰	-0.01 0.12	0.24 0.16	-0.52 0.14	-0.30 0.24	-0.45 0.22	-0.14 0.20	-0.08 0.18	-0.17 0.24	-0.18 0.12	-0.13 0.10	0.24 0.26	0.24 0.14	0.27 0.30	0.23 0.12	0.31 0.18	0.11 0.18	0.69 0.18	0.02 0.16	0.39 0.16	1.42 0.18	1.38 0.16	0.08 0.16	0.06 0.20
	Sr	µg/l	47.9	68	99.1	55.8	47	29.7	114.7	165.2	26.8	90.4	35.6	50.3	88.3	101	120.9	99.4	85.9	28.5	50	35.1	55.3	4.3	2.1
	B	mg/L	21.71	82.13	17.50	28.22	26.31	23.06	35.10	36.91	25.53	38.40	32.06	24.81	23.39	25.97	24.06	22.44	20.75	24.79	31.10	37.34	42.28	21.80	17.90
	Li	mg/L	10.32	12.74	7.62	9.53	8.84	8.27	8.38	6.77	8.10	10.10	9.30	13.28	9.78	12.12	12.58	9.48	9.03	16.88	19.94	4.46	4.97	5.67	5.18
	SiO2	mg/L	432	967	384	709	662	792	664	554	877	984	1169	709	602	677	477	455	521	799	777	822	803	1006	1027
	Br	mg/L	2	3.7	1.9	3.1	3	2.7	3.2	3.7	2.8	4.4	4.3	5.4	4.5	5.8	5.6	4.6	4.3	5.6	5.9	2	2	1.3	1.1
	SO4	mg/L	19	4	74.8	58.8	60.5	43.5	12.4	3.9	41.3	10.2	9.5	34.1	35	34	33.7	44.4	34.9	11.7	6.6	26.4	14.6	2.6	3.9
	Cl	mg/L	728	1475	666	1080	1009	928	1173	1308	1017	1485	1508	2295	1882	2295	2229	1846	1731	2398	2834	797	932	764	674
	Ca	mg/L	0.106	2.143	1.647	1.758	1.859	1.624	2.146	12.625	1.872	2.686	3.003	10.923	20.434	27.667	20.911	21.261	20.247	7.965	12.086	1.891	2.468	1.632	1.232
	Mg	µg/L	< d.l.	8.4	8.3	5.3	15.4	< d.l.	5.1	25.4	< d.l.	8.9	27.9	< d.l.	8.6	7.5	9.2	16.7	7.2	14.6	27.2	< d.l.	6.9	< d.l.	< d.l.
	K	mg/L	94	206.5	72.2	137.4	110.2	121.4	102.2	81.1	136.1	199.2	220.5	217.4	146.1	213.2	175.6	139.3	128.2	288.6	338.2	92.9	99.4	126.7	112.1
	DATE Na	mg/L	07/10/2004 813	04/10/2005 1108	03/04/2006 739	03/04/2006 864	23/04/2007 825	01/05/2007 773	28/08/2007 839	31/08/2007 799	28/09/2007 780	23/07/2008 941	03/09/2008 892	07/09/2008 1264	08/10/2008 1080	05/03/2007 1302	04/03/2008 1250	14/04/2008 1060	30/07/2008 1007	24/02/2009 1246	25/02/2009 1477	29/01/2009 586	28/01/2009 665	08/04/2009 422	07/04/2009 388
	temperature	°C	215	280	205	249	248	261	243	226	265	280	297	251	229	251	220	207	219	300	290	270	260	315	320
	deep
	borehole depth	m	680	1270	1000	670	1000	1430	760	590	1360	1800	1880	800	600	2300	700	600	600	1679	2252	1306	1108	2500	2783
		longitude	176°18'43''	176°18'37''	176°18'24''	176°18'34''	176°18'36''	176°18'35''	176°19'06''	176°19'19''	176°17'53''	176°17'54''	176°18'17''	176°02'59''	176°04'22''	176°03'20''	176°03'23''	176°04'09''	176°04'18''	175°56'23''	175°55'34''	176°43'34''	176°43'07''	176°11'25''	176°11'40''
			E	E	E	E	E	E	E	E	E	E	E	E	E	E	E	E	E	E	E	E	E	E	E
		latitude	38°31'58''	38°32'16''	38°31'11''	38°31'10'	38°31'18''	38°31'06''	38°31'53''	38°37'06''	38°31'06''	38°31'10''	38°31'01''	38°37'04''	38°37'22''	38°36'55''	38°36'53''	38°37'32''	38°37'21''	38°31'34"	38°31'36''	38°03'52''	38°03'32''	38°36'33''	38°36'31''
			S	S	S	S	S	S	S	S	S	S	S	S	S	S	S	S	S	S	S	S	S	S	S
	Well ID		BR25	BR49	BR9	BR22	BR20	BR54	BR44	BR14	BR48	BR60	BR56	WK245	WK70	WK247	WK235	WK28	WK24	MK3	MK7	KA37	KA19	RK14	RK5
	Geothermal system		Ohaaki	Ohaaki	Ohaaki	Ohaaki	Ohaaki	Ohaaki	Ohaaki	Ohaaki	Ohaaki	Ohaaki	Ohaaki	Wairakei	Wairakei	Wairakei	Wairakei	Wairakei	Wairakei	Mokai	Mokai	Kawerau	Kawerau	Rotokawa	Rotokawa

Coordinates

734

Acknowledgements

This work was funded by the Research Division of the BRGM. This work benefited from the collaboration of BRGM Chemistry laboratories for the major and trace elemental analyses: J.P. Ghestem, T. Conte and C. Crouzet are thanked for their help, as well as M. Robert for her help in the Neptune laboratory. C. Guerrot (TIMS analytic sector, BRGM) is acknowledged for B and Sr isotope data. We cordially thank Contact Energy Limited, the Tuaropaki Power Company, the Rotokawa Joint Venture and Mighty River Power Limited for providing the samples. Ed Mroczek (GNS Science) and Tom Powell (Mighty River Power Limited) are acknowledged for critical comments. E. Petelet-Giraud is also thanked for fruitful discussions. A. Hegan was supported by a fellowship from AquaTRAIN MRTN (Contract No. MRTN-CT-2006-035420) funded by the European Commission FP6 Marie Curie Actions -Human Resources and Mobility Activity Area, Research Training Networks. RM is particularly grateful to C. Fouillac for the opportunity to work on this subject. R. Rudnick and an anonymous reviewer are acknowledged for providing helpful reviews of this manuscript. L.

Aquilina is also thanked for editorial handling and constructive comments. This is BRGM contribution n° XXXX.