Introduction

Gaussian process models (a.k.a. Kriging models) are one of the preferred choices for metamodeling nowadays [START_REF] Marrel | An efficient methodology for modeling complex computer codes with gaussian processes[END_REF], competing with a few others like polynomials, splines, generalized linear models and neural networks [START_REF] Beezer | The elements of statistical learning: Data mining, inference and prediction[END_REF][START_REF] Fang | Design and modeling for computer experiments[END_REF]. While sometimes similar in prediction accuracy to the other cited methods, Gaussian processes present among other advantages: (i) the capability to reproduce complex (and a priori unknown) nonlinear input-output relationships, (ii) the ability to interpolate the observations, and (iii) the interpretability of predictions which include an estimation of the uncertainty at each prediction point.

Gaussian process metamodels were originally developed for scalar inputs, but are now also available for functional inputs. The extension to functional inputs gives rise to a number of questions about the proper way to represent them in the metamodel: (i) which functional inputs are worth keeping as predictors, (ii) which dimension reduction method (DR) is ideal to use (e.g., B-splines, PCA, PLS), (iii) which is a suitable projection dimension, and given our choice to work with Gaussian process metamodels, also the question of (iv) which is a convenient distance to measure similarities between functional input points within the kernel function. Some of these characteristics of the model, along with the type of kernel (e.g., Gaussian, Matérn 5/2) are often chosen a priori, either based on the familiarity of the modeler with a given choice or in results of other metamodeling experiences. However, as shown in [START_REF] Betancourt | Gaussian process metamodeling of functional-input code for coastal flood hazard assessment[END_REF], the configuration of these characteristics, hereon called structural parameters of the model, has a strong impact on its prediction capability and the ideal setting might not be trivial. In this report, we introduce ACO-Gp, an Ant Colony based algorithm for the structural optimization of Gaussian process models with scalar and/or functional inputs. The architecture of ACO-Gp is tailored to the Gaussian process model, however, its principles are easily extensible to other model selection frameworks.

ACO-based structural optimization

Ant colony optimization (ACO) encompasses a large variety of optimization metaheuristics derived from the seminal work of Dorigo et al. in the early 90s [START_REF] Dorigo | Positive feedback as a search strategy[END_REF][START_REF] Dorigo | Optimization, learning and natural algorithms[END_REF]. Since then, ACO based heurstics have been proved to give remarkable results in a wide range of optimization problems, including DNA sequencing [START_REF] Blum | An ant colony optimization algorithm for dna sequencing by hybridization[END_REF], scheduling [START_REF] Li | Cloud task scheduling based on load balancing ant colony optimization[END_REF], protein-ligand docking [START_REF] Korb | An ant colony optimization approach to flexible protein-ligand docking[END_REF], assembly line balancing [START_REF] Simaria | 2-antbal: An ant colony optimisation algorithm for balancing two-sided assembly lines[END_REF] and packet-switched routing [START_REF] Zhao | An improved ant colony optimization for communication network routing problem[END_REF]. ACO has been recognized as one of the most successful research lines in the area of swarm intelligence [START_REF] Bonabeau | Swarm intelligence: from natural to artificial systems[END_REF][START_REF] Bonabeau | Inspiration for optimization from social insect behaviour[END_REF], and always seats beside evolutionary algorithms, iterated local search, simulated annealing, and tabu search among the top metaheuristic techniques [START_REF] Blum | Ant colony optimization: Introduction and recent trends[END_REF].

Biological inspiration

ACO algorithms work based on stigmergy, a mechanism for indirect inter-agent communication through traces left in the environment. Ants employ this type of communication to find efficient routes during foraging1 . The way it works is traditionally explained through a picture similar the one displayed in Figure 1. The frame sequence in alphabetical order illustrates the variant of the double bridge experiment performed by Goss et al. in the 80s [START_REF] Goss | Self-organized shortcuts in the argentine ant[END_REF]. In the experiment, the nest of a colony of ants was connected to a food source by two paths, one significantly longer than the other (frame A). At first the ants began to explore the environment by randomly distributing themselves in the two paths (frame B). Along its way, each ant left pheromone trails noticeable by its mates. As expected, the ants that took the shortest path met the food before (frame C). Most part of the ants that started first the way back to the nest, perceived the larger load of pheromones in the shortest path and went through it. The shortest path kept receiving pheromones at a incrementally higher rate than the longer one, gradually reducing the chances of an ant taking this last (frame D). After some time, the whole colony converged to towards the use of the shortest path (frame E).

Figure 1: Stigmergy used by ants to find efficient paths towards a food source. Red ants represent the ones going back from the food source to the nest.

General Ant Colony System

In ACO algorithms, a colony of artificial ants evaluate solutions to the optimization problem at hand. The quality of those solutions is informed to the colony through virtual pheromone trails which help the algorithm to converge towards a high quality solution. To this day there is more than a dozen ACO metaheuristics and probably hundreds of ACO based heuristics2 . The algorithm proposed here is inspired in the ant colony system (ACS), introduced by Dorigo and Gambardella in the late 90s [START_REF] Dorigo | Ant colony system: a cooperative learning approach to the traveling salesman problem[END_REF]. In this report we proceed directly with the explanation of our heuristic for model selection. The aforementioned reference is recommended to the reader interested in the original version of ACS.

The ACS operates over a decision network that must be defined in advance based on the structure of the solution space. To start, each ant is located on a base spot and the pheromone value of each link is initialized. Each ant generates a solution to the optimization problem by adding nodes of the decision network to its path according to a pseudo-random system of rules biased by the pheromone loads in the network. Each time an ant traverses a link, a local pheromone update takes place; the pheromone load of the link is slightly reduced in order to foment the diversification of solutions (principle of exploration). Once all ants have made a complete solution, the quality of each solution is evaluated and a global pheromone update occurs; the pheromone load of the links traversed by the best ants is increased in proportion to the quality of their corresponding solutions, striving for convergence towards high quality solutions (principle of exploitation). The process is then iterated until some stopping conditions are reached. A pseudocode for the ACS is presented in Algorithm 1.

Algorithm 1 Generic

The ACO-Gp algorithm for model selection Decision network

Considering the framework described in Section 1 for a set of ds scalar inputs and df functional inputs our optimization problem consists on making the following decisions:

• State of the i-th scalar input, from {inactive, active};

• State of the j-th functional input, from {inactive, active};

• Projection basis for the j-th functional input, from {B 1 , . . . , B n b };

• Projection dimension for the j-th functional input, from {0, . . . , q j };

• Distance for the j-th functional input, from {D 1 , . . . , D n d };

• Kernel type, from {C 1 , . . . , C nc }, with i ∈ {1, . . . ds}, j ∈ {1, . . . df } and q j the original dimension of input j. The sets {B 1 , . . . , B n b }, {D 1 , . . . , D n d } and {C 1 , . . . , C nc } correspond to the basis, distance and kernel families to be considered, in that order. The projection dimension 0 denotes no projection. In order to find a suitable combination of the parameters listed above, we let our artificial ants to move through a network with a structure similar to the one depicted in Figure 2. Such a structure prevents the constitution of senseless solutions (e.g., an input being both, inactive and active) and helps to keep the network data structures considerably simple by only defining strictly necessary links.

Transition rules

ACO is an iterative algorithm. At each iteration, a group of artificial ants are located at a base spot. Each ant builds a feasible model structure by walking from node to node, always respecting the direction of the links. At each step, an ant selects the next node based on a pseudo-stochastic mechanism defined by ( 1):

rule =    Rule 1 if q ≤ q 0 , Rule 2 otherwise, (1)
with q 0 a parameter ∈ [0, 1] and q a random value from U (0, 1). For ant r located at node a, our transition rules are defined as follows:

Rule 1. Move to the feasible neighbor node with greater pheromone load. Mathematically, it is to move to the node s specified by

s = arg max b∈Jr(a) τ ab , (2)
with J r (a) the set of feasible neighbor nodes for ant r located at node a, and τ ab the pheromone load of link (a, b).

Rule 2. Pick the next node based on a probability distribution proportional to the pheromone load of the feasible neighbor nodes. Formally, this can be expressed as moving to node b with probability

P (b) = τ ab b∈Jr(a) τ ab , ( 3
)
with J r (a) and τ ab interpreted as in rule 1.

Note that ACO-Gp does not make use of the heuristic visibility value considered in [START_REF] Dorigo | Ant colony system: a cooperative learning approach to the traveling salesman problem[END_REF] and often present in ACO algorithms. The role of this value, typically denoted η ab for the link (a, b), is to introduce a priori information about the potential benefit of including each link in the solution. For many optimization problems like routing-oriented and scheduling-oriented ones, the visibility of a link is naturally set to be a function of its inverse generalized cost (see e.g., [START_REF] Das | The minimum power broadcast problem in wireless networks: an ant colony system approach[END_REF] and [START_REF] Garmsiri | Resource leveling scheduling by an ant colony-based model[END_REF]). For optimization problems involving categorical variables (such as model selection problem at hand), this set up is often less intuitive since the order of preferences over different levels of the same factor and the relative degree of preference of each level are hard to estimate. For instance, consider the questions of: (i) what would be the order of preference over a set of 5 types of basis families, and (ii) how preferred each basis family should be. Since there is no published evidence of a general superiority of a given basis family irrespective of the regression instance at hand, it could be hard to answer these two questions. In order to include η ab in the algorithm, the same two questions would have to be answered for the distance and the kernel family, and every other feature of the regression model. Hence, we decided to remove the visibility from our model, while keeping in mind that a priori information can also be introduced in the algorithm through the initial pheromone load. This possibility is discussed at the end of Section 2.3.

Pheromone update

As explained in Section 2.2, ACS implements two pheromone update mechanisms -local and global -responsible for the diversification of solutions and the exploitation of acquired knowledge about the structure of high quality solutions. The local pheromone update is triggered each time an ant ads a note to its sequence. The pheromone load of the traverse link is slightly reduced and as a consequence, other ants are less motivated to use the same link in further decisions. In ACO-Gp, the local pheromone update operates on the link (a, b) based on the assignment

τ ab ← (1 -ρ l ) • τ ab + ρ l • τ 0 , ( 4
)
where τ ab is the current pheromone load of the link, τ 0 is its initial pheromone load, and ρ l ∈ [0, 1] is a parameter that can be interpreted as the pheromone evaporation rate.

On the other hand, global pheromone update takes place each time that a colony becomes complete, i.e., each time all ants in a colony complete a solution. This time, the pheromone load of links belonging to the best ants is increased in proportion to the quality of the corresponding solutions. For each high quality ant, the global pheromone operates on the link (a, b) based on the assignment

τ ab ← (1 -ρ g ) • τ ab + ρ g • ψ, ( 5
)
where τ ab is interpreted as in (4), ψ is a measure of the quality of the solution, and ρ g ∈ [0, 1] is a parameter that can be interpreted as the learning reinforcement rate. If multiple ants are used in the global update, (4) is applied in an iterated manner over the set of best ants.

Initial pheromone load

The ACS initiates with a base pheromone load on every link. This quantity is modified by the virtual ants during the optimization to communicate actions and learning with their peers. The initial pheromone load must be set with caution, since this value will be determinant on the proper functioning of the algorithm;

• If it is set too low, the heuristic will be prematurely and irreversibly biased towards the best solution of the first iteration, breaking down the learning capability of the system. In addition, the first iteration does not count with any learned information. Thus, the solution at which the system would get stuck might be of regular quality.

• If it is conversely set too high, the system will struggle to converge (if it manages to do so). High quality solutions will not drag the attention they deserve and ants will not be able to focus their exploration around them.

In ACO-Gp, pheromones are implicitly configured to take values exclusively in [0, 1]. Thanks to the structure of the updating queries ( 4) and ( 5), this is easily achieved by just setting the initial pheromone load of every link in the range [0, 1] and using a quality measure ψ for the solutions in that same range. For the model selection problem, one can pick for instance the Leave-one-out (LOO) cross-validated squared correlation coefficient Q 2 loocv (6), or alternatively its hold-out analogous, the predictive squared correlation coefficient Q 2 hout (see e.g., [START_REF] Nilsson | Multiway calibration in 3d qsar[END_REF]).

Q 2 loocv :=          1 - n.tr i=1 (y i -ŷi,-i ) 2 n.tr i=1 (y i -ȳ) 2          0 , (6)
with (y i ) i=1,...,n.tr the vector of observed output values, ȳ the average of that vector, ŷi,-i the LOO estimation of y i , and the operator • l defined as:

x l =    x if x ≥ l, l otherwise.
Several numerical trials allowed us to identify an initial pheromone load of 0.1, combined with populations of size 10, as a suitable configuration taking into account the setup described in the previous paragraphs. During those trials, we observed the affectation of the pheromone level and the overall behavior of the algorithm along the iterations. In addition, we were able to corroborate the drawbacks associated to excessively low or high τ 0 values, explained in the introductory paragraph of this subsection.

A special treatment for the assignment of the initial pheromone load was given to the links connecting a distance type with a projection dimension (see Figure 2). The framework presented here is general enough to account for any number and type of distance families, however, for practical purposes we adopted the norms • D,θ f and • S, θf defined in [START_REF] Betancourt | Gaussian process metamodeling of functional-input code for coastal flood hazard assessment[END_REF] as a baseline. The selection of one of these two norms not only might have a relevant impact on the predictability of the model, but also on its tractability. The norm • D,θ f requires a single length-scale parameter per functional input, indifferently of its projection dimension. In contrast, for each functional input, the norm • S, θf requires as many length-scale coefficients as projection terms. That means that the optimization of the hyperparameters of the model will almost always involve a larger number of decisions variables when using the norm • S, θf . This norm must be considered as an option as it could be the optimal choice in terms of predictability. However, if the projection dimension is allowed to take too high values, the norm • S, θf may imply substantially harder and more time consuming hyperparameters learning sessions than the norm • D,θ f . As a mechanism of regularization, we set up the initial pheromone load for the links pointing out to the norm • S, θf , based on a loss function of the form

τ f (x; k j , τ 0 , δ j , w j ) =                  τ 0 exp - |k j -1| -δ j 2σ 2 if x = 0, τ 0 if 0 < x ≤ δ j , τ 0 exp - |x -1| -δ j 2σ 2 otherwise, (7)
where k j is the original dimension of input j, x takes integer values corresponding to its possible projection dimensions, τ 0 denotes the general initial pheromone load of the heuristic, and

σ 2 = - w 2 j 2 log(.5
) .

The parameters δ j and w j draw the shape of the loss function by specifying the extension of its flat section and the smoothness of its decreasing section (see Figure 3). The first condition in ( 7) controls the case where the projection dimension is equal to zero, which as stated earlier in the document denotes no projection. For every functional input, the values of δ j and w j can be tuned by observing the normalized loss curve with sum equal to 1. This curve matches the values of the probability pie defined in transition rule 2 (see (3)), and will be used by the ants during the first iteration of the algorithm, each time that this rule is implemented. Some examples of the normalized loss function are provided in Figure 4, where we illustrate the effect of w and δ on the shape of the curve. Suitable w and δ values will vary depending on particularities of the regression task. For instance, if there are many candidate functional inputs, both values could be set relatively low in order to prevent the heuristic from building too heavy models.

As a closing remark, the initial pheromone value can be used to induce preferences on the behavior of the ants. Whenever there is some hint that one level of some feature will perform better than the others, this information can be directly placed in the initial pheromone loads of that feature. By doing so, the ants will be stimulated to test more often configurations including the expected best performing level of the feature. The advantage of specifying this information through the pheromones and not through a visibility value (see the discussion on the transition rules) is that in case the induced bias was erroneous, the ants will be able to systematically remove it through the local pheromone update. Conversely, the ants will be able to reinforce the bias through the global pheromone update if they find it fruitful.

Analytic test cases

Let us now check the performance of ACO-Gp. We set ourselves in a metamodeling framework where an expensive-to-evaluate computer code is to be substituted by a light-to-run statistical model (see e.g., [START_REF] Marrel | Calculations of sobol indices for the gaussian process metamodel[END_REF] or [START_REF] Betancourt | Gaussian process metamodeling of functional-input code for coastal flood hazard assessment[END_REF] for a wider definition of a general metamodeling problem). We consider three analytic black-box functions and we undertake the model selection problem for each of them. We proceed below with the definition of our black-box functions.

Black-box functions

Let F be the set of continuous functions from [0, 1] to R. Let x = (x 1 , x 2 ) ∈ [0, 1] 2 and f = (f 1 , f 2 ) ∈ F 2 be the vectors of scalar and continuous functional inputs in that order. Then, consider the black-box computer codes specified by the following analytic functions:

• Analytic black-box 1 G 1 : [0, 1] 2 × F 2 → R, (x, f ) → x 1 sin(x 2 ) + x 1 1 0 f 1 dt 1 -x 2 2 (max T 2 f 2 ) -(min T 2 f 2 ) . ( 8
)
• Analytic black-box 2

G 2 : [0, 1] 2 × F 2 → R, (x, f ) → x 1 sin(x 2 ) + 1 0 exp(x 1 t 1 )f 1 dt 1 -x 2 2 1 0 f 2 t 2 dt 1 . (9)
• Analytic black-box 3

G 3 : [0, 1] 2 × F 2 → R, (x, f ) → x 1 -2x 2 + 4 1 0 t 1 f 1 dt 1 + 1 0 f 2 dt 2 . ( 10
)
Note that the three functions above are totally independent and the only purpose of resolving the model selection problem for all of them is to check the robustness of our proposed algorithm to variations on the structure of the underlying input-output true model.

Data generation and heuristic setup

Here, we focus on the solution of the model selection problem and we keep the generation of synthetic input data (the experimental design) very simple. For all the three analytic functions we use the same input and output data. We generate the scalar part of the design from a grid over [0, 1] 2 . We assume that the functional inputs f 1 and f 2 are represented by vectors of size 10 and 22, respectively. We made this choice in order to include functional inputs with heterogeneous discretization in the experiment. We sampled all the values of each function randomly from U (0, 1). In total, we generated an arbitrary number of 100 input points. For each of them, we computed the corresponding output values using ( 8), ( 9) and [START_REF] Simaria | 2-antbal: An ant colony optimisation algorithm for balancing two-sided assembly lines[END_REF].

As mentioned earlier in the discussion about the initial pheromone load (Section 2.3), a series of numerical trials allowed us to identify the combination τ 0 = 0.1 with a population of size 10 to work properly for our application. During these tests, we also allowed the decision probability boundary q 0 , the evaporation rate ρ l and the learning reinforcement rate ρ g vary. More specifically, we tried values in the vicinity of (0.9, 0.1, 0.1), the values used for those parameters (in the corresponding order) by Dorigo et al. in [START_REF] Dorigo | Ant colony system: a cooperative learning approach to the traveling salesman problem[END_REF]. We found the triple (0.95, 0.1, 0.1) to offer a good trade-off between solution discovery and convergence. Finally, we fixed the regularization parameters δ and w at 2 and 1.4, respectively for both functional inputs. This setup produces the normalized loss functions displayed in Figure 5.

Results

In this section we check the performance of ACO-Gp from three different perspectives:

1) Absolute quality of the selected model;

2) Relative quality of the selected model; 3) Evolution of the solutions found throughout the iterations.

These are better described in the upcoming three subsections. An ideal general number of iterations that will work well for any model selection problem may likely not exist. Since our purpose is merely to show the performance of the heuristic, we used and arbitrary number of 20 iterations. For actual applications where the interest is to find the best possible solution the algorithm can give, we recommend instead using a stopping condition based on processing time. This way, ACO-Gp will search the best possible solution given the time budget.

Absolute model quality

The absolute quality of the model refers to its predictability, regardless of the quality of the other explored models. This dimension of quality can be assessed, for instance, by means of the Q 2 loocv or Q 2 hout statistics (see the discussion about the initial pheromone load in Section 2.3). Those are conveniently interpretable measures thanks to the fact that they get values in [0, 1] (for any worthwhile model), with a value of 1 indicating perfect fitting and 0 indicating a very poor one. Here we optimize the model structure in terms of the Q 2 loocv . Figure 6 displays the calibration plot and the Q 2 loocv for the model selected for each of the black-box functions. Based on the calibration plots, all the models are good in terms of accuracy and precision.

No evident fitting problems (e.g., skewness, heavy tails) are present in any of the models.

Relative model quality

The relative quality of the selected model is assessed by comparing it to the other explored models. This dimension of quality allows to see how special our model is. For instance, a Q 2 loocv of 0.95 is not that especial if the vast majority of models had a metrics over 0.92, and similarly, a Q 2 loocv of 0.58 is not that bad if most of the models had metrics below 0.35. In In the three cases an important portion of the explored models reported a Q 2 loocv that fell well below that of the selected model. Two important takeaways from these plots are:

• The arbitrary selection of the structural parameters is an unsafe move that may be the difference between a poorly performing model and a model of high prediction quality;

• Even if there are multiple relatively good models (as for black-boxes 1 and 2), ACO-Gp will go one step forward and deliver simply the best one.

Evolution

Finally, the analysis of the evolution of the solutions found by ACO-Gp along the iterations helps to verify that its learning mechanism is working properly. In Figure 8 we plot the Q 2 loocv of the models explored at each iteration of the heuristic during the optimization for each black-box function, along with the maximum and median Q 2 loocv values by iteration. In all cases, the maximum and median Q 2 loocv progressively improve, and the models gradually converge towards the best solution found. The sporadic drops of the median are just an effect of the randomness and the degree of exploration requested via the parameter q 0 (see [START_REF] Marrel | An efficient methodology for modeling complex computer codes with gaussian processes[END_REF]). This phenomenon is by no means a bad thing, since it prevents the algorithm from getting trapped in local optima and allows it to keep improving the best solution as we see in the plots. The optimizations presented here ran in 77.7, 30.4 and 25.8 seconds, for black-boxes 1, 2 and 3, respectively.

Conclusions

This technical report introduces the ACO-Gp Ant Colony based algorithm for structural configuration of Gaussian process regression models. After almost 30 decades of their introduction to the scientific community, the ACO algorithms remain powerful optimization tools for concurrent research problems of notable relevance. In this report, we validate the ability of our algorithm to develop high quality regression models through three analytic test cases. In all of them the results were highly satisfactory. Some tests have already been done with ACO-Gp applied to the RISCOPE data (see [START_REF] Betancourt | Gaussian process metamodeling of functional-input code for coastal flood hazard assessment[END_REF] and [START_REF] Idier | Toward a user-based, robust and fast running method for coastal flooding forecast, early warning, and risk prevention[END_REF] for more details on the data) and results seem promising as well. ACO-Gp was implemented in the funGp R package [START_REF] Betancourt | funGp: Gaussian Process Models for Scalar and Functional Inputs[END_REF] for Gaussian process regression with scalar and functional inputs, which was developed as part of the RISCOPE project. The funGp package is currently downloadable from GitHub and CRAN .

Aknowledgments

This study was conducted in the frame of the RISCOPE project, funded by the French Agence Nationale de la Recherche (ANR). We thank the ANR for this support. The application case addressed in RISCOPE motivated the development of the algorithm presented here. We are grateful to Déborah Idier and Jérémy Rohmer from BRGM for their contributions regarding the metamodeling aspects of RISCOPE.

Figure 2 :

 2 Figure 2: Prototype of the network used in ACO-Gp.

Figure 3 :

 3 Figure 3: Loss function used for regularization in dimension reduction for functional inputs using the norm • S, θf . Illustrated for an hypothetical functional input of dimension 10.

( a )

 a Effect of δ on the normalized loss function. (b) Effect of w on the normalized loss function.

Figure 4 :

 4 Figure 4: Normalized loss function for hypothetical input of dimension 10.

Figure 5 :

 5 Figure 5: Normalized loss functions used for the analytic cases.

Figure 6 :

 6 Figure 6: Calibration plot for selected models.

Figure 7

 7 Figure 7 we display the Q 2 loocv of all the models explored during the optimization for each of the black-box functions.

  (a) Black-box 1. (b) Black-box 2. (c) Black-box 3.

Figure 7 :

 7 Figure 7: Relative quality of selected models.

Figure 8 :

 8 Figure 8: Evolution of the heuristic.

ACS structure 1: while <stopping conditions remain unsatisfied> do

	2:	create a new population of ants
	3:	for <i=1:Psize> do
	4:	locate ant i at its base spot
	5:	tag ant i as partial
	6:	end for
	7:	while <there are still partial ants> do
	8:	randomly pick a partial ant
	9:	apply transition rule to select next node in its sequence
	10:	reduce pheromone load of chosen link
	11:	end while
	12:	evaluate the solution made by each ant
	13:	increase pheromone load of links in best solutions
	14:	P * ← best solution so far
	15: end while
	16: return P *

Foraging: searching for wild food resources.

Metaheuristic vs. heuristic: a metaheuristic is a generic solution technique that can be applied to a broad set of problems; an heuristic is a solution technique designed to resolve a particular problem. An heuristic might be an adaption of an metaheuristic to the particularities of the optimization task at hand.

RISCOPE research project (RISCOPE.fr ) is funded by the French Agence Nationale de la Recherche (ANR) for the period 2017-2021 (ANR, project No. 16CE04-0011). This publication reflects the views only of the author's, and the ANR cannot be considered liable for any use that may be made of the information contained therein.


Geeth Manthilake
email: geeth.manthilake@uca.fr

Federica Schiavi

Chengcheng Zhao

Mainak Mookherjee

Mohamed Ali Bouhifd

Laurent Jouffret

The electrical conductivity of liebermannite: Implications for water transport into the Earth's lower mantle

 The hopping of K + ions in liebermannite results in high electrical conductivity of more than 1 S/m.  Water can be present as both molecular H2O and hydroxyl (OH -) groups in liebermannite  High H2O contents of ocean island basalts can be related to deeply subducted hydrous liebermannite

Plain Language Summary

Liebermannite (KAlSi3O8, formally K-hollandite) is an important mineral phase in deeply subducted continental and oceanic crusts. Liebermannite exhibits high ionic conductivity and could explain low resistivity in the Earth's mantle transition zone and the upper part of the lower mantle, hinting towards the subduction of hydrated crusts into the lower mantle. The libermannite samples synthesized at these conditions indicate the presence of both molecular H2O and hydroxyl (OH -) groups in its crystal structure, suggesting that liebermannite could be an important host of H2O in the Earth's deep mantle. Compared to the basalts from mid-ocean ridges that exhibit trace quantities of water, ocean island basalts show higher water contents. It has been speculated that ocean basalts are derived from crustal components that were deeply subducted. It is also known that the lower mantle is relatively dry. Therefore, it is difficult to reconcile the higher water contents in OIBs that are derived from the deep mantle. Here, we demonstrate that the source crustal components for OIBs were likely to be hydrated. Libermannite is a key mineral in deeply subducted crusts and can efficiently host water.

Introduction

Trace quantities of H2O stored in the Earth's mantle significantly influence the physical properties of mantle rocks [START_REF] Smyth | Nominally anhydrous minerals and Earth ' s deep water cycle[END_REF]. Although H2O plays an important role by affecting both mineral properties and mantle dynamics of the upper mantle, it is well known that the Earth's lower mantle is relatively dry, as indicated by the low (< 1000 ppm wt.) H2O solubility in major lower mantle mineral phases [START_REF] Bolfan-Casanova | Water partitioning at 660 km depth and evidence for very low water solubility in magnesium silicate perovskite[END_REF][START_REF] Fu | Water Concentration in Single-Crystal (Al,Fe)-Bearing Bridgmanite Grown From the Hydrous Melt: Implications for Dehydration Melting at the Topmost Lower Mantle[END_REF][START_REF] Hirschmann | Water, Melting, and the Deep Earth H2O Cycle[END_REF]). Yet some oceanic island basalts (OIB), which are characterized by relatively high H2O contents up to ~ 2.0 wt. %, are fed by plumes rising from the lower mantle [START_REF] Deschamps | A deep mantle origin for the primitive signature of ocean island basalt[END_REF]. These plumes are thought to originate from enriched reservoirs composed of recycled crustal materials that subducted into the deep mantle [START_REF] Willbold | Trace element composition of mantle end-members: Implications for recycling of oceanic and upper and lower continental crust[END_REF]. Hydrous phases associated with subducting slabs descending into the lower mantle [START_REF] Van Der Hilst | Evidence for deep mantle circulation from global tomography[END_REF] may provide a plausible mechanism for transporting H2O into the deep mantle. Dense hydrous magnesium silicates (DHMS), such as phase H and phase D, which are known to be stable in the depleted portion of the subducting slab at lower mantle conditions, are likely to transport H2O into the lower mantle and subsequently hydrate the subducted crustal components [START_REF] Nishi | Stability of hydrous silicate at high pressures and water transport to the deep lower mantle[END_REF][START_REF] Pamato | Lower-mantle water reservoir implied by the extreme stability of a hydrous aluminosilicate[END_REF].

While the dominant component of all OIB is likely to be recycled basaltic oceanic crust, both chemical and isotopic signatures of OIB indicate the mixing of different source materials [START_REF] Hart | Heterogeneous mantle domains : signatures , genesis and mixing chronologies[END_REF]. Particularly the Enriched Mantle II (EMII)-type OIB, which have been found in the Society Islands, the Marquesas and Samoa, suggest a 5-10 % of continental and sedimentary components [START_REF] Hart | Heterogeneous mantle domains : signatures , genesis and mixing chronologies[END_REF]. Water contents up to 1.5 wt. % and low H2O/Ce (<150) have been reported in such EMII-type OIB [START_REF] Cabral | Volatile cycling of H2O, CO2, F, and Cl in the HIMU mantle: A new window provided by melt inclusions from oceanic hot spot lavas at Mangaia, Cook Islands[END_REF]. The estimates based on geochemical analyses indicate that H2O concentrations in EMII sources are about 400 ppm wt. [START_REF] Dixon | Recycled dehydrated lithosphere observed in plume-influenced mid-ocean-ridge basalt[END_REF]. The low H2O/Ce measured in EM-type magmas indicates up to 92 % H2O loss from the recycled continental components with progressive dehydration [START_REF] Dixon | Recycled dehydrated lithosphere observed in plume-influenced mid-ocean-ridge basalt[END_REF]. Despite intense dehydration, the deeply subducted slab appears to retain or rehydrate by the dehydration of dense hydrous phases such that a sufficient amount of H2O is available in the deep OIB sources [START_REF] Dixon | Recycled dehydrated lithosphere observed in plume-influenced mid-ocean-ridge basalt[END_REF].

Liebermannite, formerly known as Liebermannite, is the high-pressure polymorph of KAlSi3O8-feldspar and is likely to be stable in deeply subducted (> 300 km) continental crust and mid-oceanic ridge basalt compositions enriched in potassium [START_REF] Schmidt | Experimental Constraints on Recycling of Potassium from Subducted Oceanic Crust[END_REF]. The crystal structure of liebermannite exhibits tunnels that are formed between the assemblies of double chains of edge-sharing (Si, Al)O6 octahedral units. liebermannite (lowP) crystallizes in a tetragonal (I4/m) space group symmetry, whereas above these pressures, the tunnels are slightly squeezed, and the high-pressure polymorph of liebermannite (hiP) crystallizes in a monoclinic (I2/m) space group symmetry [START_REF] Ferroir | Equation of state and phase transition in KAlSi3O8 hollandite at high pressure[END_REF][START_REF] Nishiyama | Stability and P -V -T equation of state of KAlSi3O8 -hollandite determined by in situ Xray observations and implications for dynamics of subducted continental crust material[END_REF]. These tunnels accommodate incompatible cations such as Na, K, Sr, Ba, La, and Pb. A naturally occurring liebermannite containing Ba 2+ and K + is reported to contain up to 5 wt. % H2O in its crystal structure [START_REF] Miura | The crystal structure of hollandite[END_REF], indicating that liebermannite could be an important host for H2O.

The hiP polymorph is observed to be stable down to the Earth's lowermost mantle conditions [START_REF] Kawai | First-principles study on the high-pressure phase transition and elasticity of KAlSi3O8 hollandite[END_REF][START_REF] Sueda | A new high-pressure form of KAlSi3O8 under lower mantle conditions[END_REF]. At mantle transition zone conditions, phase equilibria studies indicate that average continental crust is composed of a mixture of liebermannite (~32 vol. %), stishovite (24 vol.%), majoritic garnet (30 vol. %) and minor amounts of calcium aluminosilicate (CAS) phases [START_REF] Irifune | Subduction of continental crust and terrigenous and pelagic sediments: an experimental study[END_REF].

The presence of tunnels in the crystal structure of liebermannite facilitates the motion of K + ions along the channel, attaining superionic conduction at high temperatures [START_REF] Furusawa | Quasielastic light scattering in the onedimensional superionic conductor hollandites[END_REF][START_REF] Khanna | Thermally Activated Microwave Conductivity in the Superionic Conductor Hollandite (K1.54Mg0.77Ti7.23O16)[END_REF][START_REF] Yoshikado | Ionic conductivity of hollandite type compounds from 100 Hz to 37.0 GHz[END_REF]. More recently, the electrical conductivity of liebermannite, relevant for the deeply subducted crust, has been examined using molecular dynamics simulations. The predicted electrical conductivity for liebermannite with a 12.5 % vacancy in the tunnel is ~20 S/m at 10 GPa and 1600 K [START_REF] He | First-principles prediction of fast migration channels of potassium ions in KAlSi3O8 hollandite: Implications for high conductivity anomalies in subduction zones[END_REF]. The magnitude of the predicted electrical conductivity resulting from ionic conduction of K + ions in the tunnels is significantly greater than that of the electrical conductivity due to extrinsic defects. Hence, owing to its high electrical conductivity, liebermannite potentially represents an ideal tracer mineral for tracking subducted pathways of continental crustal components into the Earth's lower mantle.

The unique crystal structure of liebermannite and the ability to accommodate large alkali ions and possibly H2O in its tunnels highlight its importance as a possible H2O repository in the subducted crustal components and a possible source for OIB. However, so far, studies on liebermannite have been restricted to the equation of state, elasticity [START_REF] Caracas | Elasticity of (K,Na)AlSi3O8 hollandite from lattice dynamics calculations[END_REF][START_REF] Kawai | First-principles study on the high-pressure phase transition and elasticity of KAlSi3O8 hollandite[END_REF][START_REF] Mookherjee | Detecting deeply subducted crust from the elasticity of hollandite[END_REF], and electrical conductivity [START_REF] He | First-principles prediction of fast migration channels of potassium ions in KAlSi3O8 hollandite: Implications for high conductivity anomalies in subduction zones[END_REF] of the dry variety.

In this study, we first conducted high-pressure and temperature experiments to constrain the electrical conductivity of hydrous liebermannite at the upper mantle, transition zone, and lower mantle conditions, i.e., at pressures of 12, 15, and 24 GPa respectively. Then we examined the water incorporation mechanisms in liebermannite. Finally, we explore whether the high H2O contents of OIB could have derived from hydrated crustal components of deeply subducted slabs penetrated down to the lower mantle.

Materials and Methods

The starting powder was prepared using high purity oxide mixtures of K2CO3, Al2O3, and SiO2. The mixture was slowly heated to 1000 °C overnight in a platinum crucible for decarbonate K2CO3. The decarbonated powdered mixture was then heated to 1500 °C and kept 10 minutes to obtain homogeneous glass. The time was optimized to avoid possible K loss at high temperatures.

High pressure-high temperature experiments were conducted using a 1500-ton multi-anvil apparatus. The pre-synthesis of K-feldspar for electrical conductivity measurements were conducted at 2 GPa and 1200 K in a 14/8 assmbly (Fig. 1a). To synthesize hydrous K-feldspar we add Al (OH)3 to the glass mixture to obtain ~ 100 wt. ppm H2O. The powdered mixture was packed into a capsule made of rhenium (Re), and kept ant high-pressure and high-temperaure for durations of 4 hours to obtain the cylindrical samples for electrical conductivity measurements. The water contents of K-feldspar is measured to be 257±29 wt. ppm.

The pressure and temperature conditions for the electrical conductivity measurements were carefully chosen based on the stability liebermannite [START_REF] Nishiyama | Stability and P -V -T equation of state of KAlSi3O8 -hollandite determined by in situ Xray observations and implications for dynamics of subducted continental crust material[END_REF] and represent the upper mantle (12 GPa), mantle transition zone (15 GPa), and the upper part of the lower mantle (24 GPa). Electrical conductivity measurements were performed using octahedral pressure medium composed of MgO and Cr2O3 (5 wt. %) in a 14/8 assembly (octahedron edge length/anvil truncation edge length) at 12 GPa, in a 14/6 assembly at 15 GPa, and a 10/4 assembly at 24 GPa (Fig. 1b). MgO ceramic sleeves insulated the electrode wires from the furnace. All ceramic assembly parts, including the pressure media, were baked at 1000 °C for more than 12 hours and tthen stored at 125 °C in high-vacuum furnaces (< 100 mTorr) before assembling. This step reduces the exposure of assembly components to atmospheric moisture and other impurities.The sample temperature was measured using a type-C tungsten-rhenium (W95Re5-W74Re26) thermocouple junction, placed at one end of the sample. One cable formed the thermocouple and a separate W95Re5 cable placed at the opposite side of the sample connected to the impedance spectroscopy for the electrical conductivity measurements.

The electrical conductivity of liebermannite samples was determined using the impedance spectroscopy method using the Modulab MTS Impedance gain-phase analyzer. The frequency range for the analyses was 10 6 -10 1 Hz. After compressing the assembly to the desired pressure, samples were heated to 500 K and while maintaining a temperature of 500 K, the electrical resistance of the sample was measured regularly until the sample resistance reached a stable value.

The decrease of electrical resistance observed at this step corresponds to the removal of the absorbed moisture in the sample capsule and the surrounding area, which could be incorporated into the sample at higher temperatures [START_REF] Manthilake | Electrical conductivity of lawsonite and dehydrating fluids at high pressures and temperatures[END_REF]. During the measurements, sample resistance was determined in several heating and cooling cycles until the electrical resistance of heating and cooling paths become reproducible. For the discussion, we have used the data from the last cooling cycle, which minimizes the uncertainty in electrical conductivity. The sample resistance can be obtained by fitting the impedance spectra with an appropriate equivalent circuit.

For polycrystalline samples, as we use in this study, the equivalent circuit can be illustrated by a combination of resistor-capacitor/constant phase element (R-C/CPE) components, connected in parallel, series, or a series-parallel combination. Once the electrical resistivity of the samples was determined from the fit of the impedance spectra, the electrical conductivity was determined using the radius and the axial length of the cylindrical sample measured after each experiment. The inherent assumption is that the sample geometry remained constant during the experiment.

The chemical compositions the experimental run products were investigated using a Cameca SxFiveTactis electron probe microanalyzer (EPMA) in wavelength dispersive (WDS) mode operating at an accelerating voltage of 15 kV and 4 nA beam current. The low accelerating voltage was used to avoid the K loss in liebermannite.

Powder X-ray diffraction (PXRD) analyses of experimental run products were carried out to identify the crystal structures of liebermannite and the phases present in quenched samples at 15 and 24 GPa. The PXRD was obtained for ~1 mm 3 of crushed samples that were placed on zerobackground Si sample holders. The PXRD was recorded using a Philipps X-Pert Pro diffractometer with CuKα radiation source (λ = 0.15405 nm). The PXRD patterns were recorded over the 10-90° (2θ) range in steps of 0.0167° with a counting time of 700 s per step. Lattice parameters were obtained from LeBail fits using the FULLPROF program [START_REF] Rodriguez-Carvajal | Recent advances in magnetic structure determination by neutron powder diffraction[END_REF]. The background can be obtained by linear interpolation between twenty-six points, and the peak profiles were modeled using a pseudo-Voigt function.

Raman spectroscopy was used to identify liebermannite and the minor mineral phases present in the experimental run products. Raman spectra were collected in a back-scattered geometry using an InVia confocal Raman micro-spectrometer, equipped with a 532 nm diode laser, a Peltier-cooled CCD detector, a Rayleigh rejection edge filter [START_REF] Schiavi | Water quantification in silicate glasses by Raman spectroscopy: Correcting for the effects of confocality, density and ferric iron[END_REF]. The laser power of 1-8 mW; the slit aperture of 65 μm, and a grating of 2400 l/mm were used for the present analyses. These conditions result in lateral and axial spatial resolutions of ~1 and 3 μm and a spectral resolution of 1 cm -1 . Acquisition times were 30-60 s and 60-240 s for the low-wavenumber and high-wavenumber region, respectively.

The quantification of H2O in liebermannite was performed using the Fourier Transform Infrared (FTIR) spectroscopy. Unpolarized FTIR spectra of liebermannite were acquired using a Vertex70 Bruker spectrometer equipped with a Globar light source, a KBr beamsplitter, and an MCT (Mercury-Cadmium-Tellurium alloy) detector. The beam size for the analyses was ~ 30 µm.

The spectra were obtained through a CaF2 window with a resolution of 4 cm -1 . About 2000 scans were accumulated for a single spectrum. Thin slabs were cut perpendicular to the compression axis of the experimental run products using a high-precision diamond wire saw and polished using diamond mats without acetone or bond materials, such as crystal-bond or orthodontic adhesive paste to avoid IR signal contaminations. Being the samples polycrystalline, it was not possible to perform polarized measurements on single crystals; therefore, we estimated the H2O content of each sample based on several spectra that were taken over many randomly oriented crystals. The spectra were integrated between 2600 and 3800 cm -1 , an average of several spectra was calculated for each sample and multiplied by a factor 3. The OH content was quantified using the Beer-Lambert law:

CH2O = 18.02 × A × 10 6 / ε × t × ρ
where CH2O is the concentration of hydrous species expressed in ppm H2O by weight, 18.02 is the molecular weight of H2O; A is the total integrated area of bands in the region of interest, t is the thickness of the sample in cm, ε is the molar absorption coefficient, and ρ is the density of the sample in g/L. For the estimation of H2O content in liebermannite, we used the theoretically predicted integral molar absorption coefficient of 250000 L/(mol.cm 2 ) (Koch-Müller & Rhede, 2010).

Results

The electrical conductivity of liebermannite samples at 12, 15, and 24 GPa are summarized in Fig. 2. The discontinuous increase in electrical conductivity at a temperature above 660 K during the first heating cycle of the 12 GPa experiment coincides with the transformation of K-Feldspar to liebermannite with increasing temperature (Fig. 2a). This transformation accounts for more than one order of magnitude increase of conductivity for liebermannite. After the complete transformation, the electrical conductivity of libermannite increases with increasing temperature from 700 K to 1500 K and reached to about 1 S/m at 12 GPa and 1500 K. The increase of pressure appears to have an inverse effect on the electrical conduction in liebermannite (Fig. 2b). The impedance spectra of the sample undergoing K-feldspar to libermannite transition with increasing temperature are shown (Fig. 3).

The activation enthalpy (H) of each conduction mechanism can be calculated by fitting the data to equation, 𝜎 = 𝜎 0 𝑒𝑥𝑝 (-∆𝐻|𝑅𝑇) , where σ is the electrical conductivity (S/m), T is the absolute temperature (K), σ0 is the pre-exponential factor (S/m), and R is the gas constant (J/K mol). Here, we analyzed the activation enthalpy, first by manually fitting different temperature segments to the fitting equation and, second, by fitting the equation to the entire temperature range [START_REF] Fei | The Effect of Water on Ionic Conductivity in Olivine[END_REF]. The estimated activation enthalpies for different conduction mechanisms obtaine using two fitting methods are listed in Table 1.

The EPMA, SEM, and Raman grid-analyses indicate that liebermannite is the principal phase present in samples synthesized at 12, 15, and 24 GPa (Fig. 4). However, minor amounts of wadeite, kyanite, and stishovite presence in our samples, an observation consistent with a previous study [START_REF] Urakawa | Synchrotron Radiation Study on the High-Pressure and High-Temperature Phase Relations of KAISi3O8[END_REF] (Table . 2) (Fig. 5). The grid analyses with a step size of 5 μm also confirmed the absence of melt in our samples. The minor phases mostly occur as isolated grains within the libermannite matrix (Fig 5a ) and the volume fraction of minor phases appeared to decrease with increasing pressure. Chemical analyses of minor phases at the 24 GPa sample were hampered by either absence of minor phases (e.g kyanite) at high-pressure samples or if present, their extremely small grain sizes (< 5μm).

The FTIR analyses of libermannite samples after the electrical conductivity measurements indicate water contents of 1044±87 wt. ppm at 12 GPa, 885±78 wt. ppm at 15 GPa and 1104±18

wt. ppm at 24 GPa (Fig 6). The high water contents compared to the K-feldspar starting material suggest that the absorbed moisture in the assembly parts may have been incorporated into the sample during the K-feldspar-libermannite transformation at high-temperature, yielding high water contents in the liebermannite samples (i.e. higher than the amount of water produced by the decomposition of Al(OH)3). We were not able to produce water-free liebermannite samples even with dry pressure materials.

The infrared spectra provide information on water incorporation in liebermannite (Fig. 6).

Several O-H vibrational modes are observed: (1) the stretching vibrations of weakly hydrogen bound hydroxyl units produce intense bands between 3700 and 3500 cm -1 ; these modes are more pronounced in the 12 GPa sample. (2) The stretching vibrations of both H2O molecules and more strongly hydrogen bound hydroxyl units occur in the 3425-3100 cm -1 wavenumber range; they become dominant in the 15 and 24 GPa samples (Fig. 6a). (3) The O-H stretching and bending vibrations of the hydronium (H3O + ) ion are observed near 2900 cm -1 and 1665 cm -1 ; these modes disappear at the highest pressure. In the absence of organic contaminants, such as acetone or bond materials, the coexistence of these two modes confirms the presence of H3O + in lowP liebermannite. ( 4) The H-O-H bending vibration of the water molecule slightly shifts from 1610 to 1630 cm -1 with increasing pressure (Fig. 6b). (5) The combination (stretching/bending) modes of both H2O molecules (near 5200 cm -1 ) and hydroxyl units (near 4500 cm -1 ) display continuous shift to lower frequencies with increasing pressure (Fig. 6c). This shift is compatible with the increase of the strength of hydrogen bonds with pressure, as usually observed in minerals [START_REF] Cynn | High-pressure IR spectra of lattice modes and OH vibrations in Fe-bearing wadsleyite[END_REF]. The observed spectral features suggest the presence of H2O and H3O + substituting for K + in the tunnel [START_REF] Bethell | The Infrared Spectrum of the H3O + Ion in Acid Hydrates[END_REF] and the substitution of hydrogen (H + ) in octahedral sites where it forms O-H bonds with one of the six possible oxygen ions in the vacant octahedral site.

The powder X-ray diffraction analyses were conducted on the samples synthesized at 15 GPa and 24 GPa. Both patterns can be refined with parameters consistent with the I4/m phase: a = 9.3304(4)Å, c = 2.7223(2)Å. Obtaining the same lattice parameters for the samples out of the different pressures of annealing experiments implies that the liebermannite phase is elastic enough to regain its low-pressure parameters consistent with previous studies [START_REF] Nishiyama | Stability and P -V -T equation of state of KAlSi3O8 -hollandite determined by in situ Xray observations and implications for dynamics of subducted continental crust material[END_REF][START_REF] Sueda | A new high-pressure form of KAlSi3O8 under lower mantle conditions[END_REF].

Discussion and conclusions

1 Electrical conduction in liebermannite

The electrical conductivity of K-feldspar before transforming into libermannite exhibit higher electrical conductivity compared to the previous electrical conductivity measurements of F-feldspar reported at 1 GPa [START_REF] Hu | Electrical conductivity of alkali feldspar solid solutions at high temperatures and high pressures[END_REF]. The high electrical conductivity of K-feldspar observed in our study can be explained by the 0.02 wt. % H2O contents in our samples, compared to the dry verity reported by [START_REF] Hu | Electrical conductivity of alkali feldspar solid solutions at high temperatures and high pressures[END_REF] study. Similar electrical behavior has been observed in dry albite and water-bearing albite [START_REF] Hu | Electrical conductivity of alkali feldspar solid solutions at high temperatures and high pressures[END_REF][START_REF] Ni | Electrical conductivity of dry and hydrous NaAlSi3O8 glasses and liquids at high pressures[END_REF].

The electrical conductivity of hydrous liebermannite increases discontinuously upon increasing temperature (Fig. 2b). At temperatures below 700 K, conductivity is likely to be dominated by molecular water that is loosely bound in the tunnel sites and the associated activation energy is less than 0.40 eV. In the temperature range between 700 K and 1100 K, the conductivity is likely largely due to extrinsic proton defects that are formed due to coupled substitutions of the octahedral sites. The associated activation energy varies between 0.64 eV and 0.78 eV. At temperatures above 1200 K, the conductivity is likely to be dominated by fast motions of K + or H + ions (from H3O + ) ions in the channel. The associated activation energy is greater than 0.82 eV.

Higher conductivity has also been observed in alkali bearing materials with the hollandite crystal structure [START_REF] Khanna | Thermally Activated Microwave Conductivity in the Superionic Conductor Hollandite (K1.54Mg0.77Ti7.23O16)[END_REF]. Such higher conductivity is often referred to as super-ionic conductivity [START_REF] He | First-principles prediction of fast migration channels of potassium ions in KAlSi3O8 hollandite: Implications for high conductivity anomalies in subduction zones[END_REF][START_REF] Khanna | Thermally Activated Microwave Conductivity in the Superionic Conductor Hollandite (K1.54Mg0.77Ti7.23O16)[END_REF]. In our study, we are documenting such effects at extreme pressures of 12 to 24 GPa in libermannite.

The electrical conductivity of minor phases, wadeite, stishovite, and kyanite should not interfere with the measured conductivity of liebermannite due to their low volume fractions (occur as isolated grains). Based on chemical analyses of individual liebermannite grains in our samples, we observe a maximum of 7 % K + site vacancy in our liebermannite after electrical conductivity measurements at 12, 15, and 24 GPa. The intrinsic defects seem to dominate the conduction mechanism at high temperatures, whereas extrinsic defects may be responsible for the observed conductivity at lower temperatures.

Water incorporation in liebermannite

The occurrence of several OH stretching bands (Fig. 4a) suggests that different substitution mechanisms take place, including H + substitution in octahedral sites and H2O/H3O + in tunnels, resulting in distinct hydrogen bonding strengths. Protons attachment to both bridging and nonbridging oxygen can also help to explain such a wide range of OH stretching vibrations. Strong pleochroism in the 3500-3600 cm -1 range of the polycrystalline sample at 12 GPa indicates the random orientation of the crystals, in contrast to higher pressure samples. Water/hydrogen could be incorporated in liebermannite structure in a variety of sites including the substitution of H + in the silicon octahedral sites, i.e., SiSi 4+ = AlSi 3+ +HSi 1+ or SiSi 4+ = 4HSi 1+ , the substitution of H + in the aluminum octahedral site AlAl 3+ = 3HAl 1+ , the substitution of H + and H3O + in the potassium site (KK 1+ = HK 1+ , or KK 1+ = (H3O)K 1+ ) and accommodation of molecular H2O in the tunnel as in KAlSi3O8 + H2O = KAlSi3O8•H2O. The chemical analyses provide crucial information on possible coupled substitution mechanisms for hydrogen incorporation in liebermannite. The slight deficiency in silicon and increase in aluminum observed in some liebermannite grains (K0.957±0.044Al1.090±0.057Si2.931±0.032 O8) is consistent with a reaction involving a substitution in the distorted octahedral with the direct replacement of Si 4+ by Al 3+ and H + or by 4H + . The chemical analyses also indicate a slight aluminum deficiency in some crystals (K1.00Al0.983Si3.01 O8). This may suggest possible proton substitution in the octahedral site, replacing Al 3+ with 3H + .

As the hiP liebermannite is considered the dominant polymorph throughout the Earth's lower mantle, knowledge on water solubility in the hiP structure becomes crucial to assess the transport of H2O to the deepest lower mantle regions through deep extending slabs. Given the positive clapeyron slope of P (GPa) = 16.6 + 0.007 K (T) between the lowP and hiP liebermannite and assuming that trace quantity of water is unlikely to affect the slope, the hydrous liebermannite sample synthesized at 24 GPa and 1500 K is expected to be the lowP polymorph with tetragonal symmetry [START_REF] Nishiyama | Stability and P -V -T equation of state of KAlSi3O8 -hollandite determined by in situ Xray observations and implications for dynamics of subducted continental crust material[END_REF]. However, upon a decrease of temperature to 1000 K, at constant pressure, the sample is likely to transform to hiP liebermannite with monoclinic symmetry [START_REF] Nishiyama | Stability and P -V -T equation of state of KAlSi3O8 -hollandite determined by in situ Xray observations and implications for dynamics of subducted continental crust material[END_REF]. Upon structural transformation from lowP to hiP liebermannite, the electrical resistivity measurements of the sample do not indicate a discontinuous decrease of the electrical resistance, which indicates the release of free fluid phase or hydrous melting in the sample [START_REF] Freitas | Electrical conductivity of hydrous silicate melts : Implications for the bottom-up hydration of Earth ' s upper mantle[END_REF]. Assuming that the water incorporation into the liebermannite structure occurs at the highest temperature, this crucial observation suggests that the hiP liebermannite structure is capable of sequestering H + into octahedral sites as molecular H2O or ionic H3O + in its structure similar to the lowP liebermannite. The incorporation of H2O has been observed in both naturally occurring tetragonal and monoclinic hollandite crystal structures with different stoichiometry [START_REF] Miura | The crystal structure of hollandite[END_REF], thus confirming our observations on liebermannite. The micro Raman grid-analyses of our samples after electrical conductivity measurements indicate the absence of melt phases, confirming that the water expulsion may not occur in our samples at high temperatures.

Due to the small grain sizes (< 5 μm), infrared analyses in these minor phases were not possible; nevertheless, absence of peaks in the high-wavenumber region of the Raman spectra of wadeite and Al-free stishovite indicates that water is not present within their crystal structure (Fig. 7). However, H2O contents up to 3 wt. % have been observed in stishovite containing aluminum [START_REF] Litasov | High hydrogen solubility in Al-rich stishovite and water transport in the lower mantle[END_REF][START_REF] Yoshino | Electrical conductivity of stishovite as a function of water content[END_REF]. IR spectral features of Al-free stishovite exhibit some similarities with those of liebermannite in the spectral range from 3111 to 3261 [START_REF] Litasov | High hydrogen solubility in Al-rich stishovite and water transport in the lower mantle[END_REF]. Particularly peak position of 3117 cm -1 in the sample synthesized at 15 GPa (Fig. 7a) overlaps the main peak position of stishovite at 3117 cm -1 observed for Al-free and Al-bearing stishovite [START_REF] Yoshino | Electrical conductivity of stishovite as a function of water content[END_REF]. In order to confirm that the spectral feature at 3117 cm -1 in libermannite, we have analyzed the O-H region of the Raman spectra of stishovite in the 15 GPa sample. The observed water-free conditions (Fig 7a) strongly suggests the spectral feature at 3117 cm -1 observed in libermannite could be unrelated to stishovite, rather a spectral feature characteristic of libermannite. We conclude that the similarities btween libermannite and stishovite may be related to similar substitution mechanisms of hydrogen in octahedral sites SiO6 or AlO6.

The preferential partitioning of water into the liebermannite phase over stishovite observed in our samples may be partly due to the low Al content in stishovite that is in equilibrium with liebermannite (Table 2). Therefore, the presence of Al-free stishovite in the studied samples should not influence the IR spectra of liebermannite. Even in the case of very low water contents of about 16-30 wt. ppm observed in Al-free stishovite [START_REF] Litasov | High hydrogen solubility in Al-rich stishovite and water transport in the lower mantle[END_REF], these should not affect water quantification of liebermannite in this study.

Geophysical implications

The upper mantle underneath the Philippine Sea [START_REF] Tarits | The heterogeneous electrical conductivity structure of the lower mantle[END_REF] and Northeastern China [START_REF] Kelbert | Global electromagnetic induction constraints on transition-zone water content variations[END_REF] is characterized by unusual electrical and seismic wave velocity structures, in particular, electrical conductivities as high as 1 S/m and high-velocity anomalies of up to +δVP 1.5 % and δVS 2.0 % have been observed in these regions. The observed anomalies cannot be readily explained by normal mantle consisting of major mantle phases such as olivine and wadsleyite [START_REF] Manthilake | Electrical conductivity of wadsleyite as a function of temperature and water content[END_REF]. If the observed positive velocity anomalies were assumed to be caused by temperature variations alone, the temperature in the uppermost region of the mantle transition zone estimated to be about 800 K below the normal mantle geotherm [START_REF] Manthilake | Electrical conductivity of wadsleyite as a function of temperature and water content[END_REF]. The corresponding electrical conductivity of hydrous olivine and wadsleyite for such low temperatures would be significantly lower; as a consequence, the proton conduction in olivine and wadsleyite cannot be responsible for the electrical conductivity observed at the upper mantle and in the mantle transition zone [START_REF] Manthilake | Electrical conductivity of wadsleyite as a function of temperature and water content[END_REF]. Similarly, the presence of melt cannot explain the observed anomalies, because even a minor fraction of melt would significantly decrease both the primary (VP) and secondary (VS) wave velocities [START_REF] Chantel | Experimental evidence supports mantle partial melting in the asthenosphere[END_REF][START_REF] Freitas | Experimental evidence supporting a global melt layer at the base of the Earth's upper mantle[END_REF][START_REF] Weidner | Ultrasonic Acoustic Velocities During Partial Melting of a Mantle Peridotite KLB-1[END_REF].

Due to superionic conduction properties and high modal abundance (> 30 vol. %) of liebermannite in deeply subducted crustal components, the electrical conduction in slabs is likely to dominate by the liebermannite. This would make such slab components highly conductive compared to slabs consisting of mafic and ultramafic lithologies [START_REF] Huang | Water content in the transition zone from electrical conductivity of wadsleyite and ringwoodite[END_REF]Yoshino et al., 2008Yoshino et al., , 2008[START_REF] Yoshino | Electrical conductivity of stishovite as a function of water content[END_REF]. However, the presence of liebermannite alone cannot account for the positive velocity anomaly observed in the upper mantle beneath the Philippine Sea [START_REF] Tarits | The heterogeneous electrical conductivity structure of the lower mantle[END_REF] and northeastern China [START_REF] Kelbert | Global electromagnetic induction constraints on transition-zone water content variations[END_REF] because the elastic wave velocities measured in liebermannite are comparable with those of common mantle mineral phases [START_REF] Caracas | Elasticity of (K,Na)AlSi3O8 hollandite from lattice dynamics calculations[END_REF][START_REF] Mookherjee | Detecting deeply subducted crust from the elasticity of hollandite[END_REF]. Instead, the seismic wave velocities of Al-poor stishovite, another principal mineral phase stable in subducting continental sediments, are notably faster than those of mantle mineral phases [START_REF] Gréaux | Sound velocities of aluminum-bearing stishovite in the mantle transition zone[END_REF]. The stishovite co-existing with liebermannite in our samples is noticeably devoid of aluminum (Table 2). Here we propose that Al-poor stishovite coexisting with liebermannite in the subducted continental sediments may be responsible for the positive velocity anomaly observed at the mantle transition zone beneath the Philippine Sea [START_REF] Tarits | The heterogeneous electrical conductivity structure of the lower mantle[END_REF] and northeastern China [START_REF] Kelbert | Global electromagnetic induction constraints on transition-zone water content variations[END_REF].

Continental crustal portion of the subducted slab, consisting of libermannite, majoritic garnet, and stishovite, are expected to have a higher 0.2 g cm -3 density than the ambient mantle, but this density relation is overturned at the 660 km discontinuity [START_REF] Irifune | Subduction of continental crust and terrigenous and pelagic sediments: an experimental study[END_REF][START_REF] Nishiyama | Stability and P -V -T equation of state of KAlSi3O8 -hollandite determined by in situ Xray observations and implications for dynamics of subducted continental crust material[END_REF]. This implies that subducting crustal components of the slabs may penetrate into the mantle transition zone, but are prevented from entering the lower mantle [START_REF] Nishiyama | Stability and P -V -T equation of state of KAlSi3O8 -hollandite determined by in situ Xray observations and implications for dynamics of subducted continental crust material[END_REF], such as those observed in the Northeastern China subduction system [START_REF] Ichiki | Water content and geotherm in the upper mantle above the stagnant slab : Interpretation of electrical conductivity and seismic P-wave velocity models[END_REF]. However, the lower mantle origin of some of the plumes [START_REF] French | Broad plumes rooted at the base of the Earth's mantle beneath major hotspots[END_REF][START_REF] Hart | Mantle Plumes and Entrainment : Isotopic Evidence[END_REF] carrying the EM-type components in OIB magmas [START_REF] Hart | Mantle Plumes and Entrainment : Isotopic Evidence[END_REF] provides strong evidence supporting the presence of recycled continental materials in the lower mantle. We propose that mixing of continental with oceanic sediments may increase the density of sedimentary components facilitating the sinking of slabs towards the core-mantle boundary. The high electrical conductivity observed at a 1200 km depths in the Earth's lower mantle beneath the Philippine sea subduction system [START_REF] Kelbert | Global electromagnetic induction constraints on transition-zone water content variations[END_REF][START_REF] Tarits | The heterogeneous electrical conductivity structure of the lower mantle[END_REF]) may be linked to such continuous subduction of continental sediments across the mantle transition zone into the Earth's lower mantle (Fig. 7).

Owing to hotter temperatures relative to the ambient mantle, the rising plume is unlikely to exchange incompatible elements, including water, with the surrounding mantle, including the hydrated MTZ and the upper mantle [START_REF] Bercovici | Whole-mantle convection and the transition-zone water filter[END_REF]. Thus, it is likely that the H2O present in plume-derived rocks is primarily linked to the deep lower mantle source. The relatively high H2O contents of OIB derived from deep mantle plumes hint towards the presence of hydrated crustal components present in the deeper part of the lower mantle (Fig. 8). While major mineral phases in pyrolitic compositions are almost devoid of H2O under lower mantle conditions, our study demonstrates that liebermannite could be an important host of H2O in the Earth's lower mantle. In this scenario, H2O in the EMII-source OIB is primarily related to deeply subducted continental materials, in which liebermannite is the principal H2O carrier. On the other hand, the subducted continental materials that may reach the core-mantle boundary would hydrate the HIMU type plume sources together with a possible contribution from dense hydrous magnesium silicates (DHMS) [START_REF] Nishi | Stability of hydrous silicate at high pressures and water transport to the deep lower mantle[END_REF][START_REF] Pamato | Lower-mantle water reservoir implied by the extreme stability of a hydrous aluminosilicate[END_REF] and primordial H2O [START_REF] Dixon | Recycled dehydrated lithosphere observed in plume-influenced mid-ocean-ridge basalt[END_REF]. [START_REF] Syracuse | The global range of subduction zone thermal models[END_REF], and along the adiabatic geotherm in the lower mantle [START_REF] Katsura | Adiabatic temperature profile in the mantle[END_REF], respectively. The red shaded triangle, square, and the circle indicate the electrical conductivity of liebermannite at 12, 15, and 24 GPa, respectively. The gray shaded boxes indicate the conductivity profiles of NE-China [START_REF] Kelbert | Global electromagnetic induction constraints on transition-zone water content variations[END_REF] and the Philippine Sea [START_REF] Tarits | The heterogeneous electrical conductivity structure of the lower mantle[END_REF]. The electrical conductivity profiles of dry (black lines) and wet (blue lines) olivine, wadsleyite, ringwoodite (Yoshino et al., 2008), stishovite with 0.1 wt% of H2O [START_REF] Yoshino | Electrical conductivity of stishovite as a function of water content[END_REF] (pink lines) and majorite garnet (Yoshino et al., 2008) (light green lines) are shown for comparison. The electrical conductivity of (Mg, Fe) bridgmanite [START_REF] Katsura | Electrical conductivity of silicate perovskite at lowermantle conditions[END_REF]) (green dashed line) and (Mg, Fe, Al) bridgmanite [START_REF] Ohta | Electrical conductivities of pyrolitic mantle and MORB materials up to the lowermost mantle conditions[END_REF][START_REF] Sinmyo | Lower mantle electrical conductivity based on measurements of Al, Fe-bearing perovskite under lower mantle conditions[END_REF][START_REF] Yoshino | Electrical conductivity model of Al-bearing bridgmanite with implications for the electrical structure of the Earth's lower mantle[END_REF]) (blue dashed line), extrapolated along the adiabatic geotherm in the lower mantle [START_REF] Katsura | Adiabatic temperature profile in the mantle[END_REF], are shown for comparison. The numbers next to conductivity lines are the water contents in wt. %. However, due to the imbalance of the mass of H2O transported to the Mantle Transition Zone (MTZ) via subduction [START_REF] Hacker | H2O subduction beyond arcs[END_REF] and the estimated mass of H2O stored in MTZ based on H2O contents of constituent phases [START_REF] Freitas | Electrical conductivity of hydrous silicate melts : Implications for the bottom-up hydration of Earth ' s upper mantle[END_REF], extraction of H2O from MTZ to the lower mantle can be expected. Because of the ability of libermannite to accommodate several wt. % of the water in its structure (both in lattice defects and the tunnel structures), we argue that libermannite could be a possible carrier of water from the MTZ to the lower mantle. In this twostep, H2O conveyor process, the extraction of H2O by rehydration of deep penetrating slabs at the MTZ would be a likely scenario. Depending on gravitational stability, the subducting continental sediment component may detach from deep-extending slabs, rise through the mantle forming EMII-source OIB, in which liebermannite is the principal water carrier. Some continental sediments may reach the CMB and provide H2O for HIMU-source lower mantle-derived plumes.

(Lib: liebermannite, Sa: K-feldspar, Mj: majorite garnet, Brg: bridgmanite, Stv; stishovite).

List of tables

Table 2. Chemical composition of the constituent mineral phases.

FiguresFig. 1 .

 1 Figures

Fig. 2 .

 2 Fig. 2. The electrical conductivity of liebermannite as a function of inverse temperature. (a)

Fig. 3 .

 3 Fig. 3. Impedance spectra of the sample at different stages of heating at 12 GPa. (a) at 500 K

Fig. 4 .

 4 Fig. 4. Raman spectra of liebermannite at 12, 15, and 24 GPa.

Fig. 5 .

 5 Fig. 5. Mineral assemblages of the recovered sample. (a) A back-scattered electron image

Fig. 6 .

 6 Fig. 6. FTIR spectra of liebermannite. (a) Near IR spectra of samples synthesized at 12 (green

Fig. 7 .

 7 Fig. 7. Raman spectra of minor phases observed at the sample synthesized at 15 GPa. The

Fig. 8 .

 8 Fig. 8. Comparison of electrical conductivities of liebermannite with global and regional

Fig. 9 .

 9 Fig. 9. The potential mechanism explaining the elevated H2O content in ocean island basalts. K-feldspar in subducting continental sediments transforms to liebermannite above 8 GPa and continues to be stable down to the core-mantle boundary conditions. Most of the H2O carried to the lower mantle by deep penetrating slabs may have been incorporated into the slab at the surface.

Table 1 .

 1 The reference pre-exponential factor Log 0 and the activation enthalpy in eV for liebermannite at 12, 15, and 24 GPa and their respective water contents.

	Pressure	Temperature	Log*	Log**	H *	H **	H2O contents
	(GPa)	(K)	(S/m)	(S/m)	(eV)	(eV)	ppm wt.
	12	550 -750	0.57	n.a	0.37	n.a.	1044 ± 87
		750 -1100	2.19	2.32	0.64	0.78
		>1100	3.14	3.88	0.85	0.82
	15	573 -683	0.58	n.a	0.4	n.a.	885 ± 78
		773 -1173	1.76	1.61	0.64	0.76
		> 1173	2.78	3.96	0.92	0.83
	24	1040 -1200	0.38	0.24	0.63	0.76	1104 ± 18
		> 1270	5	8.26	1.78	0.89
	* parameters from individual fitting
	** parameters from single fit
	n.a: not nalyzed

Table 2 .

 2 Chemical composition of the constituent mineral phases.

			Liebermannite		Wadeite	Stishovite	Kyanite
		12 GPa	15 GPa	24 GPa	12 GPa	15 GPa	12 GPa	15 GPa	12 GPa	15 GPa
	SiO2	64.44 (0.7)	63.91 (2.5)	66.70 (0.2)	73.00 (0.9)	71.82 (0.5) 100.50 (0.8) 101.13 (0.5)	38.23 (0.9)	39.11 (0.9)
	K2O	16.47 (0.4)	16.18 (1.2)	16.42 (1.8)	27.16 (0.8)	28.22 (0.6)	0.06 (0.3)	0.13 (0.07)	0.15 (0.06)	0.61 (0.4)
	Al2O3	18.98 (0.3)	20.4 (2.9)	17.21 (2.0)	0.15 (0.1)	0.6 (0.5)	0.05 (0.2)	0.07 (0.03)	62.6 (1.2)	61.35 (1.2)
	Total	99.89	100.49	100.33	100.31	100.64	100.62	101.33	100.98	101.07
	The standard deviation indicated in parenthesis.
	850

Acknowledgments

We thank Nathalie Bolfan-Casanova and Maud Boyet for constructive discussion, J-L. Devidal for his assistance with electron microprobe analyses and Emmy Voyer for her assistance with SEM analyses. The constructive and in-depth reviews by the associate Editor, and two reviewers: Lidong Dai and Hongzhan Fei substantially improved the quality of the manuscript. G.M. acknowledges funding from the(INSU-CNRS. M.M. is supported by the US National Science Foundation grants, EAR 1753125, 1763215. This research was financed by the French Government Laboratory of Excellence initiative n°ANR-10-LABX-0006, the Région Auvergne, and the European Regional Development Fund (ClerVolc contribution number xx). The authors comply with AGU's data policy, and the data of this study are available in Figshare. (https://figshare.com/articles/Water_content_calculations/12179067) and (https://figshare.com/articles/Experimental_data_xlsx/12179046).

The authors declare no competing financial interests.