Sylvain Boussier
email: sylvain.boussier@ema.fr

Michel Vasquez
email: michel.vasquez@ema.fr

Yannick Vimont
email: yannick.vimont@ema.fr

Saïd Hanafi
email: said.hanafi@univ-valenciennes.fr

Philippe Michelon
email: philippe.michelon@univ-avignon.fr

Solving the 0-1 Multidimensional Knapsack Problem with Resolution Search

Keywords: 0-1 Multidimensional Knapsack Problem, Resolution Search, Branch & Bound

We propose an exact method which combines the resolution search and branch & bound algorithms for solving the 0-1 Multidimensional Knapsack Problem. This algorithm is able to prove large-scale strong correlated instances. The optimal values of the 10 constraint, 500 variable instances of the OR-Library are exposed. These values were previously unknown.

Introduction

In this article, we present a new exact method which hybridizes the resolution search of [START_REF] Chvátal | Resolution search[END_REF] and a branch & bound algorithm inspired by a previous work from [START_REF] Vimont | Reduced costs propagation in an efficient implicit enumeration for the 01 multidimensional knapsack problem[END_REF] for the 0-1 Multidimensional Knapsack Problem (01MKP). The 01MKP is a wellknown opitmization problem which can be viewed as a resource allocation model and can be stated as follows:

(P ) Maximize n j=1 c j x j (1) subject to n j=1 a ij x j ≤ b i i = 1, ..., m (2)
x j ∈ {0, 1} j = 1, ..., n

where n is the number of items, m is the number of knapsack constraints with capacities b i (i = 1, 2, ..., m), c ∈ N n , A ∈ N m×n and b ∈ N m . Each item j (j = 1, 2, ..., n) yields c j units of profit and consumes a given amount of resource a ij for each knapsack i.

The proposed approach is centered on two main 01MKP results: (i) the consideration of a reduced costs constraint based on the reduced costs at the optimality of the problem's LPrelaxation (see [START_REF] Balas | Pivot and complement a heuristic for zero-one programming[END_REF]; [START_REF] Oliva | Constraint and linear programming : Using reduced costs for solving the zero/one multiple knapsack problem[END_REF]) and (ii) the decomposition of the search space in several hyperplane where the number of items to choose is fixed at a given integer value (see [START_REF] Vasquez | An hybrid approach for the 0-1 multidimensional knapsack problem[END_REF]; [START_REF] Vimont | Reduced costs propagation in an efficient implicit enumeration for the 01 multidimensional knapsack problem[END_REF]). Our algorithm is self-sufficient and does not require any lower bound as starting value. We show that the structure of resolution search enables to explore partially and iteratively different subproblems (hyperplanes with a fixed number of items) while keeping completeness. This way of exploration enhances the diversification of the search and permits to improve the lower bound rapidly. For each hyperplane, this lower bound, associated with the upper bound (tighter than the classical one given by the LP-relaxation), enforces the strength of the so-called reduced costs constraint widely used in the algorithm. Roughly speaking, if we consider that the variables are sorted in decreasing order of their reduced costs (basic variables at the bottom of the list), the enumeration of the first variables is carried out by resolution search where the enumeration of the remaining variables is tackled by our specific branch & bound. Our approach proved all the optimality of all the 10 constraint, 500 variable instances. These optimal solutions were previously unknown.

Resolution search

Resolution search was proposed by [START_REF] Chvátal | Resolution search[END_REF] as an alternative to branch & bound for mixed 0-1 linear programming problems. This approach, based on an original exploration of the search space, uses the information brought by the fails that have occurred during the search to progressively shrink the search tree. Each time a terminal node is reached, a minimal partial instantiation of the variables responsible for the fail is identified. This partial instantiation, which corresponds also to a terminal node, is then recorded in a specific way in order to discard the corresponding subtree from the search space and to provide the next node in the exploration. The specificity of the recording mechanism allows the algorithm to preserve space memory while keeping completeness.

Globally, resolution search is composed of two main elements: (i) a set of partial instantiations (denoted path-like family F ) corresponding to fails encountered during the search and recorded as boolean clauses and (ii) an function called obstacle which uses the information brought by the family F to explore promising parts of the search space. From a partial instantiation u(F ) derived from F , which is not an already explored instantiation, obstacle performs two different phases:

1. The waxing phase which replaces step by step the free components u(F ) j = * by 0 or 1 until a terminal node (denoted u * ) is reached (a terminal node is a partial instantiation of the variables which does not lead to any improvement of the best known lower bound).

2. The waning phase which tries to find a minimal element S of {0, 1, * } n such that S ⊑ u * and S is a terminal node (formally, let u and v be two vectors of {0, 1, * } n , we call u an extension of v and we write u ⊑ v, if v j = u j whenever u j = * ).

Once the minimal element S is identified, it is recorded in a specific way as a clause in F . The structure of F enables both to guide the search toward prosimising parts of the search space and to guaranteed the completeness using the resolution / refutation principle. Indeed, if the minimal element S correspond to the instantiation ( * , * , ..., * ) then the best known solution is an optimal soltion of P . The algorithm 1. represents the function Resolution Search which takes as parameter a lower bound (LB) of the problem.

General approach

The 01MKP is tackled by decomposing the problem in several subproblems where the number of items to choose is fixed at a given integer value. Considering that LB is a lower bound of the problem, we can define

X = {x | A • x ≤ b, c • x ≥ LB + 1, x ∈ {0, 1} n }
as the set of feasible solutions strictly better than LB and consider the following problems:

P + : maximize {1 • x | x ∈ X} and P -: minimize {1 • x | x ∈ X}
where 1 is the vector of dimension n with all its components equal to one and

X = {x | A • x ≤ b, c • x ≥ LB + 1, x ∈ [0, 1] n }.
In what follows, we note v(P ) the optimal value of a given

(P ) = max{v(P k ) | k min ≤ k ≤ k max } where P k : maximize {cx | x ∈ X, 1 • x = k}.
Solving the 01MKP by tackling separately each of the subproblems P k for k = k min , ..., k max appeared to be an interesting approach [START_REF] Vasquez | An hybrid approach for the 0-1 multidimensional knapsack problem[END_REF]; [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF]; [START_REF] Vimont | Reduced costs propagation in an efficient implicit enumeration for the 01 multidimensional knapsack problem[END_REF]) particularly because the additional constraint (1 • x = k) provides tighter upper bounds than the classical LP-relaxation. This approach takes also into consideration a constraint based on the reduced costs of the non-basic variables. Let us consider the upper bound UB = c • x, where x is the optimal solution of the LP-relaxation. Let c be the vector of the reduced costs and N the indexes of the non-basic variables. If we know a lower bound LB ∈ N of the problem, then each better solution x must satisfy the following constraint:

j∈N |x j =0 |c j |x j + j∈N |x j =1 |c j |(1 -x j ) ≤ UB -LB (4)
(the reader is referred to [START_REF] Balas | Pivot and complement a heuristic for zero-one programming[END_REF]; [START_REF] Oliva | Constraint and linear programming : Using reduced costs for solving the zero/one multiple knapsack problem[END_REF] and [START_REF] Vimont | Reduced costs propagation in an efficient implicit enumeration for the 01 multidimensional knapsack problem[END_REF] for more details on the subject). The use of the constraint (4), also called the reduced costs constraint, is twofold: (i) it enables us to identify the unfeasibility (with respect to the constraint c • x ≥ LB + 1) of some partial solutions and (ii) it allows us to fix those variables x j with |c j | > UB -LB at their optimal value xj . Both the tight upper bound and the good lower bound enhance the efficiency of the reduced costs constraint.

To summarize, our approach consists of decomposing the search space in several hyperplanes then solving the LP-relaxation for each one in order to generate the reduced costs constraint. The exploration is then carried out partially and iteratively for each hyperplane with resolution search until all the search space is explored.

Indeed, thanks to the structure of resolution search, P can be solved by progressively exploring each of the hyperplanes (1 • x = k, k ∈ N). Let CFamily[k] be the family of clauses associated to resolution search for the subproblem P k . At each step of the search, CFamily[k] provides all the information about the state of the search: the terminal nodes recorded at this stage and the next node to explore. It is possible to execute some iterations of resolution search at a given P k , then continue to another P k ′ and go on with the subproblem P k again without any loss of information. The algorithm 2 detail the Iterative RS algorithm which corresponds to resolution search limited to a given number of iterations and the algorithm 3 shows the hyperplane exploration using Iterative RS. This way of exploration enforces the diversification and the convergence of the search. Note that, in algorithm 3, greedy is a simple greedy function used for providing a first lower bound. Starting with the node u(F ) given by the path-like family F , obstacle replaces stepby-step the components * of u(F ) by 0 or 1 which constructs the node u + . If u + is a terminal node, the function provides a minimal clause S such that S ⊑ u * and S is a terminal node. Our implementation of obstacle is based on the reduced-costs constraint (4) presented in section 3. At the beginning of the procedure, the LP-relaxation of the problem is solved for each available hyperplane in order to give us the information needed for the reduced costs constraint. Then obstacle proceeds to the following steps:

• The first step, called consistency phase, consists of checking the feasibility of u(F ).

Initially, we define a value gap = UB -LB. If a constraint is violated, the descent phase stops and the corresponding partial instantiation is recorded as the clause S in F . At the same time, the reduced costs constraint ( 4) is checked and the gap value is updated: for each non-basic variable set at the opposite of its optimal value (1xj ), its reduced costs (c j ) is subtracted from the gap. If it happens that gap < 0 then the current partial solution is a terminal node. In this case, S is only composed of the variables set at the opposite of their optimal value in u(F ).

• If u(F ) is feasible, we go to the next step which we call implicit waning phase. This phase consists of branching on all the remaining free variables with a reduced cost greater than gap. Those variables must be fixed at their optimal value for satisfying the reduced costs constraint (4). The branching decisions taken in this phase are just a consequence of the instantiated variables in u(F ). Consequently they can be removed from S.

• Then the algorithm starts the so-called waxing phase which consists of assigning values to free variables. The chosen strategy is to select the free variable with the greater absolute reduced cost value and to assign its optimal value xj to it. Obviously, each time a branching is done, the feasibility of the current partial solution is verified and in case of fail, the waxing phase stops and the corresponding clause S is added to F . Note that the variables set in the implicit waning phase are still not taken into account in S.

• Once the number of remaining free variables is less than or equal to a given number spb size, the waxing phase stops and the corresponding subproblem is solved with a branch & bound algorithm. This subproblem includes the free variables with the lowest reduced cost and the basic variables. Obviously, since the branch & bound algorithm explores the whole subtree corresponding to these variables, the clause S does not contain any branching choices made during this phase. Only the branching decisions taken during the consistency phase and/or the waxing phase are considered.

• The algorithm used to enumerate the variables with the lowest reduced cost, and the basic variables, is widely inspired by a previous one published by [START_REF] Vimont | Reduced costs propagation in an efficient implicit enumeration for the 01 multidimensional knapsack problem[END_REF]. As it is represented in figure 4 (which summarizes the exploration process), this algorithm embeds a specific Depth First Search (dfs) procedure for solving small subproblems with the last 20 variables.

Computational results

Our algorithm has been experimented on the well known OR-Library set of 01MKP instances proposed by [START_REF] Beasley | Or-library: Distributing test problems by electronic mail[END_REF]. Each instance of this set is denoted cbm.n r where m is the number of constraints, n the number of variables and r the instance's number. The 01MKP OR-Library benchmark is composed of 5, 10 and 30 constraint instances and 100, 250 and 500 variables instances. There are 9 set of 30 instances with each n × m combination.

We obtained better proof time than the exact published approaches of [START_REF] James | Enumeration methods for repeatedly solving multidimensional knapsack sub-problems[END_REF], [START_REF] Vimont | Reduced costs propagation in an efficient implicit enumeration for the 01 multidimensional knapsack problem[END_REF] and the commercial software CPLEX 9.2, for the 10 constraint, 250 variable instances and the 5 constraint, 500 variable instances. We observed also that for the 5 constraint, 500 variable instances, the times required for obtaining the best solutions are clearly better than the ones provided by the best known heuristics on these instances [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF] and [START_REF] Wilbaut | New convergent heuristics for 0-1 mixed integer programming[END_REF]). Moreover, our algorithm proved the optimality of all the cb10.500 instances. The corresponding optimal values, which were previously unknown, are exposed in table 1. The column z opt is the optimal value of the instance, opt. (h) is the time in hours required for obtaining the optimal value and proof (h) is the time in hours for proving the optimality of the value z opt . The column z optz corresponds to the gap between the optimal value and the previously best known solution: (vv) indicates that z was first found by [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF] and (wh) indicates that z was first found by [START_REF] Wilbaut | New convergent heuristics for 0-1 mixed integer programming[END_REF].

Conclusion

Although our implementation of resolution search is quite far from the original one proposed by Chvátal, especially because we have exploited specific structures of the 01MKP (hybridization with branch & bound, 1.x = k hyperplane decomposition, reduced cost constraint to generate implicitly the partial instantiations responsible for the fails, etc.), we showed that resolution search is a promising framework for designing efficient algorithms. Since the proof times are long for the 10 constraint, 500 variable instances, we plan to improve again our algorithm in order to accelerate the resolution process in the hope of maybe closing the OR-Library 01MKP benchmark by solving the 60 last 30 constraint 250/500 variable instances.

Figure 1 :

 1 Figure 1: Resolution Search Algorithm Resolution Search(LB) { F = ∅; u = ( * , * , ..., * ); while(1) { try = obstacle(u,LB,S); if(try > LB) LB = try; add S to F and update F ; u = u(F ); if(( * , * , ..., * ) ∈ F ) Break; } } problem P . Let k min = ⌈v(P -)⌉ and k max = ⌊v(P + )⌋, then we have v(P ) = max{v(P k ) | k min ≤ k ≤ k max } where

Figure 2 :

 2 Figure 2: Iterative resolution search

Figure 4 :

 4 Figure 4: Global view of the exploration process of an hyperplane

Table 1 :

 1 Results obtained on the 10 constraint, 500 variable instance of the OR-Library

	Instance	z opt	opt. (h)	proof (h)	z opt -z	Instance	z opt	opt. (s) proof (h)	z opt -z
	cb10.500 0	117821	24,5	567,2	+10(vv)	cb10.500 15	215086	0	43,9	0
	cb10.500 1	119249	68,4	272,9	+17(vv)	cb10.500 16	217940	13,4	36,1	0
	cb10.500 2	119215	18,6	768,3	0	cb10.500 17	219990	150,8	348,8	0
	cb10.500 3	118829	47,4	89,6	+4l(wh)	cb10.500 18	214382	12,7	57,8	+7(vv)
	cb10.500 4	116530	86,1	2530,3	+16(wh)	cb10.500 19	220899	0,2	21,3	0
	cb10.500 5	119504	2,3	188	0	cb10.500 20	304387	6,6	8,2	0
	cb10.500 6	119827	2,7	128	0	cb10.500 21	302379	0	8,4	0
	cb10.500 7	118344	161,7	179,6	+11(wh)	cb10.500 22	302417	67,2	105,5	+1(vv)
	cb10.500 8	117815	86,3	219,9	0	cb10.500 23	300784	0,9	3,8	0
	cb10.500 9	119251	3,1	354,9	0	cb10.500 24	304374	0,1	16,8	0
	cb10.500 10	217377	0	515,8	0	cb10.500 25	301836	29,7	30,9	0
	cb10.500 11	219077	0,5	437,6	0	cb10.500 26	304952	0	18,5	0
	cb10.500 12	217847	0	5,5	0	cb10.500 27	296478	1,1	9,3	0
	cb10.500 13	216868	0	104,4	0	cb10.500 28	301359	8,1	39,1	0
	cb10.500 14	213873	59,4	1382,1	+14(vv)	cb10.500 29	307089	1,2	4,4	0


Damien Daval
email: damien.daval@univ-grenoble-alpes.fr

Gaël Choblet

Christophe Sotin

François Guyot

Theoretical considerations on the characteristic timescales of hydrogen generation by serpentinization reactions on Enceladus

Keywords: serpentinization, Fe-bearing silicate weathering, iron oxidation, icy moon, hydrogen production, kinetic modeling

Key points:

-The dissolution of Fe-silicates in the Enceladus' core takes place at far-from-equilibrium conditions, resulting in elevated reaction rates.

-The duration of H 2 generation resulting from serpentinization of Enceladus' core does not exceed 500 Myr.

-Alternative processes must be explored to explain the observed rate of H 2 emission in 26 Enceladus' plume. 27

Plain language summary

Enceladus is a small icy moon orbiting Saturn that presents a striking activity with jets emanating from parallel fractures in the tectonized south polar terrain. The amplitude of the moon's rotation indicates the presence of a global internal ocean, and chemical analyses by the Cassini spacecraft of the materials ejected from the interior suggest the ongoing occurrence of high temperature water-rock interactions, possibly in the porous rocky core. In particular, this conjecture stems from the detection of molecular hydrogen (H 2 ), a well-known by-product of the aqueous alteration of iron-bearing silicates. The lifetime of such reactions is key to assess both the habitability potential of Enceladus and to constrain plausible durations of the present active stage in a context where the evolution of the moon is debated. Here, we modeled such lifetimes, exploring various scenarios for the rate-limiting mechanisms of the process. Overall, we found that the timescales of hydrogen production never exceed 500 Myr, and are ultimately controlled by the water supply in the core. We thus suggest that either the hydrothermal activity has developed recently on Enceladus, or that alternative processes including microbial activity must be tested to explain the observed H 2 flux in Enceladus' plume.

Introduction

Many planetary bodies in the outer solar system are now considered to be ocean worlds, i.e., they harbor global layers of liquid water beneath the ice crust (see for instance the synthesis by [START_REF] Nimmo | Ocean worlds in the outer solar system[END_REF]; also [START_REF] Hendrix | The NASA roadmap to ocean worlds[END_REF]). In most cases, inner oceans are believed to be in direct contact with rocky cores: only in the case of the largest moons Ganymede, Titan and Callisto, can a deep mantle composed of high-pressure ice polymorphs be present (e.g., [START_REF] Hussmann | Interiors and evolution of icy satellites[END_REF]). While chemical reactions between rocks and liquid water would raise the astrobiological potential of such habitats ( [START_REF] Lammer | What makes a planet habitable?[END_REF]), very little conclusive evidence exist. Although alteration by water may occur at low-temperature (e.g. [Zolotov, 2009]), higher temperatures are generally held as more favorable for chemical weathering. The observation of plumes at the surface of Enceladus ( [START_REF] Porco | Cassini observes the active south pole of Enceladus[END_REF]) or Europa [START_REF] Roth | Transient water vapor at Europa's south pole[END_REF]) is interpreted as evidence of warmer interiors. Interior models indicate plausible contexts that differ for the two moons: magmatism may exist in the sub-surface of the relatively large rocky core (about 1450 km in radius) of Europa [START_REF] Běhounková | Tidally Induced Magmatic Pulses on the Oceanic Floor of Jupiter's Moon Europa[END_REF]), and warm convection of liquid water throughout the porous small core (less than 200 km in radius) of Enceladus [START_REF] Choblet | Powering prolonged hydrothermal activity inside Enceladus[END_REF]). In both cases, tidal heating modulates the thermal state and controls the efficiency of water-rock interactions. As tidal dissipation is itself strongly coupled with the orbital evolution of these moons (including resonances with neighbors), non-monotonic evolutions can be expected.

The specific case of Enceladus is the most emblematic. Clues indicative of present-day hydrothermalism have been witnessed by the Cassini spacecraft: nanoscale silica grains in the E-ring originating in Enceladus ( [START_REF] Hsu | Ongoing hydrothermal activities within Enceladus[END_REF][START_REF] Sekine | High-temperature water-rock interactions and hydrothermal environments in the chondrite-like core of Enceladus[END_REF]) as well as the identification of an endogenous source of molecular hydrogen (H 2 ) ( [START_REF] Waite | Cassini finds molecular hydrogen in the Enceladus plume: evidence for hydrothermal processes[END_REF]) and the subsequent detection of complex organics in the ejected icy grains ( [START_REF] Postberg | Macromolecular organic compounds from the depths of Enceladus[END_REF]) all point to ongoing water-rock interactions at temperatures higher than 50 °C ( [Hsu et al., 2015;[START_REF] Sekine | High-temperature water-rock interactions and hydrothermal environments in the chondrite-like core of Enceladus[END_REF]).

Because of the chondritic nature of Enceladus' core, the aqueous alteration of Fe-bearing minerals (such as olivine or orthopyroxenes) is often considered as a model reaction that can be responsible for the production of molecular hydrogen. The most generic chemical equation of this process, without any preconception on the nature of the secondary hydrated silicates that may form, combines the reaction of olivine dissolution:

Mg 2-x Fe x SiO 4 + 2H 2 O = SiO 2 (aq) + (2-x)Mg 2+ + xFe 2+ + 4OH - (1) to that of the production of molecular hydrogen resulting from Fe 2+ oxidation:

3Fe 2+ + 6OH -= Fe 3 O 4 + 2H 2 O + H 2 (2)
which ultimately leads to the following overall reaction:

A similar overall reaction can be retrieved by using a generic orthopyroxene as a reactant instead of olivine: The observation of such an active stage at present for Enceladus has strong implications.

The timescales during which H 2 can be produced by rock alteration are key aspects because (i) abiotic hydrogen can fuel life ( [START_REF] Ménez | Abiotic hydrogen and methane: fuels for life[END_REF]) and (ii) these would place major constraints on the plausible evolution of the moon in a context where the general history of the moon, let alone its age, is debated (see [START_REF] Ida | The origin of Saturn's rings and moons[END_REF]).

Based on 1D-reactive transport simulations conducted with the PHREEQC code ( [START_REF] Parkhurst | Description of input and examples for PHREEQC version 3-a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations[END_REF]), [START_REF] Zandanel | Short lifespans of serpentinization in the rocky core of Enceladus: Implications for hydrogen production[END_REF] recently suggested short timescales for H 2 generation, never exceeding 300 Myr at most. Their simulations supposed that olivine dissolution was the rate-limiting step of the serpentinization process (and therefore, of the associated production of H 2 ) and relied on parameters taken from the literature for fluid flow, olivine grain size, core porosity and olivine dissolution kinetics. Regarding this latter point, although the simulations were based on state-of-the-art dissolution kinetics parameters for the impact of e.g., pH or temperature, they were also based on rate laws derived from the transition state theory (TST). Whether such TST-based rate laws can be used to provide reliable estimates of the actual duration and intensity of fluid-mineral reactions in the field has been a source of lively debate in the geochemical community for over 40 years. Since the pioneering work of [START_REF] Paces | Rate Constants of Dissolution Derived from the Measurements of Mass Balance in Hydrological Catchments[END_REF], a significant body of literature has demonstrated that such rate laws usually overestimate the reactivity of silicates in the field by several orders of magnitude (e.g. [START_REF] Daval | Time-dependent feldspar dissolution rates resulting from surface passivation: Experimental evidence and geochemical implications[END_REF][START_REF] Maher | Rates of silicate dissolution in deep-sea sediment: In situ measurement using U-234/U-238 of pore fluids[END_REF][START_REF] Taylor | Relation between Soil Age and Silicate Weathering Rates Determined from the Chemical Evolution of a Glacial Chronosequence[END_REF][START_REF] Velbel | Constancy of silicate-mineral weathering-rate ratios between natural and experimental weathering: implications for hydrologic control of differences in absolute rates[END_REF][START_REF] White | Chemical weathering rates of a soil chronosequence on granitic alluvium .1. Quantification of mineralogical and surface area changes and calculation of primary silicate reaction rates[END_REF][START_REF] White | The effect of time on the weathering of silicate minerals: why do weathering rates differ in the laboratory and field?[END_REF]Art F. White et al., 2017;Zhu, 2005]). Such a well-known "field-lab discrepancy" is often bypassed in modeling exercises by decreasing the so-called "reactive" or "effective" surface area term by several orders of magnitude (e.g., [START_REF] Aradóttir | Multidimensional reactive transport modeling of CO 2 mineral sequestration in basalts at the Hellisheidi geothermal field, Iceland[END_REF][START_REF] Montes-H | Modelling of geochemical reactions and experimental cation exchange in MX80 bentonite[END_REF]). However, this strategy is not devoid of risks, because it does not account for the mechanistic origins of the "field-lab discrepancy". With this in mind, one can wonder the extent to which the short timescales of molecular hydrogen generation by serpentinization reactions on Enceladus suggested in previous studies (e.g., [START_REF] Zandanel | Short lifespans of serpentinization in the rocky core of Enceladus: Implications for hydrogen production[END_REF]) could not originate from the use of TST-based rate laws, which do not account for the rate-controlling role exerted by the spontaneous modifications of reacted mineral surface features (e.g., etch pit nucleation or development of passivating layers), ultimately leading to underestimated timescales of serpentinization.

Therefore, a detailed kinetic analysis is required to evaluate if the conclusions drawn in such previous studies may be a consequence of the assumptions they made to describe olivine dissolution kinetics.

In that respect, the objective of the present paper is to widen the spectrum of this previous study by investigating extensively all the known alternative possible rate-limiting mechanisms (and correspondingly, alternative rate laws) governing fluid-rock interactions in conditions relevant to Enceladus. A stepwise approach is described below, starting from the well-known 'lifetime diagram' approach developed by [START_REF] Olsen | Using a mineral lifetime diagram to evaluate the persistence of olivine on Mars[END_REF]. We show that although this approach may seem incomplete (as it does not account for the impact of the Gibbs free energy of reaction with respect to olivine (or orthopyroxene) dissolution), practically, it actually describes the intrinsic reactivity of primary minerals following the same rates as reactive transport simulations where secondary phases are supposed to precipitate at equilibrium (as in the work of [START_REF] Zandanel | Short lifespans of serpentinization in the rocky core of Enceladus: Implications for hydrogen production[END_REF]). Then, based on assumptions consistent with the literature regarding the fluid chemistry of Enceladus, we screen alternative ratelimiting mechanisms (and corresponding rate laws) such as the impact of (i) etch pit nucleation; (ii) surface layers; (iii) diffusion in nanoporous secondary assemblages; (iv) fluid flow. Importantly we restrict our analysis to fluid compositions buffered by mineral assemblages resulting from serpentinization reactions, as previously shown for high temperature experiments (T ≥ 120 °C) by e.g. [START_REF] Hsu | Ongoing hydrothermal activities within Enceladus[END_REF] or [START_REF] Sekine | High-temperature water-rock interactions and hydrothermal environments in the chondrite-like core of Enceladus[END_REF]. We finally follow an inverse problem strategy to relate today's H 2 flux to the lifetime of primary Fe-bearing silicates. In a nutshell, this in-depth analysis confirms the initial estimates by [START_REF] Zandanel | Short lifespans of serpentinization in the rocky core of Enceladus: Implications for hydrogen production[END_REF] based solely on transition state-derived rate laws; the reasons for this rather unexpected agreement will be thoroughly discussed.

The main assumptions of the present work can be summarized as follows: (i) the release of Fe from primary Fe-bearing silicates drives the rate of the whole process; (ii) (in connection with (i)) the fluid composition is supposed to be controlled by an equilibrium with respect to secondary phases consistent with the literature; (iii) reaction-induced cracking, which may occur during serpentinization because of the positive volume variation associated to serpentinization (e.g. [START_REF] Klein | Quantifying the volume increase and chemical exchange during serpentinization[END_REF] and references therein), is omitted, as it would result in a decrease of the lifetime of the dissolving Fe-bearing silicates and (iv) finally, for the sake of simplicity, the kinetic parameters related to the pure Mg endmembers forsterite (Mg 2 SiO 4 ) and enstatite (MgSiO 3 ) will be considered for the calculations, because of the lack of reliable dissolution kinetics studies of Fe-rich olivine and pyroxenes in the basic pH range.

Notwithstanding, it is generally admitted that the reactivity of Fe-bearing olivine and pyroxene is greater than that of their magnesian counterparts (e.g. [START_REF] Casey | CONTROL OF DISSOLUTION RATES OF ORTHOSILICATE MINERALS BY DIVALENT METAL OXYGEN BONDS[END_REF][START_REF] Daval | Fayalite (Fe2SiO4) dissolution kinetics determined by X-ray absorption spectroscopy[END_REF]), ensuring that this choice does not artificially result in an underestimation of olivine and enstatite lifetimes (and associated timescales of H 2 production).

Far-from-equilibrium olivine and enstatite dissolution rates

One of the simplest approaches to evaluate the duration of the dissolution of a given mineral is to consider that this mineral consists of spherical particles fully immersed in a large excess of water of homogeneous composition. In fact, considering a high ratio between the volume of water and the mass of mineral also corresponds to the shortest durations for the process, as it neglects the approach of the fluid towards close-to-equilibrium conditions, where the dissolution rate eventually drops to zero.

To the best of our knowledge, such an approach has been initially developed by [START_REF] Olsen | Using a mineral lifetime diagram to evaluate the persistence of olivine on Mars[END_REF] in the context of the Martian surface to estimate the maximum contact times for olivine grains with water before they dissolve away completely. Applying such a reasoning to the context of Enceladus would not seem absurd, since the very low pressure in the core (~ 40 MPa) suggests it is highly porous (up to 20-30%) and permeable to liquid water [START_REF] Choblet | Powering prolonged hydrothermal activity inside Enceladus[END_REF], so that it can be modeled as an aggregate of individual spherical particles representing the grains that make up the core. The core is assumed to have a chondritic composition with olivine, pyroxenes, native iron, and iron-nickel alloys being the main phases likely to contribute to H 2 release during water-rock interaction [START_REF] Holm | Serpentinization and the formation of H2 and CH4 on celestial bodies (planets, moons, comets)[END_REF]. The corrosion rate of iron and iron alloys being orders of magnitude greater than the dissolution rate of silicates (e.g., [START_REF] Mundhenk | Corrosion of carbon steel and the passivating properties of corrosion films formed under high-PT geothermal conditions[END_REF]), these phases will not be considered in the kinetic treatment described below.

[ [START_REF] Olsen | Using a mineral lifetime diagram to evaluate the persistence of olivine on Mars[END_REF] showed that the lifetime (t) of spherical grains can be modeled following:

𝑡 = 𝑟 𝑔𝑟𝑎𝑖𝑛 𝑉 𝑖 ̅ . 𝑅 𝑖 (5)
where r grain [cm] is the radius of the individual grains, 𝑉 𝑖 ̅ and R i are the molar volume [cm 3 /mol] and dissolution rate [mol/cm²/s] of mineral i (olivine or enstatite), respectively (see Appendix A for the full list of parameters used in the present study). The average size of chondrules in a collection of carbonaceous chondrites was found to be 74 µm in radius, ranging between 25 and 500 µm [START_REF] Rubin | Size-frequency distributions of chondrules in CO3 chondrites[END_REF]. Supposing that the size of chondrules is indicative of the mineral grain size, in the calculations detailed below, we propose to apply a two-fold factor on the lower and upper bounds of this distribution (i.e., grain sizes ranging from 10 µm to 1 mm radius) to make sure that we cover more than the entire range of expected grain sizes in a typical chondrite.

In the basic pH range, the far-from-equilibrium dissolution rate laws of enstatite and olivine can be written:

𝑅 𝑖 = 𝑘 0 . exp(-𝐸 𝑎 𝑅𝑇 ⁄ ) . (H + ) 𝑛 (6)
where k 0 [mol/cm²/s] is the dissolution rate constant, E a [J/mol] is the activation energy, R is the gas constant [mol/J/K], T the temperature [K], (H + ) is activity of protons and n the reaction order. The pH of Enceladus' ocean was suggested to range between 9 and 11 ([ [START_REF] Glein | The Geochemistry of Enceladus: Composition and Controls, Enceladus and the Icy Moons of Saturn[END_REF][START_REF] Glein | The pH of Enceladus' ocean[END_REF]) and the temperature of downwelling fluids may range from ~ 0 °C at the seafloor to up to 200 °C at the core center ([C Glein et al., 2018;[START_REF] Sekine | High-temperature water-rock interactions and hydrothermal environments in the chondrite-like core of Enceladus[END_REF]). Using the rate parameters derived from [START_REF] Rimstidt | Systematic review of forsterite dissolution rate data[END_REF] for olivine, and [START_REF] Palandri | A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modelingRep[END_REF] for enstatite, the lifetime of olivine and enstatite crystallites can be estimated to range between a few days and a few Myrs (Fig. 1). The longest durations are obtained for the larger enstatite grains, the more alkaline pH and low temperatures. As a matter of fact, these values argue for very short durations of serpentinization in the rocky core of Enceladus -and correspondingly, to short timescales of hydrogen generation through this process (see Fig. 1 and Eq. ( 4)).

C
Fig. 1. Estimated lifetime of olivine and enstatite grains at conditions relevant for fluid-rock interactions in Enceladus' core following the approach of [START_REF] Olsen | Using a mineral lifetime diagram to evaluate the persistence of olivine on Mars[END_REF]. Top: Olivine lifetime, at pH 9 (left) and pH 10 (right). Bottom: Olivine lifetime at pH 11 (left); Enstatite lifetime at pH 9-11 (the pH-dependence of enstatite dissolution in basic conditions is not known but assumed to be not significant in most studies dedicated to pyroxene dissolution).

In the next section, we explore how these rates are affected when more realistic fluid compositions are considered, so that the dependence of silicate dissolution rate on the Gibbs free energy of dissolution (∆G) can be accounted for.

Fluid composition dictated by an equilibrium with Talc + Serpentine (Chrysotile)

In order to evaluate whether the dissolution reaction takes place at far-from-equilibrium conditions, some constraints regarding the fluid composition are required to compute the Gibbs free energy of olivine and enstatite dissolution, which is the objective of the present section. Here we assume that the fluid composition is dictated by an equilibrium with talc (Tlc) and serpentine (chrysotile, Ctl), consistent with the first scenario described in [C Glein et al., 2018].

The chemical equation corresponding to chrysotile dissolution reads:

Ctl + 6H + = 2SiO 2 (aq) + 3Mg 2+ + 5H 2 O (7)
for which the associated equilibrium constant can be expressed following:

𝐾 𝐶𝑡𝑙 𝑇 = (𝑀𝑔 2+ ) 𝑒𝑞 3 . (𝑆𝑖𝑂 2 (𝑎𝑞)) 𝑒𝑞 2 (𝐻 + ) 𝑒𝑞 6 (8)
Similarly, the chemical equation corresponding to talc dissolution reads:

Tlc + 6H + = 4SiO 2 (aq) + 3Mg 2+ + 4H 2 O (9)
for which the associated equilibrium constant can be expressed following:

𝐾 𝑇𝑙𝑐 𝑇 = (𝑀𝑔 2+ ) 𝑒𝑞 3 . (𝑆𝑖𝑂 2 (𝑎𝑞)) 𝑒𝑞 4 (𝐻 + ) 𝑒𝑞 6 (10)
If one supposes that the fluid is at equilibrium at the same time with respect to talc and serpentine, then, Eq. ( 11) corresponding to (7) -( 9) yields:

Ctl + 2SiO 2 (aq) = H 2 O + Tlc (11) for which the associated equilibrium constant can be expressed following:

𝐾 𝐶𝑡𝑙 𝑇 𝐾 𝑇𝑙𝑐 𝑇 = 1 (𝑆𝑖𝑂 2 (𝑎𝑞)) 𝑒𝑞 2 (12a)
This imposes the following condition for the activity of SiO 2 (aq) (note that we are considering SiO2(aq) here, and not all dissolved Si species (ΣSi), which is very different at basic pH):

(𝑆𝑖𝑂 2 (𝑎𝑞)) 𝑒𝑞 = √ 𝐾 𝑇𝑙𝑐 𝑇 𝐾 𝐶𝑡𝑙 𝑇 (12b)
Substituting this value in e.g. Eq. ( 8) allows for the determination of the activity ratio between Mg 2+ and H + that verifies that the fluid is at equilibrium with respect to serpentine:

𝐾 𝐶𝑡𝑙 𝑇 2 𝐾 𝑇𝑙𝑐 𝑇 = (𝑀𝑔 2+ ) 𝑒𝑞 3 (𝐻 + ) 𝑒𝑞 6 (13a)
Obviously, this result would have been equivalent if Eq. ( 10) was used instead of Eq. ( 8), since the fluid is at equilibrium with respect to both chrysotile and talc. Therefore, (Mg 2+ ) eq can be determined at any fixed pH value after re-arranging Eq.( 13a):

(𝑀𝑔 2+ ) 𝑒𝑞 = 10 (2 3 ⁄ )log(𝐾 𝐶𝑡𝑙 𝑇 )-(1 3 ⁄ ) log(𝐾 𝑇𝑙𝑐 𝑇 )-2𝑝𝐻 (13b)
Finally, by combining Eqs. ( 12b) and (13b), one can demonstrate that if the solution is at equilibrium with respect to talc and chrysotile, this fixes the Gibbs free energy with respect to olivine dissolution (∆G Fo ) to a unique value, whatever the pH. This point is very important, as if one considers this mineral assemblage, there is no need to make any assumption about the pH of the fluid circulating in the rocky core:

∆𝐺 𝑓𝑜 = 𝑅𝑇. ln ( (𝑀𝑔 2+ ) 𝑒𝑞 2 . (𝑆𝑖𝑂 2 (𝑎𝑞)) 𝑒𝑞 𝐾 𝐹𝑜 (𝐻 + ) 𝑒𝑞 4 ) = 𝑅𝑇. ln (𝐾 𝐶𝑡𝑙 𝑇 5/6 . 𝐾 𝑇𝑐𝑙 𝑇 -1/6 . 𝐾 𝐹𝑜 𝑇 -1 ) (14)
The same holds true with respect to enstatite:

∆𝐺 𝑒𝑛𝑠 = 𝑅𝑇. ln ( (𝑀𝑔 2+ ) 𝑒𝑞 . (𝑆𝑖𝑂 2 (𝑎𝑞)) 𝑒𝑞 𝐾 𝑒𝑛𝑠 (𝐻 + ) 𝑒𝑞 2 ) = 𝑅𝑇. ln (𝐾 𝐶𝑡𝑙 𝑇 1/6 . 𝐾 𝑇𝑐𝑙 𝑇 1/6 . 𝐾 𝐸𝑛𝑠 𝑇 -1 ) (15)
The values of ∆G fo and ∆G ens were calculated following Eq. ( 14) and ( 15) using the thermodynamic equilibrium constants tabulated in the database of the CHESS code ([van der [START_REF] Van Der Lee | CHESS Tutorial and Cookbook[END_REF]). The values of ∆G range between -15.5 and -10.4 kJ/mol for enstatite, and between -32.5 and -18.5 kJ/mol for forsterite for temperatures ranging between 0 and 200 °C (Fig. 2a).

In the conventional framework of mineral dissolution kinetics, silicate dissolution rates (R silicate ) are supposed to obey the following TST-derived relation, which is implemented in classical reactive transport codes following (see [START_REF] Eyring | The activates complex in chemical reactions[END_REF]A.C. Lasaga, 1981]):

𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 = 𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 𝑝𝐻,𝑇 . (1 -exp ( ∆𝐺 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 𝑅𝑇 )) (16)
where 𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 𝑝𝐻,𝑇 corresponds to the dissolution rate constant of the considered silicate at given pH and T conditions. Practically, 𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 𝑝𝐻,𝑇 corresponds to the far-from-equilibrium dissolution rates determined in the previous section following Eq. ( 6). Using the values of ∆G fo and ∆G ens calculated using Eq. ( 14) and ( 15), the 𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 /𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 𝑝𝐻,𝑇 ratio can be calculated to estimate how the distance from equilibrium affects the dissolution rates for a fluid composition dictated by an equilibrium with talc and chrysotile. As can be seen in Fig. 2b, the 𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 / 𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 𝑝𝐻,𝑇 ratio remains close to 1 for all temperatures. In other words, with the assumptions detailed in the present section, the attenuation of the dissolution rate due to approaching equilibrium is negligible if the dissolution rate laws are supposed to obey the transition state theory. As can be seen, the dissolution rates remain within 90% of their far-from-equilibrium values.

Fluid composition dictated by an equilibrium with Talc + Calcite + Tremolite + Serpentine (Chrysotile)

In their study, [C Glein et al., 2018] suggested that the composition of aqueous fluids circulating in the core of Enceladus may be further refined based on constraints on the suspected chlorinity and carbonate concentration of the ocean waters. They suggested that downwelling ocean waters bring their chlorinity and carbonate concentration in the rocky core, and included calcite (Cal) and tremolite (Tr) as additional saturating phases. Below, we

show that such assumptions fix the activity of CO 2 (aq) and pH of the fluids. Then, we derive an expression for the activities of Mg 2+ and SiO 2 (aq), making it possible to constrain the fluid composition and ultimately, the Gibbs free energy of forsterite and enstatite dissolution.

First, [C Glein et al., 2018] proposed that the following equation buffers the activity of CO 2 (aq):

3Tr + 2Ctl + 6CO 2 (aq) = 6Cal + 7Tlc (17)
The corresponding equilibrium activity of CO 2 (aq) is therefore set by: Note that all equilibrium constants are directly available in the thermodynamic database of CHESS, except for K Cal , which can be calculated by combining the following equations:

(𝐶𝑂 2 (𝑎𝑞)) = 𝐾 𝐶𝑎𝑙 𝑇 . 𝐾 𝑇𝑙𝑐 𝑇 7/6 𝐾 𝑇𝑟 𝑇 1/2 . 𝐾 𝐶𝑡𝑙 𝑇 1/3 (18)
Cal + H + = Ca 2+ + HCO 3 - (19) CO 2 (aq) + H 2 O = H + + HCO 3 - ( 20
)
which yields:

Cal + 2H + = CO 2 (aq) + H 2 O + Ca 2+ (21)
The corresponding equilibrium constant can then be calculated from the equilibrium constants associated with Eq. ( 20) and ( 21), which were retrieved from the thermodynamic database of CHESS:

𝐾 𝐶𝑎𝑙 𝑇 = 𝐾 19 𝑇 𝐾 20 𝑇 ( 22
)
Substituting 𝐾 𝐶𝑎𝑙 𝑇 in Eq. ( 18) using Eq. ( 22) makes it possible to calculate the equilibrium activity of CO 2 (aq) for temperatures ranging between 0 and 200 °C (see Fig. 3a), which is required to ultimately calculate the pH of the downwelling waters. To reach this goal, the laws of mass action of the carbonate system can be combined, with a fixed constraint on the carbonate alkalinity. The law of mass action associated to Eq. ( 20) is:

𝐾 20 𝑇 = (𝐻 + ) 𝑒𝑞 . (𝐻𝐶𝑂 3 -) 𝑒𝑞 (𝐶𝑂 2 (𝑎𝑞)) 𝑒𝑞 (23)
The law of mass action associated to the protonation of carbonate species is:

𝐾 24 𝑇 = (𝐻𝐶𝑂 3 -) 𝑒𝑞 (𝐻 + ) 𝑒𝑞 . (𝐶𝑂 3 2-) 𝑒𝑞 (24)
Finally, [C Glein et al., 2018] recommended to use a total carbonate alkalinity of either 31 or 56 meq/kg H 2 O inherited from the ocean (C ocean ), depending on the ocean's pH:

[𝐻𝐶𝑂 3 -] + 2[𝐶𝑂 3 2-] = 𝐶 𝑜𝑐𝑒𝑎𝑛 (25)
Combining Eqs. (23-25) makes it possible to derive an equation to estimate [HCO 3 -]:

2𝛾 𝐻𝐶𝑂3 2 𝐾 20 𝑇 . 𝐾 24 𝑇 . 𝛾 𝐶𝑂3 . (𝐶𝑂 2 (𝑎𝑞)) 𝑒𝑞 [𝐻𝐶𝑂 3 -] 2 + [𝐻𝐶𝑂 3 -] -𝐶 𝑜𝑐𝑒𝑎𝑛 = 0 (26)
For which the solution is: For temperatures ranging between 0 and 200 °C, the values of [HCO3 -] that satisfy this equation are shown in Fig. 3b, supposing an ionic strength fixed to 0.13 molal (here corresponding to [Cl -] = 0.1 M and carbonate alkalinity of 0.03 M).

[𝐻𝐶𝑂 3 -] = (
Finally, re-arranging Eq. ( 23) allows for the calculation of pH following Eq. ( 28), where the values of (CO 2 (aq)) and [HCO 3 -] were calculated using respectively Eq. ( 18) and Eq. ( 27)

(see Fig. 3a-b):

pH = log(𝛾 𝐻 . [𝐻𝐶𝑂 3 -]) -log 𝐾 20 𝑇 -log((𝐶𝑂 2 (𝑎𝑞)) 𝑒𝑞 ) (28)
As can be seen in Fig. (3c), the pH of the circulating fluids ranges between 12.0 and 8.7 for temperatures ranging between 0 and 200 °C.

Finally, because the fluid is in equilibrium with talc and chrysotile, the equilibrium activity of SiO 2 (aq) is simply set by Eq. ( 12b), and that of Mg 2+ can ultimately be calculated by re-arranging Eq. ( 10): Overall, supposing that the fluid chemistry is dictated by pre-determined saturating phases typical of serpentinization reactions (as suggested in, e.g., [C Glein et al., 2018;[START_REF] Zandanel | Short lifespans of serpentinization in the rocky core of Enceladus: Implications for hydrogen production[END_REF] and references therein) virtually corresponds to far-from-equilibrium dissolution conditions if one considers the rate-∆G relation resulting from the transition state theory, which is implemented in most reactive transport codes. As a consequence, and in spite of the simplistic calculations conducted in Section 2, the conclusions remain unchanged, and such calculations already provide a rather good picture of the mineral reactivity modeled using classical reactive transport codes simulating permanent fluid-rock contact. This results directly stems from the implementation of TST-based rate laws in reactive transport models.

(𝑀𝑔 2+ ) 𝑒𝑞 = ( 𝐾 𝑇𝑙𝑐 𝑇 . (𝐻 + ) 𝑒𝑞
In the next section, we explore the extent to which alternative rate laws that were not [START_REF] Maher | Rates of silicate dissolution in deep-sea sediment: In situ measurement using U-234/U-238 of pore fluids[END_REF][START_REF] Paces | Rate Constants of Dissolution Derived from the Measurements of Mass Balance in Hydrological Catchments[END_REF][START_REF] Velbel | Constancy of silicate-mineral weathering-rate ratios between natural and experimental weathering: implications for hydrologic control of differences in absolute rates[END_REF][START_REF] White | The effect of time on the weathering of silicate minerals: why do weathering rates differ in the laboratory and field?[END_REF]Art F. White et al., 2017] and references therein). As a consequence, one can wonder whether the short lifetimes determined above could only be another illustration of such a well-known discrepancy existing between field-derived and laboratory-derived estimates.

The pioneering work of [Antonio C. [START_REF] Lasaga | Surface chemistry, etch pits and mineral-water reactions[END_REF] expanded later by [Antonio C. [START_REF] Lasaga | Variation of Crystal Dissolution Rate Based on a Dissolution Stepwave Model[END_REF] suggested that the nucleation of etch pits at screw dislocation outcrop could play a critical role in the kinetics of mineral dissolution, resulting in rate-∆G relations that radically deviate from TST-based rate laws (Eq. ( 16)). Several studies later demonstrated that the use of such alternative dissolution rate laws decrease modeled mineral dissolution rates, which ultimately results in a better agreement between the outputs of simulations and field estimates (e.g., [K. [START_REF] Maher | The role of reaction affinity and secondary minerals in regulating chemical weathering rates at the Santa Cruz Soil Chronosequence, California[END_REF]Jacques Schott et al., 2012;[START_REF] Wild | In-situ dissolution rates of silicate minerals and associated bacterial communities in the critical zone (Strengbach catchment, France)[END_REF] and references therein). The present section aims at evaluating the consequence of the use of such an alternative conceptual framework.

The model developed by [Antonio C. Lasaga and Luttge, 2001] is based on the following concept: at far-from-equilibrium conditions, the nucleation of etch pits at line defects surface is spontaneous. The pit walls then represent a continuous source of steps which travel throughout the crystal surface, eventually leading to its overall retreat (dissolution).

Conversely, if the energetic threshold to nucleate etch pits (often referred to as the critical free energy of etch pit nucleation, or ∆G crit ) is exceeded, etch pit nucleation at dislocation outcrops is no longer spontaneous, and the dissolution is essentially driven by point defects and dissolution on pre-existing corners and edges at a much slower rate. This results in two very distinct dissolution regimes, for which the dissolution rate at ∆G < ∆G crit may be up to two to three orders of magnitude greater than the rate at ∆G > ∆G crit . The direct observation of a dramatic drop of the dissolution rate associated to mineral surfaces free of etch pits has been well documented for many minerals including albite ( [START_REF] Arvidson | Mineral dissolution kinetics as a function of distance from equilibrium -New experimental results[END_REF][START_REF] Beig | Albite dissolution kinetics as a function of distance from equilibrium: Implications for natural feldspar weathering[END_REF]), K-spar ([Pollet-Villard et al., 2016]) and calcite ( [START_REF] Bouissonnié | From mixed flow reactor to column experiments and modeling: Upscaling of calcite dissolution rate[END_REF]), and may be universal, although to the best of our knowledge, it has not been verified directly for olivine and enstatite. However, [Antonio C. [START_REF] Lasaga | Surface chemistry, etch pits and mineral-water reactions[END_REF] provided a theoretical relation that can be used to calculate ∆G crit for etch pit nucleating at screw dislocation for any mineral:

∆𝐺 𝑐𝑟𝑖𝑡 = -2𝑉 ̅ 𝜋 2 𝜎 𝑠 2 𝜇𝑏 2 (30)
where σ the surface energy is not unique, and depends on the crystallographic face that is considered.

On the same footing, several Burgers vectors exist for a given crystal structure. As a consequence we will only consider here, as a first approximation, the combination of b and σ values which result in the lowest ∆G crit estimates reported in [Antonio C. [START_REF] Lasaga | Surface chemistry, etch pits and mineral-water reactions[END_REF], ensuring that we do not artificially force olivine and enstatite reactivity to occur in the fastest dissolution regime. As shown in [Antonio C. [START_REF] Lasaga | Surface chemistry, etch pits and mineral-water reactions[END_REF], the lowest expected value of ∆G crit is -12.0 kJ/mol for forsterite, and -7.5 kJ/mol for enstatite. Such values remain above the Gibbs free energy values determined in Section 3 (see Fig. 2). As a consequence, etch pit nucleation should remain spontaneous in the conditions of fluid-rock interactions in the Enceladus' core, corresponding to the far-from-equilibrium dissolution regime already described in Sections 2 and 3. Stated in other words, considering that the reactivity of olivine and enstatite is driven by etch pit nucleation does not modify the expected lifetime of these minerals and the associated timespan of hydrogen generation if one supposes permanent fluid-rock contact in Enceladus' core.

Reactivity of olivine and enstatite controlled by amorphous silica rich surface layers (ASSLs)

As mentioned above, knowledge of a systematic kinetic 'field-lab discrepancy' has fostered the search for alternative reaction mechanisms and corresponding rate laws that deviate from the simple TST-based rate equations. More specifically, it has been previously

shown that the surface of weathered olivines and pyroxenes spontaneously evolves towards a configuration with a lesser reactivity. This process is often referred to as "surface aging" [START_REF] Eggleston | Sample Preparation and Aging Effects on the Dissolution Rate and Surface-Composition of Diopside[END_REF]. The formation of nm-thick amorphous layers enriched in silica that cover the pristine surface of the dissolving substrate as a result of their incongruent dissolution represents one of such aging mechanisms. Such layers are known to have a much lower reactivity than the underlying surface, and the evolution of their transport properties may gradually isolate the dissolving pristine substrate. This mechanism has been clearly demonstrated for olivine [START_REF] Bearat | Carbon sequestration via aqueous olivine mineral carbonation: Role of passivating layer formation[END_REF][START_REF] Daval | Influence of amorphous silica layer formation on the dissolution rate of olivine at 90 degrees C and elevated pCO(2)[END_REF][START_REF] Johnson | Olivine dissolution and carbonation under conditions relevant for in situ carbon storage[END_REF]Kate Maher et al., 2016]) and pyroxenes [START_REF] Daval | Linking nm-scale measurements of the anisotropy of silicate surface reactivity to macroscopic dissolution rate laws: New insights based on diopside[END_REF][START_REF] Schott | X-Ray Photoelectron Studies of the Mechanism of Iron Silicate Dissolution during Weathering[END_REF]) in the acidic pH range. Whether such a mechanism also applies in the basic pH range remains an open question. Whereas the formation of ASSLs in the basic pH range has been well documented for the weathering of borosilicate glasses (e.g., [Stephane Gin et al., 2020;[START_REF] Gin | Atom-Probe Tomography, TEM and ToF-SIMS study of borosilicate glass alteration rim: A multiscale approach to investigating rate-limiting mechanisms[END_REF]), very few studies have investigated the chemical composition at the nm-scale of the near-surface region of dissolved olivine and pyroxenes in the basic pH range. Worthy of note, the spectroscopic characterizations of the surface composition of altered olivine performed by [START_REF] Pokrovsky | Forsterite surface composition in aqueous solutions: A combined potentiometric, electrokinetic, and spectroscopic approach[END_REF] suggested that ~2 nm-thick ASSLs could be evidenced up to a pH of ~ 9. Therefore, we believe that determining the lifetime of olivine and enstatite for a dissolution regime controlled by the formation of ASSLs may be worth considering.

In appearance, for silicate minerals passivated by ASSLs, the driving force of the reaction is determined by the thermodynamic properties of ASSLs, such that the chemical affinity that should be considered for the dissolution reaction has to be defined with respect to the surface layers themselves. This results in the following rate equation (e.g. [START_REF] Daval | Linking nm-scale measurements of the anisotropy of silicate surface reactivity to macroscopic dissolution rate laws: New insights based on diopside[END_REF][START_REF] Daval | Influence of amorphous silica layer formation on the dissolution rate of olivine at 90 degrees C and elevated pCO(2)[END_REF]):

𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 = 𝑘 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 𝑝𝐻,𝑇 . (1 -exp ( ∆𝐺 𝑆𝑖𝑂2 𝑅𝑇 )) (31)
where ∆G SiO2 represents the Gibbs free energy with respect to a Si-rich surface layer with an apparent solubility close to SiO 2 (am) for olivine [START_REF] Daval | Influence of amorphous silica layer formation on the dissolution rate of olivine at 90 degrees C and elevated pCO(2)[END_REF]), or cristobalite-α for pyroxene [START_REF] Daval | Linking nm-scale measurements of the anisotropy of silicate surface reactivity to macroscopic dissolution rate laws: New insights based on diopside[END_REF]). Following the calculations detailed in Section 3, the activity of SiO 2 (aq) is fixed by the equilibrium between talc and chrysotile. Then Eq. (12b) makes it possible to calculate ∆𝐺 𝑆𝑖𝑂2 at any T following:

∆𝐺 𝑆𝑖𝑂2 = 𝑅𝑇. ln ( √𝐾 𝑇𝑙𝑐 𝑇 𝐾 𝐶𝑡𝑙 𝑇 ⁄ 𝐾 𝑆𝑖𝑂2 𝑇 ) (32)
where 𝐾 𝑆𝑖𝑂2 𝑇 stands for the solubility constant of the ASSL, which may range between that of SiO 2 (am) and cristobalite-α. Using the value of ∆𝐺 𝑆𝑖𝑂2 calculated with Eq. ( 32) and shown in Fig. 4a, the 𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 /𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 𝑝𝐻,𝑇 ratio was calculated to estimate how accounting for the formation of ASSLs may impact the dissolution rate. As can be seen in Fig. 4b, the 𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 / 𝑅 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 𝑝𝐻,𝑇 ratio remains close to 1 for all temperatures. Stated in other words, even under the hypothesis that ASSLs control olivine and pyroxenes dissolution at pH up to 11-12, the solution composition corresponds to fluid highly undersaturated with respect to ASSLs if one assumes that the fluid is buffered by the mineral assemblages proposed by [C Glein et al., 2018]. As a consequence, the mechanisms explored in the present section lead to lifetimes of enstatite and olivine as short as those estimated in Sections 2 and 3.

Transport-controlled serpentinization timespans

In the previous sections, it was assumed that the rate-limiting step of the serpentinization reaction was the dissolution of the considered silicate, supposed to be surface-controlled.

Although it has been suggested for over 40 years that the weathering of silicates on Earth was actually surface-controlled (e.g., [START_REF] Berner | Dissolution of pyroxenes and amphiboles during weathering[END_REF]), we propose here to evaluate the timespans of serpentinization for reaction regimes that would be hypothetically transportlimited, in order to see whether this would significantly extend the expected timespan of serpentinization.

Dissolution rates controlled by sluggish diffusion in nanoporous media

There is significant evidence to suggest that replacement reactionsto which serpentinization reactions may be assimilatedoccur through an interface-coupled dissolution-precipitation mechanism rather than solid-state diffusion ( [START_REF] Plumper | The interface-scale mechanism of reaction-induced fracturing during serpentinization[END_REF][START_REF] Putnis | Mineral Replacement Reactions[END_REF]). Depending on whether the reaction is isovolumetric or not, serpentinization may either result in pseudomorphic replacement or grain fracturing. Whatever the case, a permanent contact between the bulk fluid and the primary silicate interface must be maintained for the reaction to proceed. Several studies suggest that this mass transfer may occur through the nanoporosity generated in the reaction products [START_REF] Cailleteau | Insight into silicate-glass corrosion mechanisms[END_REF][START_REF] Daval | Dynamics of altered surface layer formation on dissolving silicates[END_REF][START_REF] Daval | Mechanism of wollastonite carbonation deduced from micro-to nanometer length scale observations[END_REF][START_REF] Plümper | Fluid escape from subduction zones controlled by channel-forming reactive porosity[END_REF]). The diffusivity of such secondary assemblages is very poorly constrained, and the transport of dissolved species might even not follow simple Fickian diffusion ([Stéphane Gin et al., 2018;[START_REF] Plümper | Fluid escape from subduction zones controlled by channel-forming reactive porosity[END_REF]). The few studies that were intended to characterize nanoporous media in terms of their diffusivity argue for a large range of possible diffusion coefficient values. For instance, molecular dynamics simulations performed by [START_REF] Bourg | Molecular Dynamics Simulations of Water Structure and Diffusion in Silica Nanopores[END_REF] suggest that as long as the pore dimensions remain larger than 2 nm, water keeps its bulklike properties, and the diffusion coefficients of aqueous species in nanopores remain similar to their counterparts in bulk water (~2×10 -5 cm²/s at 25 °C). Conversely, direct attempts to estimate the diffusivity of nm-to µm-thick secondary coatings through the use of conservative tracers showed that it could range between 10 -10 to 10 -17 cm²/s at 25 °C [START_REF] Daval | Dynamics of altered surface layer formation on dissolving silicates[END_REF].

[ [START_REF] Daval | Influence of amorphous silica layer formation on the dissolution rate of olivine at 90 degrees C and elevated pCO(2)[END_REF] suggested that the diffusivity of surface layers formed on olivine may be as low as 5×10 -18 cm²/s at 90 °C. This latter value will be used to run the calculations below, although it may be considered as a lower bound, as it corresponds to the diffusivity in the putative porosity of an ASSL, which is unlikely to be stable in the conditions relevant to the serpentinization of Enceladus' core (see Section 5). More likely, the secondary coatings should be made of an assemblage of nm-to µm-scale crystallites, whose diffusivity was previously demonstrated to be much higher than that of ASSLs (see [START_REF] Saldi | Mineralogical evolution of Fe-Si-rich layers at the olivine-water interface during carbonation reactions[END_REF] for details).

The diffusion-controlled lifetime of a grain of olivine or enstatite can be approximated following:

𝑡 ≈ 𝑟 𝑔𝑟𝑎𝑖𝑛 ² 𝐷 𝑇 = 𝑟 𝑔𝑟𝑎𝑖𝑛 ² 𝐷 0 . exp (-𝐸 𝑎,𝑑𝑖𝑓𝑓 𝑅𝑇) ⁄ (33)
where D T [cm²/s] is the diffusion coefficient of aqueous species in the nanoporous medium at the considered temperature, D 0 [cm²/s] is the intrinsic diffusion constant, and 𝐸 𝑎,𝑑𝑖𝑓𝑓 is the corresponding activation energy. To the best of our knowledge, the value of this latter parameter has never been determined in previous studies, but should be similar to that of water self-diffusion (~10 kJ/mol). In their study, [START_REF] Rotenberg | Water and ions in clays: Unraveling the interlayer/micropore exchange using molecular dynamics[END_REF] calculated that the activation energy of the exchange of water and ions between clay interlayers and microporosity (≤ 2 nm) ranged between ~0 and ~23 kJ/mol for cations and anions, respectively. Considering this process as a reasonable analogue of the migrations of ions in the clayey coatings covering olivine grains, we set 𝐸 𝑎,𝑑𝑖𝑓𝑓 = 15 kJ/mol in the calculations.

Considering grain sizes comprised between 10 µm and 1 mm in radius (see Section 2) results in lifetimes ranging between ~ 2000 years and 326 Myr (Fig. 5). Such estimates are about two orders of magnitude larger than the longest durations derived from Section 2 assuming surface-controlled rates, although they remain relatively short at geological timescales (and likely overestimate the actual duration of the process, as lower bounds of diffusion coefficients were considered).

Serpentinization rates controlled by slow fluid velocities in Enceladus' core

Serpentinization reactions are known to consume water. This aspect may become ratelimiting if the dissolution of olivine in the core occurs faster than the supply of water from the ocean, which is a likely configuration, as fluid velocity may be as low as 1 cm/yr in the core ( [START_REF] Choblet | Powering prolonged hydrothermal activity inside Enceladus[END_REF]). To test this scenario, let us imagine that (i) downwelling waters are gradually percolating through the initially fluid-saturated porosity of the core and (ii) serpentinization occurs as soon as water contacts olivine (consistent with the assumption that serpentinization is fast compared to water supply). In the following, we consider that the rocky core is initially entirely made up of fluid-saturated pores and olivine.

As shown in Section 2, the lifetime of olivine in a water-saturated environment never exceeds 200 kyrs at conditions relevant for fluid-rock interactions in the Enceladus' core (see Fig. 1). After the consumption of the initial amount of water that was filling up the pores, the proportion of olivine in the core that will remain to be altered through water circulation corresponds to 1 -𝛷 -𝑄 𝑑 𝛷, where 𝛷 is the porosity and Q d is a parameter that accounts for the water to olivine molar ratio requested to complete the serpentinization reaction (Y), as well as the discrepancies between the molar volumes of olivine and water. This parameter can therefore be expressed following:

𝑄 𝑑 = 𝑉 ̅ 𝑜𝑙 𝑉 ̅ 𝑤𝑎𝑡𝑒𝑟 . 𝑌 (34)
In the above equation, 𝑉 ̅ 𝑜𝑙 and 𝑉 ̅ 𝑤𝑎𝑡𝑒𝑟 are the molar volumes of olivine (45 cm 3 /mol) and water (18 cm 3 /mol), respectively. The water to olivine molar ratio requested to complete the reaction is given by the stoichiometric coefficient of water in the corresponding serpentinization reaction Y (e.g., (2-2x/3) in Eq. ( 3)). Several values can be found in the literature for Y (e.g. [START_REF] Malvoisin | Serpentinization of oceanic peridotites: 1. A high-sensitivity method to monitor magnetite production in hydrothermal experiments[END_REF] and references therein) depending on the nature of the secondary phases and iron content in olivine, and are tested below (Fig. 6).

For a fully connected porous medium, the longest lifetime of olivine corresponds to the time required for water to reach the center of the core. This duration is proportional to the ratio between the remaining volume of olivine to be serpentinized before fluid circulation to the volume of pores, accounting for the Q d factor described above:

𝑡 = 𝑟 𝑐𝑜𝑟𝑒 𝜐 ( 1 -𝛷 -𝑄 𝑑 𝛷 𝑄 𝑑 𝛷 ) = 𝑟 𝑐𝑜𝑟𝑒 𝜐 ( 𝑉 ̅ 𝑤𝑎𝑡𝑒𝑟 . 𝑌. (1 -𝛷) 𝑉 ̅ 𝑜𝑙 𝛷 -1) (35)
where 𝑟 𝑐𝑜𝑟𝑒 [cm] is the radius of the core (1.9×10 7 cm) and 𝜐 [cm/yr] the fluid velocity. We stress that, in the above equation, the temporal decrease of porosity due to the positive difference of molar volumes between the solid products and reactants, which can range between ~ +20% and ~ +50%, is neglected. Instead, we calculated t for 𝛷 ranging between its upper bound (~ 30%) and its lower bound (~ 2%), the latter corresponding to the residual porosity expected after the full serpentinization of the core using an initial porosity of ~30% and the prescribed value of volume increase of ~ +40% proposed in the recent work of [START_REF] Klein | Quantifying the volume increase and chemical exchange during serpentinization[END_REF]. Note that without the presence of organics, such a low value of porosity is difficult to reconcile with the estimated density of the satellite's core [START_REF] Mckinnon | Effect of Enceladus's rapid synchronous spin on interpretation of Cassini gravity[END_REF].

Finally, [START_REF] Mason | Olivine composition in chondrites[END_REF] showed that the value of x for olivine in chondrites can range from ~0.28 to ~0.62. Considering the slowest fluid velocity suggested by [START_REF] Choblet | Powering prolonged hydrothermal activity inside Enceladus[END_REF], the lifetime of serpentinization ranges between < 1 Myr and ~ 500 Myr, depending on the porosity, the Fe content in olivine and chemical reaction considered (Fig. 6). Note that these values are in excellent agreement with the outputs of 1D-reactive transport simulations of [START_REF] Zandanel | Short lifespans of serpentinization in the rocky core of Enceladus: Implications for hydrogen production[END_REF] conducted with a similar set of input parameters for fluid flow rate and core porosity. Therefore, this result independently confirms that the rate-limiting step in their corresponding simulations is the water supply, consistent with the sharp serpentinization fronts reported in their Fig. 3. Importantly, such results directly stem from the intrinsically high dissolution rate of olivine (relative to fluid velocity), which itself is a consequence of a fluid composition that is locked at far-from-equilibrium conditions (see Sections 3 and 4).

Fig. 6. Lifetime of olivine (and therefore, production of molecular H 2 ) for a process controlled by slow fluid velocity in Enceladus' core. The x parameter refers to the Fe content in olivine as reported in Eq. ( 1). Two different olivine alteration reactions have been tested, namely (i) a reaction where the secondary phases are Fefree serpentine, Fe-free brucite, magnetite and SiO 2 (am): Ol + (1.5-1.25x/3)H 2 O = (0.5-0.25x) Sp + (0.5-0.25x)

Brc + (x/3) Mt + (x/2) SiO 2 (am) + (x/3)H 2 ; (ii) Eq. ( 2) from [START_REF] Malvoisin | Serpentinization of oceanic peridotites: 1. A high-sensitivity method to monitor magnetite production in hydrothermal experiments[END_REF], where the secondary phases are supposed to be Fe-free serpentine, Fe-free brucite and magnetite, i.e.: Ol + (1.5-2x/3) H 2 O = 0.5 Sp + (0.5-x)

Brc + (x/3) Mt + (x/3) H 2 . Note that if reactions involving the serpentine-talc buffer are considered, e.g.: Ol + (1.5-2x/3) H 2 O + 0.5 Sp = (2-x) Mg 2+ + (4-2x) OH -+ 0.5 Tlc + (x/3) Mt + (x/3) H 2 , the stoichiometric coefficient of water is the same as the one found in in the equation proposed by [START_REF] Malvoisin | Serpentinization of oceanic peridotites: 1. A high-sensitivity method to monitor magnetite production in hydrothermal experiments[END_REF], so that the estimated olivine lifetime remains unchanged. Overall, considering a serpentinization regime that is controlled either by the water velocity in the core of Enceladus or by aqueous diffusion in nanoporous secondary minerals provides serpentinization timespans that are at least two orders or magnitude greater than those estimated for surface-controlled reactions. Nonetheless, these results suggest that hydrothermal activity in the core developed only recently and no earlier than 500 Myr ago, in agreement with the results of [START_REF] Zandanel | Short lifespans of serpentinization in the rocky core of Enceladus: Implications for hydrogen production[END_REF].

Using Cassini spacecraft data to infer the onset and timescales of serpentinization on Enceladus

The present section aims at evaluating the rate parameters that would be required to match the measurements of the Cassini spacecraft, supposing that olivine dissolution is either (i) surface-controlled (section 7.1) or (ii) diffusion-controlled (section 7.2). For each scenario, we briefly discuss the physical acceptability of the range of values resulting from such a treatment.

Effective size of silicate aggregates in Enceladus' core

In the present section, we estimate the effective size that individual aggregates of Febearing silicates would have to match the flux of hydrogen (F H2 ; [mol/s]) measured from the plume gas by the Cassini spacecraft (cf. [START_REF] Waite | Cassini finds molecular hydrogen in the Enceladus plume: evidence for hydrothermal processes[END_REF]). In that respect, a shrinking sphere model of the aggregates that make up Enceladus' core is described below. It is supposed that the core is made of a close-packing of equal spherical and non-porous aggregates. Therefore, all porosity is located outside the aggregates, and amounts 26% due to the specific geometric configuration that has been chosen. We suppose that that the dissolution of primary Fe-bearing silicates is the rate-limiting step of H 2 generation.

First, for given physicochemical conditions, the radial dissolution rate [cm/s] of an aggregate can be derived from Eq. ( 5):

𝑑𝑟 𝑎𝑔𝑔 𝑑𝑡 = -𝑉 𝑖 ̅ . 𝑅 𝑖 (36)
The negative sign of the derivative ensures that the aggregates shrink with time. Integrating Eq. ( 36) makes it possible to calculate the evolution the radius of a mineral aggregate as a function of time:

𝑟(𝑡) 𝑎𝑔𝑔 = 𝑟(𝑡 = 0) 𝑎𝑔𝑔 -𝑉 ̅ 𝑖 . 𝑅 𝑖 . 𝑡 (37)
The surface area of each individual aggregate is given by:

𝑆(𝑡) 𝑎𝑔𝑔 = 4𝜋. 𝑟(𝑡) 𝑎𝑔𝑔 2 (38)
Recalling that the volume of Enceladus'core can be written following:

𝑉 0 = (4𝜋 3 ⁄ ). 𝑟 𝑐𝑜𝑟𝑒 3 (39)
This parameter can be related to the volume of individual aggregates (V agg ) that make up the core following:

𝛷 = 𝑉 0 -𝑁 𝑎𝑔𝑔 𝑉 𝑎𝑔𝑔 𝑉 0 = 1 -𝑁 𝑎𝑔𝑔 𝑟(𝑡 = 0) 𝑎𝑔𝑔 3 𝑟 0 3 ( 40
)
where Φ is the porosity, N agg is the number of aggregates, and r agg their radius. Therefore, the total number of aggregates can be calculated following:

𝑁 𝑎𝑔𝑔 = (1 -𝛷) 𝑟 0 3 𝑟(𝑡 = 0) 𝑎𝑔𝑔 3 (41)
Expressing the total surface area of the aggregates (S TOT ) as a function of the surface area of individual spherical aggregates and replacing N agg by its expression taken from Eq. ( 41) yields:

𝑆(𝑡) 𝑇𝑂𝑇 = 𝑁 𝑎𝑔𝑔 . 4𝜋. 𝑟(𝑡) 𝑎𝑔𝑔 2 = (1 -𝛷). 4𝜋. 𝑟 0 3 𝑟(𝑡 = 0) 𝑎𝑔𝑔 3 . 𝑟(𝑡) 𝑎𝑔𝑔 2 (42)
If the core is now decomposed into n concentric spherical shells, then the surface area of the aggregates located in each infinitesimal shell j can be expressed following:

𝑆(𝑡) 𝑇𝑂𝑇,𝑗 = (1 -𝛷) 𝑟(𝑡 = 0) 𝑎𝑔𝑔 3 (𝑟 𝑗-1 3 -𝑟 𝑗 3 )4𝜋. 𝑟(𝑡) 𝑎𝑔𝑔
Finally, if one supposes that the flux of hydrogen (F(t) H2 ; [mol/s]) measured from the plume gas by the Cassini spacecraft chiefly results from the chemical equations ( 3) or (4), the flux of H 2 resulting from the release and oxidation of Fe 2+ can be written following:

𝐹(𝑡) 𝐻2 = (𝑥 3 ⁄ ). ∑ 𝑅 𝑖,𝑗 . 𝑆(𝑡) 𝑇𝑂𝑇,𝑗 𝑛 𝑗=1 (44)
where 𝑅 𝑖,𝑗 is the dissolution rate at the temperature in the j th infinitesimal shell. It is supposed that the temperature of the core observes a linear gradient ranging between 0 °C at the interface with the ocean and 200 °C at the center of the core. Combing Eq. ( 43) with Eq. ( 44) yields:

𝐹(𝑡) 𝐻2 = 4𝜋. (1 -𝛷). (𝑥 3 ⁄ ) 𝑟(𝑡 = 0) 𝑎𝑔𝑔 3 ∑ 𝑅 𝑖,𝑗 . (𝑟 𝑗-1 3 -𝑟 𝑗 3 ). 𝑛 𝑗=1 (𝑟(𝑡 = 0) 𝑎𝑔𝑔 -𝑉 ̅ 𝑖 . 𝑅 𝑖,𝑗 . 𝑡) 2 (45)
Practically, the shrinking rate of the aggregate will be controlled by the dissolution rate of the mineral with the slowest dissolution rate, i.e., orthopyroxene (see Section 2). The calculations below will therefore be conducted using R enstatite as the dissolution rate, supposing a composition of En 50 (x = 0.5). We suppose that the porosity of the core at the formation of Enceladus is identical to today's porosity (i.e., 26%). This value may be considered as an upper bound, since the variation of molar volume associated to serpentinization is usually positive. Worthy of note, changing the initial porosity would not drastically change the results of simulations. The evolution of F H2 is reported in Fig. 7 for timescales of H 2 production that are compatible with the estimations of today's flux measured by the Cassini spacecraft and the age suspected for Enceladus (here, between 100 Myr and 4.5 Gyr, the former value corresponding to the controversial young age of formation proposed by e.g. [Ćuk et al., 2016]). Such constraints require the aggregates to be on the order of 1 cm in radius for an age of 100 Myr, and between 1 and 5 m for an age of 4.5 Gyr. Very importantly, aggregate sizes on the order of decimeters or even more do not make any sense in classical reactive transport simulations. Usually, the grain size is supposed to be close to that of crystallites, i.e., up to one millimeter or so in our case. Assuming that the aggregate size is within the abovementioned range actually refers to fluid flow in a fractured medium, and not in a homogeneous porous medium. In addition, such estimates are at odds with the aggregate size required to support the tidal heating model developed in [START_REF] Choblet | Powering prolonged hydrothermal activity inside Enceladus[END_REF]. One could also conceive that the hydrothermal activity presently suspected in Enceladus is much more recent than the moon itself: more reasonable aggregate sizes provide lifetimes of about 1 Myr (0.1 mm) or 10 Myr (1 mm). We note however that the slopes of the observed curves for such grain sizes are extremely steep in the region corresponding to the value of F H2 expected for Enceladus at present, which would imply the unlikely notion that the Cassini spacecraft witnessed an exceptional instant for the moon.

As a consequence, and in agreement with all of the scenarios mentioned above, (i) either hydrothermal circulation and subsequent core serpentinization are very recent and the present time is a very unusual period for Enceladus' activity or (ii) serpentinization is not the main source of H 2 on Enceladus.

Effective diffusivity of surface layers developed on silicate grains

Here we suppose that the rate of H 2 generation is controlled by the slow diffusion of reactants across the nanoporosity of coatings covering the primary silicates that make up Enceladus core (cf. Section 6.1 for details). The present section aims at estimating the diffusion coefficient that would be required to match today's H 2 flux measured by the Cassini spacecraft for serpentinization durations compatible with the age of Enceladus.

Following the abovementioned assumptions, the evolution of the grain radius can be approximated following:

𝑟(𝑡) 𝑔𝑟𝑎𝑖𝑛,𝑗 = 𝑟 (0)𝑔𝑟𝑎𝑖𝑛,𝑗 -√𝐷 0 . exp ( -𝐸 𝑎,𝑑𝑖𝑓𝑓 𝑅𝑇 𝑗 ) . 𝑡 (46)
where the subscript j indicates that the considered grains j are located at a depth where the temperature T j prevails. Over a given time interval ∆t, the flux of H 2 resulting from the dissolution of one single grain can therefore be written following: If the core is sliced into n concentric spherical shells (see Section 7.1), then the total flux of H 2 generated in the j th increment (located at a distance from the center of the core comprised between r j and r j-1 ) can be calculated by combining Eq. ( 41) and (47):

𝐹(𝑡) 𝐻2,𝑗 = 𝑁 𝑔𝑟𝑎𝑖𝑛𝑠,𝑗 𝐹(𝑡) 𝐻2,𝑔𝑟𝑎𝑖𝑛,𝑗 = (𝑟 𝑗 3 -𝑟 𝑗-1 3 ) (1 -𝛷) 𝑟(0) 𝑔𝑟𝑎𝑖𝑛 3 𝐹(𝑡) 𝐻2,𝑔𝑟𝑎𝑖𝑛,𝑗 (48)
Finally, the total flux of H 2 coming from the whole core over a given time interval ∆t can be simply written following:

𝐹(𝑡) 𝐻2 = ∑ 𝐹(𝑡) 𝐻2,𝑗 𝑛 𝑗=1 (49)
As can be seen in Fig. 8, considering an activation energy of 15 kJ/mol (see Section 6.1)

and the maximum grain size of 1 mm to allow for the longest lifetime of primary silicate grains, the value of the intrinsic diffusion constant reported in Section 6.1 is too high to support serpentinization durations of several Gyrs. Actually, we found that D 0 would have to be decreased by two orders of magnitude to get H 2 fluxes compatible with that measured by the Cassini spacecraft and serpentinization durations of several Gyr. Although such low values of diffusivities seem fairly unlikely (for the record, the diffusion coefficient selected in Section 6.1 already falls in the lower range of known diffusion coefficients in nanoporous media), knowledge of the transport properties of secondary coatings formed through an interface-coupled dissolution-precipitation mechanism is still in its infancy, and there is no obvious theoretical limitations that would prevent apparent diffusivities lower than currently known in the literature. As a consequence, we suggest that future studies should be directed to ascertain whether under the conditions relevant for Enceladus (high pH, temperature range comprised between 0 °C and 200 °C), the alteration of primary Fe-bearing silicates could result in the formation of secondary coatings with unexpectedly low diffusivity, making olivine serpentinization a reliable mechanism for the observed H 2 flux by the Cassini spacecraft.

Concluding remarks

The observation of H 2 in Enceladus plume poses an interesting problem to understand the formation and evolution of this tiny moon of Saturn. In the framework of serpentinization reactions, various rate-limiting mechanisms have been explored in the present paper, resulting in the characteristic timescales shown in Fig. 9. As can be seen in Fig. 9, the lifetime of H 2 generation following the commonly invoked rate-limiting steps of Fe-bearing silicate serpentinization ranges between a few millions and a few hundreds of millions years. Outlier scenarios explored in Sections 7.1 and 7.2 include (i) a possible characteristic aggregate size of up to a few meters in the core, and (ii) the very slow diffusion of reactants across nanoporous secondary coatings covering individual olivine (or enstatite) grains. Regarding (i), this scenario would actually result in weathering textures similar to those observed on partially serpentinized ophiolites on Earth. However, it is not supported by any observation of carbonaceous chondrites or other primitive silicate meteorites considered to be building blocks of icy moons such as Enceladus, and it would be hardly compatible with the grain size required to support the tidal heating mechanism

proposed by [START_REF] Choblet | Powering prolonged hydrothermal activity inside Enceladus[END_REF] that is suspected to occur on Enceladus. Regarding (ii), current knowledge on the characteristic diffusivity of secondary coatings formed on Fe- bearing silicates (see [START_REF] Daval | Influence of amorphous silica layer formation on the dissolution rate of olivine at 90 degrees C and elevated pCO(2)[END_REF]Saldi et al.;[START_REF] Saldi | The role of Fe and redox conditions in olivine carbonation rates: An experimental study of the rate limiting reactions at 90 and 150 °C in open and closed systems[END_REF][START_REF] Sissmann | The deleterious effect of secondary phases on olivine carbonation yield: Insight from time-resolved aqueous-fluid sampling and FIB-TEM characterization[END_REF]) is not compatible with the intrinsic diffusion constant required to extend serpentinization duration up to several Gyrs. Therefore, alternative scenarios must be invoked to explain the observation of molecular hydrogen in Enceladus plume at present time.

A first avenue to develop in the future could be to relax some of the assumptions made in this study regarding the chemical composition of the solution, which controls the degree of undersaturation with respect to enstatite, olivine and SiO 2 (am), and represents a fundamental rate-controlling parameter of olivine and enstatite dissolution. Our calculations were made according to the assumption that the chemical composition of the fluid was buffered by stable mineral assemblages typical of serpentinization reactions, which is consistent with the theoretical considerations of [C Glein et al., 2018;[START_REF] Glein | The pH of Enceladus' ocean[END_REF], and confirmed experimentally at T ≥ 120 °C by e.g., [START_REF] Hsu | Ongoing hydrothermal activities within Enceladus[END_REF] or [START_REF] Sekine | High-temperature water-rock interactions and hydrothermal environments in the chondrite-like core of Enceladus[END_REF]. In this framework, the corresponding fluids have low [Mg 2+ ]/[SiO 2 (aq)] ratios, and the dissolution of olivine and enstatite proceeds at their maximum far-from-equilibrium dissolution rate.

However, the use of stable mineral assemblages to simulate the water composition on Enceladus at low temperature has recently been questioned: for instance [START_REF] Fukushi | Situ Formation of Monohydrocalcite in Alkaline Saline Lakes of the Valley of Gobi Lakes: Prediction for Mg, Ca, and Total Dissolved Carbonate Concentrations in Enceladus' Ocean and Alkaline-Carbonate Ocean Worlds[END_REF] suggested that the solubility of Mg may instead be controlled by metastable phases such as amorphous Mg-carbonates, based on the study of alkaline lakes from the Valley of Gobi (Mongolia), considered as relevant terrestrial analogues of ocean worlds. Similarly, [Sekine et al., 2015] suggested that the ocean water located at the boundary with the seafloor could be at saturation with respect to SiO 2 (am). Considering fluid compositions that are enriched in Mg and SiO 2 (aq) compared to the fluid compositions that we have considered here clearly has the potential to lengthen the lifetime of olivine and enstatite. A key question however would be to determine whether both elevated concentrations of Mg and SiO 2 (aq) can be maintained in strongly alkaline solutions (pH 9 to 11) at low temperatures, without precipitation of Mg-phyllosilicate that would ultimately bring back the solution to a composition typical of the mineral assemblages buffers considered in the present study. Strengthening this idea, [START_REF] Zeyen | Integrative analysis of the mineralogical and chemical composition of modern microbialites from ten Mexican lakes: What do we learn[END_REF] showed that when alkaline lakes also contain elevated concentrations of SiO 2 (aq), the Mg concentration is no longer controlled by amorphous Mg carbonate solubility, and leads to the precipitation of several Mg-phyllosilicates. As a consequence, the solid phases that control the chemical composition of fluids at low temperature remain a source of debate. Additional work, outside the goal of the present study, is needed to predict the Mg and SiO 2 (aq) concentrations in Enceladus' ocean in a consistent chemical model.

Another possibility that would reconcile the short timescales of H 2 production with observations could be that Enceladus is a very young moon ( [Ćuk et al., 2016]) only ~ 100 million years old. However, this explanation was challenged by a subsequent study ( [START_REF] Fuller | Resonance locking as the source of rapid tidal migration in the Jupiter and Saturn moon systems[END_REF]) which took into account the frequency dependence of the dissipation factor leading to an age consistent with 4.5 Gyr.

Alternatively, an explanation could be related to the fact that gravity vanishes at the center of Enceladus (it roughly varies linearly with radius), advection velocities are expected to decrease accordingly [START_REF] Le Reun | Internally Heated Porous Convection: An Idealized Model for Enceladus' Hydrothermal Activity[END_REF]). Even if permeability in this central region is identical to the one in the rest of the core, once the present water is consumed by the reaction, extremely reduced water supply is expected. This would lead to the preservation of an inner core experiencing basically no fluid circulation. This could result in the preservation of almost undisturbed olivine or organic matter at the center of Enceladus despite the hot temperatures, whose limited amounts of reaction could explain the present day production of H 2 as detected by Cassini.

Another hypothesis could be that the present day source of H 2 might not be due to serpentinization, which would have only occurred early on to explain the presence of nanoparticles of SiO 2 in the E-ring (note that those have not been found in the plume). The Zhu, C. (2005), In situ feldspar dissolution rates in an aquifer, Geochim Cosmochim Ac, 69(6), 1435-1453, doi:DOI 10.1016/j.gca.2004.09.005. Zolotov, M. Y. (2009), On the composition and differentiation of Ceres, Icarus, 204 (1), 183-193. Figure 5 1.E+03

Figure 8

Figure9.

Figure 9

1.E-10

1.E-08

1.E-06

1.E-04

1.E-02

1.E+00

1.E+02

1.E+04

Lifetime [Myr]

  Mg 2-x Fe x SiO 4 + (2-2x/3)H 2 O = (x/3)Fe 3 O 4 + (x/3)H 2 + SiO 2 (aq) + (2-x)Mg 2+ + (4-2x)OH -

Mg 1 -

 1 x Fe x SiO 3 + (4x/3-1)H 2 O + (2-2x)H + = (x/3)Fe 3 O 4 + SiO 2 (aq) + (1-x)Mg 2+ + (x/3) H 2 (4)

  The mineral assemblages that control the composition of the fluids circulating in the Enceladus core are a matter of debate, and will remain an open question as long as sampling the corresponding fluids will remain out of reach. In the following subsections, we estimate the reactivity of olivine and enstatite for fluid compositions considering two possible fluid compositions suggested in[C Glein et al., 2018].

Fig. 2 .

 2 Fig. 2. (a) Gibbs free energy of reaction associated to forsterite and enstatite dissolution supposing that the fluid composition is dictated by an equilibrium with respect to talc + chrysotile. As can be seen, these values never exceed -10 kJ/mol, and remain below the expected value of the Gibbs free energy of etch pit nucleation on forsterite (green dashed line) and enstatite (black dashed line). See Section 4 for details regarding the energetics of etch pit nucleation. (b) Impact of the fluid composition on olivine and enstatite dissolution rates in the TST framework.As can be seen, the dissolution rates remain within 90% of their far-from-equilibrium values.

Fig. 3 .

 3 Fig. 3. Main modeling results resulting from the assumption that the composition of the fluid circulating in the core of Enceladus is dictated by an equilibrium with Talc + Calcite + Tremolite + Chrysotile. The plots display the evolution of (a) the activity of CO 2 (aq), (b) [HCO 3 -], (c) pH and (d) the estimated lifetimes of forsterite and enstatite as a function of the temperature (which can be used as a proxy for the depth within the core). The lifetime of silicates is less than 15 Myr for the conditions relevant to Enceladus.

  framework of the transition state theory would modify the expected timespan of serpentinization in the context of Enceladus.4. Reactivity of olivine and enstatite controlled by etch pit nucleationIt has long been argued that laboratory-derived dissolution rate laws overestimate the actual rates of mineral dissolution in the field by up to several orders of magnitude (see [K.

  [J/m²] is the surface free energy, 𝑉 ̅ [mol/m 3 ] is the molar volume, µ [J/m 3 ] is the bulk shear modulus and b [m] is the length of the Burgers vector of the dislocation. Whereas the values of 𝑉 ̅ and µ are well known for forsterite and enstatite, the values of σ and b that should be considered are less well constrained. As discussed at length in [Pollet-Villard et al., 2016],

Fig. 4 .

 4 Fig. 4. (a) Gibbs free energy of reaction associated to SiO 2 (am) and cristobalite-α dissolution (considered as proxies for the solubility of ASSLs formed on olivine and enstatite, respectively), supposing that the fluid composition is dictated by an equilibrium with respect to talc + chrysotile or talc + calcite + tremolite + chrysotile. (b) Impact of the fluid composition on olivine and enstatite dissolution rates for a dissolution process controlled by the formation of ASSLs. As can be seen, the dissolution rates remain within 95% of their far-fromequilibrium values.

Fig. 5 .

 5 Fig. 5. Lifetime of enstatite and olivine as a function of temperature for a process controlled by diffusion in the nanoporosity of secondary coatings covering the reacted grains.

Fig. 7 .

 7 Fig. 7. Estimated flux of molecular hydrogen resulting from the oxidation of Fe 2+ released through enstatite dissolution for grain/aggregate radii ranging between 0.1 mm and 5 m. The zone located between the green dashed lines indicates today's flux of hydrogen as inferred from the measurements realized by the Cassini

Fig. 8 .

 8 Fig. 8. Estimated flux of molecular hydrogen resulting from the oxidation of Fe 2+ released through enstatite dissolution. The simulations were conducted with an arbitrary number of 7 concentric shells (see text). The zone located between the dashed red lines delineates today's flux of hydrogen as inferred from the measurements realized by the Cassini spacecraft. The vertical dashed line indicates the suspected age of Enceladus. The solid red line represents the outputs of the simulation conducted using D 0 as the intrinsic diffusion coefficient.

Fig. 9 .

 9 Fig.9. Summary of the estimated lifetime of Fe-bearing silicates (olivine or pyroxene) following the mechanisms investigated in Sections 2 to 6.2. As can be seen, if the rate-limiting factor of serpentinization is the hydrolysis rate of olivine/pyroxene (first 5 bars in the diagram above), the corresponding duration of serpentinization is extremely short, on the order of 10 Myr at most. If the supply of reactants (water, protons) is supposed to be rate-limiting, this extends the duration of serpentinization up to a few hundreds of Myrs. Such durations remain much lower than the most generally accepted age for Enceladus (4.5 Gyr).

  Figure1.

  Figure2

  Figure 3

  Figure 4

Figure7.

  Figure7.

  Figure 7

  When the numerical value of a parameter results from a calculation involving a range of parameter values, no specific value was assigned to it (-). :https://doi.org/10.1016/j.gca.2021.04.030.

	about	their	formation?,	Geochim	20.5(100); 16.9(150); 14.2(200) Cosmochim Ac, 305,	148-184,
	equilibrium constant of enstatite dissolution	𝑇 𝐾 𝐸𝑛𝑠	logK(T [°C]) = 12.5(0); 11.3(25); 9.7(60); 8.2(100); 6.7(150); 5.5(200)
	equilibrium constant of calcite dissolution (calculated from Eq. (22))	𝑇 𝐾 𝐶𝑎𝑙	logK(T [°C]) = 8.8(0); 8.2(25); 7.6(60); 6.8 (100); 6.6(150); 6.4(200)
	equilibrium constant of tremolite dissolution	𝑇 𝐾 𝑇𝑟	logK(T [°C]) = -66.9(0); -61.2(25); -53.2(60); -45.3(100); -37.3(150); -31.0(200)
	carbonate alkalinity inherited from the ocean	C ocean	0.03 molal
	activity coefficient of aqueous species i critical free energy of etch pit nucleation	𝛾 𝑖 ∆G crit	-𝑐𝑟𝑖𝑡 = -12.0 kJ/mol ; ∆𝐺 𝑒𝑛𝑠 ∆𝐺 𝑜𝑙 𝑐𝑟𝑖𝑡 = -7.5 kJ/mol
	diffusion coefficient of aqueous species at T intrinsic diffusion constant in nanopores	𝐷 𝑇 D 0	-7.2 × 10 -16 cm²/s
	activation energy of diffusion in nanopores	𝐸 𝑎,𝑑𝑖𝑓𝑓	15 kJ/mol
		porosity of Enceladus' core		𝛷	26% (unless otherwise mentioned)
		water to olivine molar ratio in serpentinization reactions		Y	1.5-1.25x/3 or 1.5-2x/3
		radius of Enceladus' core		𝑟 𝑐𝑜𝑟𝑒	190 km
		fluid velocity in Enceladus' core		𝜐	1 cm/yr
	radius of mineral aggregates in the core	𝑟 𝑎𝑔𝑔	100 µm ≤ r agg ≤ 5 m
		surface of spherical aggregates		𝑆 𝑎𝑔𝑔	4𝜋. 𝑟 𝑎𝑔𝑔	2
		volume of Enceladus' core		𝑉 0	(4𝜋 3 ⁄ ). 𝑟 𝑐𝑜𝑟𝑒 3
	hydrogen flux released from Enceladus	F H2	-
	hydrogen flux released from a single grain	F H2,grain	-

doi

(𝑆𝑖𝑂 2 (𝑎𝑞)) 𝑒𝑞 4 )1/3(29) where (SiO 2 (aq)) eq and (H + ) eq are calculated using Eq. (8b) and (28), respectively. As the values of (SiO 2 (aq)) eq , (Mg 2+ )eq and (H + ) eq are all fixed, ∆G Fo and ∆G Ens can be calculated.The lifetimes of olivine and enstatite were then calculated by combining Eq. (5) and (16), and are represented in Fig.3d. Again, the mineral assemblage dictating the fluid composition in this configuration results in far-from-equilibrium conditions with respect to forsterite and enstatite, if one considers that the dissolution rate laws of those mineral observe the transition state theory, thus resulting in short lifetimes, as in the previous Section.

 (43)

Acknowledgements

The authors thank colleagues from ISTerre (Grenoble, France) and in particular Amber Zandanel, Roland Hellmann, Laurent Truche, Benjamin Malvoisin and Fabrice Brunet for stimulating discussions about serpentinization reactions in the Enceladus' core. CS acknowledges the support of NASA's 'Exploring Ocean Worlds' Research Coordinated Network, managed by the Woods Hole Oceanographic Institution. Finally, we thank Prof.

Sekine (AE) and an anonymous reviewer for their stimulating comments that improved the manuscript.

Data availability statement

All shapefiles produced in this work are available at [START_REF] Daval | _JGRPlanets_H2GenerationTimescalesEnceladus[END_REF] 20211206_JGRPlanets_H2GenerationTimescalesEnceladus

[Dataset]. Zenodo.

https://doi.org/10.5281/zenodo.5761723. All data used in the present study come from previously published sources and in particular from [START_REF] Rimstidt | Systematic review of forsterite dissolution rate data[END_REF] for olivine and [START_REF] Palandri | A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modelingRep[END_REF] for enstatite dissolution kinetics, [van der Lee and de [START_REF] Van Der Lee | CHESS Tutorial and Cookbook[END_REF] for mineral solubility constants, [Antonio C. [START_REF] Lasaga | Surface chemistry, etch pits and mineral-water reactions[END_REF] for etch pit nucleation energetics, [START_REF] Daval | Influence of amorphous silica layer formation on the dissolution rate of olivine at 90 degrees C and elevated pCO(2)[END_REF] and [START_REF] Daval | Linking nm-scale measurements of the anisotropy of silicate surface reactivity to macroscopic dissolution rate laws: New insights based on diopside[END_REF] for the solubility of ASSLs developed on olivine and pyroxene, respectively, [START_REF] Rotenberg | Water and ions in clays: Unraveling the interlayer/micropore exchange using molecular dynamics[END_REF] for the activation energy of diffusion in nanopores, [START_REF] Daval | Influence of amorphous silica layer formation on the dissolution rate of olivine at 90 degrees C and elevated pCO(2)[END_REF] for estimated diffusivity of ASSLs on olivine and [START_REF] Waite | Cassini finds molecular hydrogen in the Enceladus plume: evidence for hydrothermal processes[END_REF] for the flux of hydrogen in Enceladus' plume.

presence of heavy organics ( [START_REF] Postberg | Macromolecular organic compounds from the depths of Enceladus[END_REF]) raises the possibility that their pyrolysis could produce the present time H 2 [START_REF] Sotin | Carbon-rich composition models of icy moons and dwarf planets[END_REF]). It would also explain at the same time similar amounts of CH 4 and H 2 in the plumes, which are difficult to explain by serpentinization without invoking contribution of alive microorganisms [START_REF] Affholder | Bayesian analysis of Enceladus's plume data to assess methanogenesis[END_REF]). In that case, a low porosity inner-core of Enceladus made of an assemblage of hydrated silicates and organic matter, and an outer core made of porous hydrated silicates required for tidal heating could explain the Cassini data. Laboratory experiments measuring the kinetics of such process at temperatures relevant to Enceladus are necessary to validate / rule out such a scenario.

Finally, the oxidation kinetics of Fe 2+ from hydrous silicates (e.g. from Fe 2+ -serpentine minerals formed in an early serpentinization stage) and the transport properties of secondary phases formed under those conditions are poorly known, and might provide an additional avenue for explaining long H 2 production times and low fluxes.

Appendix A: Summary table of model parameters