Benjamin Malvoisin
email: benjamin.malvoisin@univ-grenoble-alpes.fr

Yury Y Podladchikov

Artem V Myasnikov

Achieving complete reaction while the solid volume increases: a numerical model applied to serpentinisation

Keywords: serpentinisation, numerical model, reaction rate, tectonic deformation, solid volume increase, fluid pathways, porosity clogging, crystallisation pressure

Solid volume increases during the hydration of mantle rocks. The fluid pathways necessary to feed the reaction front with water can be filled by the low density reaction products. As a result, the reaction front dries out and the reaction stops at low reaction progress. This process of porosity clogging is generally predicted to dominate in reactive transport models, even when processes such as reaction-induced fracturing are considered. These predictions are not consistent with observations at mid-ocean ridges where dense mantle rocks can be completely replaced by low density serpentine minerals. To solve this issue, we develop a numerical model coupling reaction, fluid flow and deformation. High extents of reaction can only be achieved when considering that the increase in solid volume during reaction is accommodated through deformation rather than porosity clogging. The model can generate an overpressure that depends on the extent of reaction and on the boundary conditions. This overpressure induces viscoelastic compaction that limits the extent of the reaction. The serpentinisation rate is therefore controlled by the accommodation of volume change during reaction, and thus by deformation, either induced by the reaction itself or by tectonic processes.

Introduction

Olivine, the main component of mantle peridotites, hydrates at temperature below ∼ 400 • C to form serpentine ± brucite ± magnetite ± hydrogen [START_REF] Evans | Metamorphism of Alpine Peridotite and Serpentinite[END_REF]. This so-called reaction of serpentinisation requires fluid transport along pathways often formed by tectonic deformation. In marine geodynamic settings, the fluid pathways include detachment faults at slow-spreading mid-ocean ridges [START_REF] Cannat | Emplacement of mantle rocks in the seafloor at mid-ocean ridges[END_REF][START_REF] Co. Escartín | Central role of detachment faults in accretion of slow-spreading oceanic lithosphere[END_REF][START_REF] Rüpke | Global rates of mantle serpentinization and H2 production at oceanic transform faults in 3-D geodynamic models[END_REF], bend-faults at ocean trenches [START_REF] Ranero | Bending-related faulting and mantle serpentinization at the Middle America trench[END_REF][START_REF] Van Avendonk | Structure and serpentinization of the subducting Cocos plate offshore Nicaragua and Costa Rica[END_REF], and extensional faults at passive margins [START_REF] Bayrakci | Fault-controlled hydration of the upper mantle during continental rifting[END_REF]. Here we focus on serpentinisation during mantle rock exhumation in the footwall of detachment faults at slow spreading mid-ocean ridges.

Serpentinisation induces changes in the behaviour of the oceanic lithosphere due to strong differences in the physico-chemical properties of the reactants and the products. Serpentine content above ∼ 10 % reduces the rock strength by decreasing the friction coefficient from 0.6 to ∼ 0.4 [START_REF] Escartín | Effects of serpentinization on the lithospheric strength and the style of normal faulting at slow-spreading ridges[END_REF] with consequences on the rheology of the detachment faults active at mid-ocean ridges [START_REF] Bickert | How do detachment faults form at ultraslow mid-ocean ridges in a thick axial lithosphere?[END_REF]. Magnetite is a ferrimagnetic mineral with a magnetisation four orders of magnitudes higher than olivine (Malvoisin et al., 2012b). Its precipitation during serpentinisation can generate magnetic anomalies on the seafloor [START_REF] Szitkar | The magnetic signature of ultramafic-hosted hydrothermal sites[END_REF]. Hydrogen is an energy source for chemosynthetic organisms and participates in the abiotic synthesis of organic compounds possibly involved in the origin of life [START_REF] Russell | Serpentinization as a source of energy at the origin of life[END_REF][START_REF] Mccollom | Serpentinites, Hydrogen, and Life. Elements 9[END_REF]. Serpentinisation is also an exothermic reaction possibly modifying the heat budget at mid-ocean ridges [START_REF] Macdonald | Rate of serpentinization in seafloor environments[END_REF]. Predicting the evolution of the oceanic lithosphere thus requires to model the rate and the extent of serpentinisation at mid-ocean rdiges.

Serpentinisation is already modelled at the kilometre scale to investigate its influence on hydrothermal circulation, heat budget and geodynamics [START_REF] Emmanuel | Suppression and stimulation of seafloor hydrothermal convection by exothermic mineral hydration[END_REF][START_REF] Delescluse | Serpentinization pulse in the actively deforming Central Indian Basin[END_REF][START_REF] Iyer | Feedbacks between mantle hydration and hydrothermal convection at ocean spreading centers[END_REF][START_REF] Bickert | How do detachment faults form at ultraslow mid-ocean ridges in a thick axial lithosphere?[END_REF]. Reaction rate is assumed to be controlled by a single process in these models. [START_REF] Emmanuel | Suppression and stimulation of seafloor hydrothermal convection by exothermic mineral hydration[END_REF] and [START_REF] Iyer | Feedbacks between mantle hydration and hydrothermal convection at ocean spreading centers[END_REF] used a law derived from experiments of olivine powder hydration [START_REF] Martin | Some experimental and theoretical observations on the kinetics of hydration reactions with particular reference to serpentinization[END_REF] whereas [START_REF] Bickert | How do detachment faults form at ultraslow mid-ocean ridges in a thick axial lithosphere?[END_REF] relate serpentinisation rate to the shear deformation. It remains unclear which assumption should be preferred to model serpentinisation at large scale. Moreover, the used models are either incompressible, considering a constant porosity or using a parametrized law for the evolution of porosity. Such approaches are justified by the need to use simplifying assumptions when modelling processes at different time and space scales. However, they cannot take into account the impact of hydration on density and porosity whereas these parameters strongly modify deformation and fluid flow.

Volume change during reaction is the key parameter for modelling the couplings between reaction, deformation and fluid flow. Serpentinisation induces a change in solid volume (∆V /V ) of ∼ 50 % in a closed system [START_REF] Macdonald | Rate of serpentinization in seafloor environments[END_REF][START_REF] O'hanley | Solution to the volume problem in serpentinization[END_REF]. Mass transfer through aqueous species transport can reduce the increase in solid volume [START_REF] Thayer | Serpentinization considered as a constant-volume metasomatic process[END_REF][START_REF] Carmichael | Induced Stress and Secondary Mass Transfer: Thermodynamic Basis for the Tendency toward Constant-Volume Constraint in Diffusion Metasomatism[END_REF][START_REF] Fletcher | Mineral growth in rocks: Kineticrheological models of replacement, vein formation, and syntectonic crystallization[END_REF]. However, Malvoisin et al. (2020b) and [START_REF] Klein | Quantifying the volume increase and chemical exchange during serpentinization[END_REF] measure ∆V /V in natural samples and in experiments close to the theoretical value for reaction in a closed system, suggesting that mass transfer is limited, at least during the incipient stage of the reaction. The serpentinisation rate directly depends on the evolution of the fluid pathways during reaction.

They initially form by thermal cracking [START_REF] Demartin | Experimental constraints on thermal cracking of peridotite at oceanic spreading centers[END_REF][START_REF] Boudier | Serpentine Mineral Replacements of Natural Olivine and their Seismic Implications: Oceanic Lizardite versus Subduction-Related Antigorite[END_REF] or tectonic deformation [START_REF] Rouméjon | Serpentinization and Fluid Pathways in Tectonically Exhumed Peridotites from the Southwest Indian Ridge (62-65 E)[END_REF], and then evolve during reaction as a consequence of ∆V /V . The precipitation of reaction products in the pores reduces the permeability by several orders of magnitude [START_REF] Godard | Incipient hydration of mantle lithosphere at ridges: A reactive-percolation experiment[END_REF][START_REF] Farough | Evolution of fracture permeability of ultramafic rocks undergoing serpentinization at hydrothermal conditions : An experimental study[END_REF]. However, a porosity is observed at the nanoscale in serpentinized peridotites [START_REF] Tutolo | Nanoscale constraints on porosity generation and fluid flow during serpentinization[END_REF]. Fluid pathways could thus be preserved during reaction, even though this should be confirmed with the acquisition of data on pore connectivity at the nanoscale. The positive change in volume during reaction also generates stresses (i.e. crystallisation pressure) sufficient to fracture the rock [START_REF] Jamtveit | Reaction induced fracturing during replacement processes[END_REF][START_REF] Kelemen | Reaction-driven cracking during retrograde metamorphism: Olivine hydration and carbonation[END_REF][START_REF] Plümper | The interface-scale mechanism of reaction-induced fracturing during serpentinization[END_REF][START_REF] Malvoisin | Control of serpentinisation rate by reaction-induced cracking[END_REF].

The positive ∆V /V is considered in some models of serpentinisation. It induces stress build-up and crack generation in the models of reaction-induced fracturing [START_REF] Rudge | A simple model of reaction-induced cracking applied to serpentinization and carbonation of peridotite[END_REF][START_REF] Ulven | Reaction-driven fracturing of porous rock[END_REF][START_REF] Shimizu | The roles of fluid transport and surface reaction in reaction-induced fracturing, with implications for the development of mesh textures in serpentinites[END_REF][START_REF] Malvoisin | Control of serpentinisation rate by reaction-induced cracking[END_REF][START_REF] Zhang | Modeling Porosity Evolution Throughout ReactionInduced Fracturing in[END_REF]. These models predict an increase in reactive surface area during reaction and reproduce the mesh microtexture observed in natural samples. However, they either do not take into account Darcian flow, or they do not conserve fluid mass, and are consequently not able to quantitatively model porosity evolution during reaction [START_REF] Evans | Phase-Field Modeling of Reaction-Driven Cracking: Determining Conditions for Extensive Olivine Serpentinization[END_REF].

Reactive transport models such as TOUGHREACT [START_REF] Xu | Toughreact version 2.0: A simulator for subsurface reactive transport under non-isothermal multiphase flow conditions[END_REF] incorporate the volumetric effect of reaction on transport properties. They predict a strong reduction in porosity during reaction inducing a positive ∆V /V , preventing the achievement of the high extents of reaction measured in natural samples [START_REF] Liu | Reactive transport modeling of mineral carbonation in unaltered and altered basalts during CO2 sequestration[END_REF][START_REF] Marini | Geological Sequestration of Carbon Dioxide: Thermodynamics, Kinetics, and Reaction Path Modeling[END_REF]. This limitation probably arises from the omission in these models of mechanisms responsible for fluid pathways generation or preservation. [START_REF] Evans | Phase-Field Modeling of Reaction-Driven Cracking: Determining Conditions for Extensive Olivine Serpentinization[END_REF] consider reaction-induced fracturing in a model combining reaction and fluid flow. They were able to model the effect on fluid pathway evolution of the competition between cracking and porosity clogging. However, high extents of reaction (> 60 %) are only achieved in the cracks with limited hydration in the rock groundmass. The main issue regarding serpentinisation modelling thus consists in being able to predict the high extents of reaction observed in the natural samples.

The present contribution aims at proposing a new model of serpentinisation addressing this issue. We derive a new set of equations to model the couplings between reaction, deformation and fluid flow. Based on observations in natural samples, we consider that the solid volume increase is accommodated through deformation rather than porosity clogging. The new model is applied to serpentinisation at mid-ocean ridges in the footwall of a detachment fault. We model fluid-rock interaction during hydrothermal circulation with a one-dimensional numerical model. The impact of the parameter values on the model predictions is explored. On the contrary to previous modelling attempts, the new approach allows high extents of reaction to be reached and reaction-induced pressure increase to be predicted. The model is finally used to discuss the processes controlling the serpentinisation rate.

Modelling the couplings between reaction, deformation and fluid flow

The model considers a porous matrix of dunite reacting with water in the stability field of serpentine + brucite. Serpentinisation is assumed to be isochemical except for the addition of water based on the study of Malvoisin et al. (2020b). Si, Mg and Fe are thus not transported as aqueous species between the model grid points. We use here a two-phase continuum medium approach by defining macroscopic properties such as porosity and extent of reaction for a mineralogical assemblage potentially composed at each grid point of both olivine and its reaction products. This approach is valid at a scale larger than the olivine grain size so that the rock properties can be averaged over a significant number of grains. Olivine grains reacting during serpentinisation have a typical size of 100 µm [START_REF] Malvoisin | Control of serpentinisation rate by reaction-induced cracking[END_REF]. The model can thus be used from a scale of ∼ 1 mm. The model parameters and their values are given in Table 1. The model is based on the model developed in [START_REF] Malvoisin | Coupling changes in densities and porosity to fluid pressure variations in reactive porous fluid flow: Local thermodynamic equilibrium[END_REF] and in [START_REF] Co. Omlin | Pore Fluid Extraction by Reactive Solitary Waves in 3-D[END_REF]. The thermodynamic admissibility of the systems of equations derived below is also verified following the procedure described in Yarushina and Podladchikov (2015) (see supplementary materials for details).

The mass conservation equations can be written for the total system as:

∂ (ρ w φ + ρ s (1 -φ)) ∂t + ∇ j ρ w φv j w + ρ s (1 -φ) v j s = 0 (1)
, and for the dry solid as:

∂ (ρ s (1 -X s ) (1 -φ)) ∂t + ∇ j ρ s (1 -X s ) (1 -φ) v j s = 0 (2)
where X s , φ, ρ, v j and ∇ j are the mass fraction of water bound to the solid, the porosity, the density, the j components of the velocity vector and the j components of the del operator, respectively. The solid properties are denoted with a subscript s, and the water properties are denoted with a subscript w. This particular choice of two independent mass conservation equations is motivated by their independence on (de)hydration reactions and the absence of reactionrelated source terms.

The conservation of fluid momentum is modelled with the Darcy's law:

φ v j w -v j s = - k φ η w ∇ j P w + ρ w ge j z ( 3
)
where k φ is the permeability, η w is the fluid viscosity, P w is the fluid pressure and e j z are the components of an upward directed unit vector. Combining equation 3 with equations 1 and 2 allows water velocity to be eliminated.

The [START_REF] Connolly | Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation[END_REF] solid density evolution:

1 ρ s dρ s dt = β s dP s dt - P -P w (1 -φ) 2 dφ dt -β sr 1 1 -X s dX s dt (4)
where β s is the solid compressibility measured in an unjacketed test, β sr is the density change associated with hydration (X s modification), and P s is the solid pressure. d dt is the material derivative of the solid phase defined as

d dt = ∂ ∂t + v j s ∇ j .
P s is combined with P w to calculate the total pressure (P ) as:

P = φP w + (1 -φ) P s (5)
The change in fluid density is calculated by introducing the fluid compressibility (β w ):

1 ρ w d w ρ w dt = β w d w P w dt (6)
where dw dt is the material derivative of the water defined as dw dt = ∂ ∂t + v j w ∇ j Finally, the change in composition (X s ) is modelled with a first-order kinetic law:

dX s dt = φ φ p 2/3 X seq -X s τ ( 7
)
where τ is the characteristic time for reaction and X seq is the water mass fraction bound to the solid at the equilibrium, that is in the serpentine + brucite assemblage. The evolution of reactive surface area during reaction is modelled as a function of porosity with the geometrical relationship proposed by [START_REF] Lichtner | The quasi-stationary state approximation to coupled mass transport and fluid-rock interaction in a porous medium[END_REF] and [START_REF] Kieffer | An experimental study of the reactive surface area of the Fontainebleau sandstone as a function of porosity, permeability, and fluid flow rate[END_REF]. φ p is the minimum porosity for which a porous rock has the same reactive surface area than a powder. We use a value of 0.38 in the following based on the experimental results of [START_REF] Llana-Fúnez | Experimental dehydration kinetics of serpentinite using pore volumometry[END_REF]. Equations 1, 2, 3, 4, 6, 5 and 7 form a system of nine equations in thirteen unknowns, 7 scalars (ρ w , ρ s , φ, X s , P w , P s and P ) and 6 vector components (v j w and v j s ). Solving this system of equations thus requires four additional equations.

The conservation of total momentum could first be used [START_REF] Malvoisin | Coupling changes in densities and porosity to fluid pressure variations in reactive porous fluid flow: Local thermodynamic equilibrium[END_REF] but this would require to solve for shear stress evolution which is beyond the scope of this study. We thus explore the behaviour of the system of equations for two cases: (i) a constant solid volume (∇ j v j s = 0); (ii) a constant total pressure (P = constant). These approximations correspond to idealized systems in which reaction-induced change in volume is either not accommodated (i) or fully accommodated (ii) by deformation. The behaviour of a real system is expected to be intermediate between these two extreme cases. Nevertheless, these approximations provide maximum estimates for total pressure variation during reaction (i) and reaction-induced deformation (ii).

The last equation necessary to close the system of equations is the equation for solid deformation. We show below why this equation plays a key role for modelling solid volume change during reaction. We first present the classical form of this equation usually used in reactive transport models. We then show how it should be modified to take into account observations at the nanoscale of limited mass transfer during serpentinisation.

Classical approach

The classical approach refers here to models in which reaction-induced deformation is not considered in the equations for solid deformation. This approach is used in most of the models considering the couplings between reaction and deformation [START_REF] Balashov | Modeling metamorphic fluid flow with reaction-compaction-permeability feedbacks[END_REF][START_REF] Xu | Toughreact version 2.0: A simulator for subsurface reactive transport under non-isothermal multiphase flow conditions[END_REF][START_REF] Tian | The impact of porosity waves on crustal reaction progress and CO2 mass transfer[END_REF][START_REF] Malvoisin | Coupling changes in densities and porosity to fluid pressure variations in reactive porous fluid flow: Local thermodynamic equilibrium[END_REF][START_REF] Tian | The Potential for Metamorphic Thermal Pulses to Develop During Compaction-Driven Fluid Flow[END_REF]Malvoisin et al., 2020a;[START_REF] Evans | Phase-Field Modeling of Reaction-Driven Cracking: Determining Conditions for Extensive Olivine Serpentinization[END_REF][START_REF] Schmalholz | 2D HydroMechanicalChemical modelling of (de)hydration reactions in deforming heterogeneous rock: The periclasebrucite model reaction[END_REF]. With the classical approach, solid deformation is modelled independently from reaction for example with a visco-elastic rheology (Yarushina and Podladchikov, 2015):

∇ j v j s = -β d dP dt + (β d -β s ) dP w dt - P -P w (1 -φ) η φ (8)
where β d and η φ are respectively, the drained compressibility and the pore viscosity.

Re-arranging equations 1, 2, 3, 4, 5, 6 and 7 with additional equation 8 provides a closed system of three equations describing the system evolution with three unknown scalar functions, porosity (φ), fluid pressure (P w ) and mass fraction of water in the solid (X s ):

                   dX s dt = φ φ p 2/3 X seq -X s τ ∇ j k φ η w ∇ j P w + ρ w ge j z = -β et dP dt + β ew dP w dt - P -P w (1 -φ) η φ + ρ s ρ w -(1 + β sr ) 1 -φ 1 -X s dX s dt 1 (1 -φ) dφ dt = β et - β s φ 1 -φ - dP dt + dP w dt - P -P w (1 -φ) η φ - 1 + β sr 1 -X s dX s dt (9)
where β et and β ew are the total and fluid effective compressibilities, respectively, defined as

β et = β d -β s , (10)
and

β ew = β et + φ (β w -β s ) . ( 11
)
Reactive terms appear in the fluid pressure and the porosity evolution equations in 9. This implies that volume change during reaction is accommodated through precipitation in the pores and fluid pressure variation. Mineralogical reactions in the presence of fluid proceed through four steps: the transport of fluid up to the primary minerals, the dissolution of the primary minerals, the transport of aqueous species and the precipitation of the secondary minerals 1); primary mineral (in yellow) dissolution (2); aqueous species transport (3); and secondary phase (in green) precipitation (4). A: the classical approach consists in considering that aqueous species are transported up to the pores where precipitation occurs. This process does not generate deformation directly. B: the new approach proposed here is based on observations at the nanoscale showing that aqueous species transport can be limited to a scale of less than 100 nm. As a result, the reaction products do not precipitate in the pores but at the surface of the primary mineral. This implies that reaction induces rock deformation (i.e. swelling; red arrows).

(Figure 1). Modelling fluid-mediated reactions thus requires to assume that aqueous species transport occurs. We do not consider this process in our equations, implying that the occurrence of aqueous species transport is restricted to a scale below the spacing of the grid points in our model (< 1 mm). However, we directly relate here porosity evolution to reaction. This implicitly requires that mass transfer efficiently occurs at least at the micrometer scale from the olivine surface up to the pores (Figure 1A). This latter assumption may not be verified during hydration reactions as suggested by observations in natural samples [START_REF] Beard | Onset and progression of serpentinization and magnetite formation in Olivine-rich troctolite from IODP hole U1309D[END_REF][START_REF] Boudier | Serpentine Mineral Replacements of Natural Olivine and their Seismic Implications: Oceanic Lizardite versus Subduction-Related Antigorite[END_REF]Malvoisin et al., 2020b). We propose in the following a way to take into account these observations by modifying equation 8.

A new approach: reaction-induced swelling

At the submicrometre scale, [START_REF] Beard | Onset and progression of serpentinization and magnetite formation in Olivine-rich troctolite from IODP hole U1309D[END_REF] and Malvoisin et al. (2020b) report reaction products dominated by a serpentine + brucite mixture with a composition similar to olivine except for the addition of water. [START_REF] Boudier | Serpentine Mineral Replacements of Natural Olivine and their Seismic Implications: Oceanic Lizardite versus Subduction-Related Antigorite[END_REF] show the presence of a ∼ 100 nm thick zone coating the olivine surface, and composed of pre-hydrated olivine. Serpentinisation is thus isochemical at least at the hundred of nanometre scale with limited aqueous species transport, probably by diffusion at grain boundaries [START_REF] Tutolo | Nanoscale constraints on porosity generation and fluid flow during serpentinization[END_REF] (Figure 1B). Such a replacement process with limited transport is not expected to significantly modify the existing pores. This implies that the positive change in solid volume occurring during reaction is not accommodated by mass transfer up to the pores but rather by deformation (Figure 1B). All these observations can be taken into account by modifying equation 8 such that porosity evolution does not depend on volume change during reaction:

∇ j v j s = -β d dP dt + (β d -β s ) dP w dt - P -P w (1 -φ) η φ + 1 + β sr 1 -X s dX s dt (12)
The last term in equation 12 is the deformation induced by solid volume change. With this new equation for solid momentum conservation, a new closed system of three equations for three unknown functions, φ, P w and X s is obtained by also using equations 1, 2, 3, 4, 5, 6 and 7:

                             dX s dt = φ φ p 2/3 X seq -X s τ ∇ j k φ η w ∇ j P w + ρ w ge j z = -β et dP dt + β ew dP w dt - P -P w (1 -φ) η φ + ρ s (1 -φ) ρ w + φ (1 + β sr ) 1 1 -X s dX s dt 1 (1 -φ) dφ dt = β et - β s φ 1 -φ - dP dt + dP w dt - P -P w (1 -φ) η φ (13)

Numerical methods, initial set-up and parameters

The numerical model aims at determining the conditions necessary to reach high reaction progress during serpentinisation at mid-ocean ridges. We build a simple idealized model for fluid-rock interaction during hydrothermal circulation in exhumed mantle rocks. We do not consider shear deformation and the strong changes in rock composition occurring during extensive alteration. The model thus rather applies to a region located outside the shear zones involved in exhumation in the footwall of a detachment fault . We consider a one-dimensional domain of length L undergoing serpentinisation at 200 • C and 50 MPa (initial fluid and total pressures). We impose a constant fluid pressure of 50 MPa at the model boundaries during simulation. The domain is assumed horizontal so that gravity terms can be neglected. It is initially exclusively composed of olivine (dunite; X s = 0). The initial porosity is of φ 0 = 0.01 to allow for fluid transport and reaction initiation. We predict the evolution of porosity, pressures and rock hydration degree with the systems of equations 9 and 13 for a wide range of model parameter values.

The numerical scheme consists in solving nondimensional versions of equations 9 and 13 with an explicit algorithm using pseudotransient iterations [START_REF] Kelley | Convergence Analysis of Pseudo-Transient Continuation[END_REF]. This algorithm is validated for the system of equations 9 by reproducing an analytical solution in [START_REF] Co. Omlin | Pore Fluid Extraction by Reactive Solitary Waves in 3-D[END_REF]. We also produce an analytical benchmark for the system of equations 13. The algorithm is validated for this system of equations by reproducing an analytical prediction (see the derived analytical solution in the supplementary materials). The non-dimensional equations are obtained by choosing the model scale as the length scale (δ = L), the inverse of the compressibility as the pressure scale (p = 1 β ef 0 ) and the characteristic time for viscous compaction as the time scale (t = η φ p ). This leads to the introduction of two nondimensional numbers describing the system behaviour : the Dahmköhler number comparing the characteristic times for viscoelastic relaxation and for reaction (Da =

β ef 0 η φ τ
); the Deborah number comparing the characteristic times for viscoelastic relaxation and for viscous compaction (De = k φ η φ ηwL 2 ). The non-dimensional systems of equations solved in the numerical models can be found in the supplementary materials. For clarity, we present the following results with dimensional quantities.

The values of the parameters used in the model are given in Table 1. Some parameters such as compressibility are well-constrained. β w is calculated with Perple X [START_REF] Connolly | Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation[END_REF] by computing the derivative of ρ w with pressure at 200 • C and 50 MPa. We also use Perple X to compute the change in solid density as a function of X s that is used to calculate β sr . For the other less constrained parameters, we provide in the following ranges of possible values which are used to compute the Da and De values to be used in the simulations.

Observations in natural samples reveal that preferred fluid pathways separate 10 cm-to 100 m-wide domains undergoing serpentinisation [START_REF] Rouméjon | Serpentinization and Fluid Pathways in Tectonically Exhumed Peridotites from the Southwest Indian Ridge (62-65 E)[END_REF]. We use this size range as the length of the model (L). Permeability measured during fluid flow experiments on peridotite takes values comprised between 10 -20 and 10 -16 m 2 [START_REF] Godard | Incipient hydration of mantle lithosphere at ridges: A reactive-percolation experiment[END_REF][START_REF] Farough | Evolution of fracture permeability of ultramafic rocks undergoing serpentinization at hydrothermal conditions : An experimental study[END_REF]. Pore viscosity spans over four orders of magnitude (10 17 to 10 20 P a.s; [START_REF] Hilairet | High-Pressure Creep of Serpentine, Interseismic Deformation, and Initiation of Subduction[END_REF]). The characteristic time for dunite serpentinisation are calculated for an initial olivine grain size of 100 µm based on kinetic laws determined with experiments on powders. The main parameter controlling reaction rate is the fluid composition. The slower rates are measured in experiments using pure water (τ = 2.10 8 s) (Malvoisin et al., 2012a) whereas the faster rates are measured with alkaline fluids (τ = 10 7 s) [START_REF] Lafay | Mineral replacement rate of olivine by chrysotile and brucite under high alkaline conditions[END_REF]. We use in the following τ values comprised between 10 7 s and 10 9 s. These estimates of parameters are used to determine Da and De values ranging from 10 -1.15 to 10 3.8 and 10 -1.5 to 10 5 , respectively. For each system of equations, we perform simulations with the two assumptions of constant solid volume (i) and constant total pressure (ii). We run 930 simulations at various De and Da values for each pair of solved equations/assumption (4 x 930 simulations in total). The simulations are stopped when the minimum porosity is below 10 -15 or when the maximum extent of reaction is above 99 %.

Porosity clogging with the classical approach

We first use the system of equations 9 to model serpentinisation. Reaction leads to fluid pressure decrease due to water consumption by reaction. This Figure 2: Results of simulations using the classical approach (equations 7 and 9). The evolutions with time of extent of reaction (A), fluid pressure (plain line) and total pressure (dashed line) (B), and porosity (C) are displayed. The models displayed with black and red lines consider a constant volume (i) or a constant total pressure (ii) for closing the system of equations, respectively. The blue line in subpanel A corresponds to the maximum extent of reaction predicted for a non-deforming solid matrix (ξmax; see text for details). The models are run with k φ = 10 -19 m 2 , L = 10m, τ = 5.10 7 s and η φ = 10 18 P a.s (De = 3; Da = 14) decrease highly depends on the De and Da numbers used in the simulations when a constant solid volume (i) is assumed (Figure S2B). At high De and low Da, P w variations are limited to several megapascals (Figure 2) whereas, at low De and high Da, it can reach several tens of megapascals (Figure S2). The water vapor pressure at 200 • C (1.6 MPa) can even be achieved. The model does not take into account the effect of vaporisation on density evolution and is thus not relevant for fluid pressures below the water vapor pressure. The fluid pressure decrease is limited to less than 3 MPa when a constant total pressure (ii) is used (Figures 2 B and S3 B). The difference between fluid and total pressure indeed drives viscoelastic compaction, limiting fluid pressure variations relative to total pressure. For all the parameter values investigated here, the model predicts a similar decrease of the reaction rate due to porosity clogging by reaction products (Figures 2 A,S2 A and S3 A). The reaction ultimately stops when the porosity and the reactive surface area are zero. The maximum extent of reaction reached in the simulations is ξ max = 3%. The maximum reaction extent can be predicted by integrating the porosity evolution equation in the system of equation 9 when assuming no deformation and reaction thus limited by the free space initially available (φ 0 ). The obtained expression is ξ max = φ0 Xseq(1+βsr) giving a value of 2.94% for φ 0 = 0.01, consistent with the model predictions. This expression of ξ max also provides the minimum initial porosity required to completely serpentinise a rock with the model based on equations 9. We calculate here a minimum initial porosity of 34 %.

Complete serpentinisation with reaction-induced deformation

The results of numerical modelling with equations 13 significantly depend on the assumption used for closing the system of equations. When assuming a constant solid volume (i), the evolutions of extent of reaction and porosity are similar to the classical approach. The extent of reaction reaches a maximum close to ξ max (Figures 3 A andS4 A) as the porosity decreases to become zero (Figures 3 C). The mechanism of porosity decrease is however different since it is due to mechanical closure here whereas it is due to the precipitation of reaction products in the simulations performed with the classical approach. The driving force for viscoelastic compaction is the total pressure increase as a result of solid volume increase during reaction (Figure 3 C). We calculate a maximum total pressure increase of 14 MPa along a segment spanning from Da = 10 and De = 10 1.5 to Da = 10 4 and De = 10 3 in Figure S4 B. This indicates that reaction-induced pressure increase is favored by a high reaction rate, a high viscosity and/or a low effective compressibility.

We observe two types of behaviour when assuming a constant total pressure (ii). At high Da and low De, porosity clogging occurs as a a result of mechanical compaction as observed with assumption (i) (Figure S5). The final extent of reaction is close to ξ max . At low Da and high De, the porosity and the fluid pressure only change during simulation by less than 0.5 % and several megapascals, respectively (Figure 3 E andF). Fluid pathways are thus preserved, and the extent of reaction reaches 100 % (Figures 3 D andS5 A). The solid volume increase is accommodated by the expansion of the rock matrix. As porosity is rather constant during reaction, the duration of the reaction only depends on τ and φ 0 according to equation 7. The simulation displayed in Figure 3 reaches for example completion after approximately 150 yr.

Discussion and conclusion

Several numerical approaches have been used to model serpentinisation. Most of them investigate the influence of positive volume increase during reaction on fluid pathways formation [START_REF] Rudge | A simple model of reaction-induced cracking applied to serpentinization and carbonation of peridotite[END_REF][START_REF] Ulven | Reaction-driven fracturing of porous rock[END_REF][START_REF] Malvoisin | Control of serpentinisation rate by reaction-induced cracking[END_REF][START_REF] Zhang | Modeling Porosity Evolution Throughout ReactionInduced Fracturing in[END_REF]. Nevertheless, the most complete models coupling reaction, stress variation, fluid flow and permeability do not produce complete serpentinisation in the entire rock [START_REF] Evans | Phase-Field Modeling of Reaction-Driven Cracking: Determining Conditions for Extensive Olivine Serpentinization[END_REF]. Reactive transport models such as TOUGHREACT [START_REF] Xu | Toughreact version 2.0: A simulator for subsurface reactive transport under non-isothermal multiphase flow conditions[END_REF]) also fail to reach high extents of reaction while the solid volume increases during reaction [START_REF] Paukert | Mineral Carbonation in Mantle Peridotite of the Samail Ophiolite , Oman : Implications for permanent geological carbon dioxide capture and storage[END_REF]. All these models are based on the same principle as the classical approach detailed here in section 2.1 with a decoupling of reaction and deformation. Positive volume increase can only be accommodated by precipitation in the pores. Mass transfer is assumed to occur at least between the surface of the reacting primary minerals and the pores. The fluid pathways are rapidly clogged, leading to the stop of the reaction if the initial porosity is below ∼ 30 % (Figure 2). The consumption of fluid by the reaction leads to a decrease in fluid pressure which can reach several tens of megapascals (Figure 2). [START_REF] Evans | Phase-Field Modeling of Reaction-Driven Cracking: Determining Conditions for Extensive Olivine Serpentinization[END_REF] predict a similar behaviour but they consider an additional mechanical feedback between fluid pressure variation and tensile failure. This latter coupling allows porosity to be generated in the fractures as fluid pressure decreases whereas the model proposed here only allows viscoelastic compaction to occur. However, the generated porosity is not sufficient to reach high extents of reaction outside the cracks [START_REF] Evans | Phase-Field Modeling of Reaction-Driven Cracking: Determining Conditions for Extensive Olivine Serpentinization[END_REF]. These model predictions are not consistent with the observations in natural samples in which extents of reaction above 60 % are commonly measured [START_REF] Oufi | Magnetic properties of variably serpentinized abyssal peridotites[END_REF].

The observations of [START_REF] Beard | Onset and progression of serpentinization and magnetite formation in Olivine-rich troctolite from IODP hole U1309D[END_REF], [START_REF] Boudier | Serpentine Mineral Replacements of Natural Olivine and their Seismic Implications: Oceanic Lizardite versus Subduction-Related Antigorite[END_REF] and Malvoisin et al. (2020b) indicate that reaction proceeds locally with mass transfer limited to the submicrometric scale. Positive volume change is thus not accommodated by precipitation, probably due to a slow transport of the aqueous species. The only other way to account for positive volume change in a closed system is to introduce solid volume change in the equation for solid deformation (equation 12). Deriving the change in porosity with this latter equation (last equation in 13) provides an equation without reactive terms, that is without porosity clogging. The modification of equation 12 thus allows for accommodating positive change in solid volume through reaction-induced deformation rather than precipitation in the pores. We solve the system of equations by considering two extreme cases, either no solid volume change (e.g. volume-controlled) or constant total pressure (e.g. stress-controlled).

Reaction is controlled by deformation in our model. If the rock can freely deform by assuming negligible shear stresses (i.e. stress-controlled simulations (ii)), a maximum serpentinisation rate ( ξmax ) is estimated since shear deformation can only delay the volumetric expansion. ξmax depends only on τ and is thus controlled by dissolution and precipitation occurring at the olivine surface (Figures 3 andS5). Complete reaction is achieved in less than 1000 yr, in agreement with previous estimates based on experimental data and reaction-induced fracturing modelling [START_REF] Malvoisin | Control of serpentinisation rate by reaction-induced cracking[END_REF]. The rock expands by approximately 50 % with consequences for surface uplift and large scale deformation at mid-ocean ridges and in subduction zones [START_REF] Germanovich | Deformation and surface uplift associated with serpentinization at mid-ocean ridges and subduction zones[END_REF].

If deformation is not allowed (i.e. volume-controlled simulations (i)), the serpentinisation rate is zero and the reaction does not reach completion (Fig- ures 3 andS4)). We refine the estimation of the minimum reaction rate ( ξmin ) in the following by considering volumetric expansion associated with tectonic deformation which is the main process of deformation at mid-ocean ridges [START_REF] Cannat | Emplacement of mantle rocks in the seafloor at mid-ocean ridges[END_REF]. Tectonic deformation is partitioned between dilatational and shear strains. Microstructural observations indicate that strain in detachment faults is accommodated by crystal-plastic flow and diffusive mass transfer in multiple shear zones [START_REF] Schroeder | Strain localization on an oceanic detachment fault system[END_REF][START_REF] Boschi | Mass transfer and fluid flow during detachment faulting and development of an oceanic core complex[END_REF]. Reaction-softening associated with talc formation promotes shear strain localisation, suggesting that shear deformation accommodates a significant proportion of the deformation. However, this proportion is difficult to quantify since most numerical models of mantle rock exhumation at mid-ocean ridges are incompressible [START_REF] Behn | Magmatic and tectonic extension at mid-ocean ridges: 1. Controls on fault characteristics[END_REF][START_REF] Bickert | How do detachment faults form at ultraslow mid-ocean ridges in a thick axial lithosphere?[END_REF]. The model developed here is one-dimensional and does not allow for shear deformation to be modelled. Moreover, volume increase during reaction can only be accommodated by expansion. As a result, we only consider dilatational strain in the following by introducing a constant rate of expansion in the model (∇ j v j s = ˙ V with ˙ V a constant dilatational strain rate imposed by tectonics). Our model is thus not able to probe all the complex interactions between tectonic deformation and serpentinisation at midocean ridges. However, it allows determining the conditions for which tectonic expansion is slower than volume increase induced by dissolution/precipitation (tectonic-controlled versus dissolution/precipitation-controlled rate).

Simulations performed with several imposed dilatational strain rates indicate that the serpentinisation rate is first equal to a high rate controlled by dissolution and precipitation and rapidly decreases with time (Figure 4). For ˙ V < 10 -13 s -1 and η φ = 10 18 P a.s, the reaction rate decreases to zero as a result of viscoelastic compaction associated with reaction-induced total pressure increase. For ˙ V > 10 -13 s -1 and η φ = 10 18 P a.s, the serpentinisation rate reaches a constant value, indicating, with the formalism adopted in equation 7, that the porosity is constant. The equation for porosity in the system of equations 13 can be solved for a constant dilatational strain rate and a constant porosity to obtain an equation for the minimum serpentinisation rate:

ξmin = - β s (1 -X s ) X seq (1 + β sr ) (β d (1 -φ) -β s ) P -P w (1 -φ) η φ + ˙ V (1 -X s ) X seq (1 + β sr ) (14)
The analytical solution provided by equation 14 fits well with the constant reaction rate achieved in the numerical simulations (Figure 4), showing that the minimum serpentinisation rate is controlled by a competition between viscous compaction (first term of right-hand side of equation 14) and tectonic expansion (second term of right-hand side of equation 14). Figure 5 provides the time necessary to reach an extent of reaction of 80 % as a function of dilatational strain rate and viscosity (see Supplementary material for details). At low viscosity and low strain rate, reaction does not reach completion as viscous compaction dominates. At moderate to high viscosity and ˙ V < 3.10 -10 s -1 , the minimum serpentinisation rate is controlled by tectonic deformation. The time to reach 80 % of reaction linearly decreases with strain rate according to the second term of right-hand side of equation 14. It is in the order of millions or thousands of years for ˙ V = 10 -13 s -1 and ˙ V = 10 -10 s -1 , respectively. At ˙ V > 3.10 -10 s -1 , the serpentinisation rate calculated with equation 14 is higher than the rate calculated with equation 7, suggesting that serpentinisation rate is controlled by dissolution and precipitation at the olivine surface.

Estimating if serpentinisation is controlled by reaction-induced deformation or by tectonics is essential to determine if ξmax or ξmin should be used to model serpentinisation at large scale. Reaction-induced deformation can only occur if the pressure generated during volume expansion is sufficient to induce deformation. The maximum total pressure increase (∆P max ) can be derived from equation 13 as:

∆P max = - 1 + β sr β d log (1 -ξX seq ) (15)
Its values for ξ max and ξ = 1 are 14 MPa and ∼ 500 MPa, respectively, in good agreement with the maximum value calculated in the volume-controlled simulations (14 MPa at ξ max ; Figure S4). Pressure increase during precipitation (i.e. crystallisation pressure; ∆P c ) has been reported in experiments [START_REF] Correns | Growth and dissolution of crystals under linear pressure[END_REF]. It is generally calculated in excess of the confining pressure as [START_REF] Kelemen | Reaction-driven cracking during retrograde metamorphism: Olivine hydration and carbonation[END_REF]:

∆P c = RT V m log Ω ( 16
)
where Ω is the supersaturation in a thin fluid film where precipitation occurs, R is the gas constant and V m is the molar volume of the precipitating crystal. ∆P c values calculated for the serpentinisation reaction can exceed 300 MPa [START_REF] Kelemen | Reaction-driven cracking during retrograde metamorphism: Olivine hydration and carbonation[END_REF], close to the value of ∆P max at ξ = 1. However, equations 15 and 16 are based on different approaches since ∆P max expression considers elastic deformation and depends on the extent of reaction whereas ∆P c equation is based on thermodynamic disequilibrium between the fluid and the solid and requires the fluid to be significantly supersaturated with respect to the precipitating mineral. Experiments of stress monitoring during hydration show that the reaction-induced pressure increase depends on the extent of reaction, that fluid pathways can shut-down as a result of reaction-induced stress, and that the generated overpressures are generally one order of magnitude smaller than the one calculated with equation 16 [START_REF] Wolterbeek | Reaction-driven casing expansion: potential for wellbore leakage mitigation[END_REF][START_REF] Skarbek | Competition Between Crystallization-Induced Expansion and Creep Compaction During Gypsum Formation, and Implications for Serpentinization[END_REF]. The predictions of our model are in agreement with these results with progressive pressure build-up during reaction, compaction and a relatively small total pressure increase (Figure S4). This suggests that equations 13 are a good alternative to crystallisation pressure calculation based on fluid composition (equation 16), particularly when the supersaturation of the solution is difficult to estimate or when other processes such as mechanical compaction have to be considered. Moreover, it indicates that, for the investigated parameters, the overpressure generated during reaction does not exceed 14 MPa due to viscous compaction (Figure S4). This value is below the yield stress for serpentinised peridotites in compression [START_REF] Escartín | Strength of slightly serpentinized peridotites: Implications for the tectonics of oceanic lithosphere[END_REF]. Our model thus predicts a limited driving force for reaction-induced deformation, which could thus only occur in the first several hundreds of metres below the surface. Figure 5 indicates that serpentinisation may be mainly controlled by tectonic deformation for dilatational strain rates below 3.10 -10 s -1 at mid-ocean ridges. This statement is in agreement with the extensive serpentinisation observed near detachments faults [START_REF] Rouméjon | Serpentinization and Fluid Pathways in Tectonically Exhumed Peridotites from the Southwest Indian Ridge (62-65 E)[END_REF]. Moreover, it validates the approach used in some tectonic models [START_REF] Bickert | How do detachment faults form at ultraslow mid-ocean ridges in a thick axial lithosphere?[END_REF] in which serpentinisation rate is related to the strain rate. However, one should keep in mind that shear deformation is not taken into account in Figure 5. We indeed use a constant dilatational strain rate that differs from the shear strain rate computed in previous simulations of mid-ocean ridge processes [START_REF] Bickert | How do detachment faults form at ultraslow mid-ocean ridges in a thick axial lithosphere?[END_REF]. As a result, we do not capture the potential evolution of strain rate during reaction possibly leading to the switch from one controlling mechanism to another. An accurate modelling of the effect of tectonic deformation on serpentinisation rate would require to determine the extent of strain partitioning during exhumation. This should be achieved in the future by solving for the conservation of total momentum in two or three dimensions and by implementing equation 14 in a compressible code. In addition, the model could be improved by considering aqueous species in the fluid to compute the additional thermodynamic force potentially driving stress build-up (equation 16). Viscous compaction is also modelled here with a formalism derived for macropores (Yarushina and Podladchikov, 2015). However, the nanopores preserved during reaction [START_REF] Tutolo | Nanoscale constraints on porosity generation and fluid flow during serpentinization[END_REF] may behave differently leading to a different evolution of the porosity through time.

The numerical approach proposed here is a critical leap forward in reactive transport modelling since it allows for complete reaction without porosity clogging. Serpentinisation but also other hydration and carbonation reactions inducing a positive change in solid volume can be modelled with the approach proposed here. As a first step, we investigate the volume-controlled and stresscontrolled cases but, as shown above, implementing shear stresses will be required in the future to produce quantitative prediction relevant for natural systems. Another improvement of the model concerns the implementation of gas formation (hydrogen and steam). The reaction indeed consumes water, decreasing the fluid pressure to a level sufficient to produce gases at high Da and low De values (Figures S4 B). Vaporisation will in turn significantly modify the fluid density and is thus coupled to fluid pressure evolution. Finally, it would be in-teresting to distinguish between nano-and macropores in the model. [START_REF] Tutolo | Nanoscale constraints on porosity generation and fluid flow during serpentinization[END_REF] indeed show that nanopores are preserved during serpentinisation whereas macropores can potentially be filled by an assemblage of reaction products and nanopores. This aspect should be considered to accurately model the evolution of the hydraulic properties during reaction and reproduce the evolution of permeability observed in experiments of serpentinisation [START_REF] Godard | Incipient hydration of mantle lithosphere at ridges: A reactive-percolation experiment[END_REF][START_REF] Farough | Evolution of fracture permeability of ultramafic rocks undergoing serpentinization at hydrothermal conditions : An experimental study[END_REF].

Rocks

Figure 1 :

 1 Figure1: The two types of models to link reaction, porosity evolution and solid deformation. Mineralogical reactions in the presence of fluid can be decomposed in four steps: fluid transport in the pore network (1); primary mineral (in yellow) dissolution (2); aqueous species transport (3); and secondary phase (in green) precipitation (4). A: the classical approach consists in considering that aqueous species are transported up to the pores where precipitation occurs. This process does not generate deformation directly. B: the new approach proposed here is based on observations at the nanoscale showing that aqueous species transport can be limited to a scale of less than 100 nm. As a result, the reaction products do not precipitate in the pores but at the surface of the primary mineral. This implies that reaction induces rock deformation (i.e. swelling; red arrows).

FFigure 3 :

 3 Figure3: Results of the reaction-induced deformation numerical model (equations 7 and 13) when assuming a constant solid volume (A, B and C) or a constant total pressure (D, E and F). The evolutions with time of extent of reaction (A and D), fluid pressure (plain line) and total pressure (dashed line) (B and E), and porosity (C and F) are displayed. The blue line in subpanels A and D corresponds to the maximum extent of reaction predicted for a nondeforming solid matrix (ξmax; see text for details). The model is run with k φ = 10 -19 m 2 , L = 10m, τ = 5.10 7 s and η φ = 10 18 P a.s

Figure 4 :

 4 Figure4: Serpentinisation rate ( ξ) as a function of time for the model considering reactioninduced deformation and an imposed tectonic expansion. The plain black lines corresponds to five simulations performed with k φ = 10 -19 m 2 , L = 10m, τ = 5.10 7 s, η φ = 10 18 P a.s, and at the indicated dilatational strain rates ( ˙ V ). The blue dashed line corresponds to the rate imposed by dissolution and precipitation. The black dashed lines corresponds to the rate calculated with equation 14 considering tectonic deformation and viscous compaction. The red dashed lines corresponds to the rate associated with tectonic deformation only and calculated with the second term of right-hand side of equation 14. Note that the serpentinisation rate tends towards the rate imposed by tectonic deformation for ˙ V > 10 -13 s -1 .

Figure 5 :

 5 Figure 5: Time necessary to reach 80 % of extent of reaction (t 80 in years) as a function of viscosity (η φ ) and dilatational strain rate ( ˙ V ). t 80 is calculated with equation 14. Three domains can be distinguished: a domain with no reaction due to viscous compaction leading to fluid pathways clogging, tectonic-controlled domain in which the serpentinisation rate can be expressed with the second term of the right-handed side of equation 14, and a dissolutionand precipitation controlled domain in which the rate determined with equation 14 is faster than the one determined with equation 7.

  With Implications for Serpentinization. Journal of Geophysical Research: Solid Earth 124, 5708-5733. URL: https://onlinelibrary.wiley. com/doi/abs/10.1029/2018JB016872, doi:10.1029/2018JB016872.

Acknowledgements B.M. acknowledges support from the Swiss National Science Foundation (Grant No. PZ00P2 168083). Y.P. and A.M. were supported by the Russian Ministry of Science and Higher Education (project No. 075-15-2019-1890). We thank Luc Lavier and Lars Rüpke for their thoughtful comments and Jean-Philippe Avouac for editorial handling. B.M. had the idea for the reactioninduced swelling approach aboard Chikyu during the Phase 2 Shipboard Science Party of the Oman Drilling Project. He thanks the members of this Science Party for the stimulating atmosphere aboard.


Thiago C Bellardi
email: thiago.craesmeyer-bellardi@inrialpes.frchristian.laugier@inrialpes.fr

Dizan Vasquez
email: vasquez@mavt.ethz.ch

Christian Laugier

C Bellardi

Frame Rate Object Extraction from Video Sequences with Self Organizing Networks and Statistical Background Detection

In many computer vision related applications it is necessary to distinguish between the background of an image and the objects that are contained in it. This is a difficult problem because of the double constraint on the available time and the computational cost of robust object extraction algorithms.

This paper builds upon former work on combining the strong theoretical foundations of clustering with the speed of other approaches. It is based on a novel Self Organizing Network (SON) which has a robust initialization schema and is able to find the number of objects in an image or grid. The main contribution of our extension is that it eliminates the use of a threshold, allowing the algorithm to work on continuous, while having a complexity that remains linear with respect to the number of pixels or cells.

I. INTRODUCTION

Image segmentation is an important and challenging problem in vision. Its goal is to identify homogeneous regions in images as distinct from the background and belonging to different objects. A common approach is to classify pixels on the basis of local features (e.g. color, position, texture), and then group them together according to their class in order to identify different objects in the scene.

For the specific problem of finding moving objects from static cameras, the traditional segmentation approach is to separate pixels into two classes: background and foreground. This is called Background Subtraction [START_REF] Koller | Towards robust automatic traffic scene analsis in realtime[END_REF] and constitutes an active research domain, the interested reader is referred to [START_REF] Piccardi | Background subtraction techniques: a review[END_REF] for an overview of the field's state of the art.

The output of most background segmentation techniques consists of a bitmap image, where values of 0 and 1 correspond to background and foreground, respectively (eg [START_REF] Friedman | Image segmentation in video sequences: A probabilistic approach[END_REF], [START_REF] Stauffer | Learning patterns of activity using real-time tracking[END_REF], [START_REF] Cucchiara | Detecting moving objects, ghosts, and shadows in video streams[END_REF]). Having such a bitmap, the next processing step consists of merging foreground pixels to form bigger groups corresponding to candidate objects, this process is known as object extraction.

One common procedure to perform object extraction consists of finding 4 or 8-connected components. This is done using efficient algorithms whose time complexity is linear with respect to the number of pixels in the bitmap [START_REF] Suzuki | Fast connected-component labeling based on sequential local operations in the course of forwardraster scan followed by backward-raster scan[END_REF], [START_REF] Chang | A linear-time component-labeling algorithm using contour tracing technique[END_REF].

A problem with this approach is that it usually produces many small regions which may correspond to noise or to larger regions which failed to merge.

One approach to dealing with this situation is to filter out regions composed of less than a given number of pixels [START_REF] Chen | Object segmentation for video coding[END_REF]. Although this approach is fast, it has the drawback of assuming that all small regions are noise, which, in many situations, is clearly not the case. A second approach consists of relaxing the neighborhood criterion by assuming, for example, that regions separated by one background pixel are still connected. The usual way of doing this is by preprocessing the bitmap image using morphological operators (eg dilation, closing), which have the effect of "thickening" the pixels and "filling in" the holes [START_REF] Meyer | Morphological segmentation[END_REF]. Two problems with this approach are the difficulty of finding the appropriate parameters for the operators and the lack of clear physical interpretation of the operators' parameters. A third approach to object extraction is the use of clustering techniques to group pixels. This opens up the possibility of choosing between a plethora [START_REF] Jain | Data clustering: A review[END_REF] of different algorithms having well understood theoretical properties. On the other hand, most of the robust clustering algorithms (eg [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF], [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]) have three problems when applied to object extraction: a) the number of objects to be found should be known beforehand, b) the algorithms' performance is strongly dependent on the initialization and c) most algorithms are just too complex to be used in systems subject to demanding real-time constraints.

In previous papers, we have proposed a novel clustering approach for object extraction based on Self Organizing Networks (SON). We have applied this algorithm to images [START_REF] Vasquez | A novel self organizing network to perform fast moving object extraction from video streams[END_REF] and occupancy grids [START_REF] Vasquez | Fast object extraction from bayesian occupancy grids using self organizing networks[END_REF], and shown that it is able to produce good results in real time. However, this approach suffered from two drawbacks: a) first, it needed a binary input, which implied the necessity of applying a threshold to the input grid, thus losing information; and b) its complexity depended on the number of input cells and the size of the SON, thus, when the number of input cells was high, the performance of the algorithm degraded. This paper presents an improved version of the algorithm which is able to process continuous input values and, whose complexity only depends on the size of the input grid.

The rest of this paper is structured as follows: in the next section, we present a general approach to object extraction using the k-means algorithm and then present the details of the improved algorithm. Sec. III presents an overview of the statistical background/foreground classifier we used. In section IV we explain our implementation of the approach and present some preliminary results against real and synthetic data. Finally our conclusions and some further research directions are presented in §V.

II. CLUSTERING-BASED OBJECT EXTRACTION

Assuming that the number of objects k in a bitmap is known, applying clustering to object extraction using the kmeans (i.e. Expectation-Maximization) [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF], [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] algorithm is relatively straightforward: 1) Initialize k cluster centers µ i with arbitrary values. 2) Assign each foreground pixel to its closest cluster center. 3) Re-estimate every cluster center µ i as the mean of the points allocated to that cluster. 4) Repeat steps 2-4 until some convergence criterion is met (e.g. minimal cluster reassignment).

However, in most cases, the value of k is unknown. Furthermore, even knowing k, the quality of the obtained clustering depends heavily on initialization, since the algorithm trends to get stuck in local minima. Finally every iteration has a cost of O(N f k) (where N f is the number of foreground pixels) and, sometimes, many iterations are needed before converging.

In order to deal with those problems, this paper proposes an object extraction approach which combines a Selforganizing Network inspired by the Growing Neural Gas [START_REF] Fritzke | A growing neural gas network learns topologies[END_REF] combined with a graph theoretic algorithm used to cut edges in the network's graph.

A. SON-based object extraction

The network is built from M = W × H nodes connected with undirected edges, arranged in a grid with H rows and W columns (fig. 1(a)). This means that, with the exception of nodes located in the borders, every node i will be connected to four other nodes or neighbors (neigh(i)), individually denoted by u(i), d(i), r(i) and l(i) for up, down, right and left, respectively. Every node i has two associated variables: its mean value µ i = (x i , y i ) and an accumulator c i . In a similar manner, for every edge connecting nodes i and j there will be an accumulator e i, j . Besides W and H, the algorithm has two other parameters: 0 < ε n < ε w ≤ 1 which are the learning rates for node mean adaption.

The following subsections describe the steps that our algorithm performs for every video frame, using the bitmap image produced by foreground/background classification as an input.

1) Initialization: The network is initialized by assigning values to all the µ i node centers in order to form a regular grid (fig. 1(a)). Also, the values of all the weights are set to zero:

{c i ← 0, e i, j ← 0 ∀ i, j | i ∈ [1, M], j ∈ neigh(i)} (1)
2) Learning: The learning stage takes as input a bitmap I, typically representing the difference between the last image comming from the camera and a background model (fig. 1(b)). Pixels from the input image, are processed starting from the upper-left corner and then sweeping every row from left to right. For every pixel i, its coordinates p i and value I(p i ) are used to update the SON in four steps: a. Find the node whose mean value is closest to p i , referenced as winner (w i ). The search is restricted to a subset of nodes surrounding the previous winner w i-1 (eg the one corresponding to the previous processed pixel). This subset, that we call the search boundary, is represented by sbound(i -1) (see fig. 2).

w i = argmin j∈sbound(i-1) p -µ j -1 (2)
Fig. 2. Neighborhood and search boundary relative to node i b. Look in neigh(w i ) for the second nearest node to p i :

s i = argmin j∈neigh(w i ) p i -µ j (3)
c. Increment the accumulators e w i ,s i and c w i by the pixel value I(p i ):

e w i ,s i ← e w i ,s i + I(p i ) (4)
and

c w i ← c w i + I(p i ) (5)
d. Adapt the mean of w i and all his neighbors:

µ w i ← µ w i + ε w I(p i ) c w i (p i -µ w i ) (6)
µ j ← µ j + ε n I(p i ) c j (p i -µ j ) ∀ j ∈ neigh(w i ) (7)
3) Relabeling nodes: As a result of the learning step, the network adapts its form to represent the objects in the bitmap (fig. 1(c)). The last step of our algorithm, identifies individual objects by assigning a discrete value to every node in the SON, so that nodes having the same label belong to the same node.

Our algorithm finds groups of nodes by merging nodes according to the weight of their common edges e i, j . The idea is that a higher value of e i, j corresponds to a higher likelihood that nodes i and j belong to the same object. Under this assumption, it is possible to compute a maximum likelihood estimation of the probability, denoted by P i, j , that two nodes "belong together" by using the Laplace law of succession 1 :

P i, j = e i, j + 1 V + (W -1)H + (H -1)W (8)
where

V = ∑ p∈I I(p) (9)
Also by using the Laplace law of succession, we calculate the value of the uniform link probability distribution, which may be seen as the maximum entropy estimate of P i, j prior to learning.

U links = 1 (W -1)H + (H -1)W (10)
In a similar fashion, the weight c i is an indicator of the likelihood that node i belongs to an object, it permits us to evaluate how likely it is that the node belongs to the background instead of an object. This is formulated as a probability P i .

P i = c i + 1 V +W H ( 11
)
1 P i, j is a notational shortcut introduced for the sake of readability, being rigorous it should be written as

P([O i = m] | [O j = m]
), where all the O i variables are binary and O i = m indicates that node i has been assigned to cluster m. A similar shortcut has been used with P i for the same reasons.

With the corresponding uniform being:

U nodes = 1 W H (12)
We use a conventional scanning algorithm to relabel the nodes. The only particularity of our approach is that we use P i, j as the region-merging criterion instead of using colors or other features. Here, we will outline the labeling algorithm, however, the presentation of the complete implementation details is beyond the scope of this paper. The reader is referred to [START_REF] Suzuki | Fast connected-component labeling based on sequential local operations in the course of forwardraster scan followed by backward-raster scan[END_REF], [START_REF] Chang | A linear-time component-labeling algorithm using contour tracing technique[END_REF] for efficient linear-time ways to implement the algorithm.

The algorithm starts from the upper-left node and proceeds by scanning from left to right and from top to bottom, for every node i the following steps are applied: a. Assign the label ∞ to i. b. If P i,l(i) > U links , assign to i the label of l(i) (merge with left region). c. If P i,u(i) > U links , assign to i the minimum between its current label and the label of u(i). Let a be that minimal label and let b be the label of u(i). Relabel nodes on the previous rows having label b to a (merge with upper region). d. If i's label is ∞ assign the next unused label to i (create a new region). At this point, it is important to highlight the fact that the labels obtained above are just identifiers used to distinguish one region (ie object) from the other and that new labels are obtained by incrementing a counter.

4) Computing cluster representations:

Having labeled the nodes, a the probability that a pixel, given by its image coordinates p belongs to a cluster m may be represented using a gaussian distribution 2 :

P * (p | m) = N (p; µ * m , S * m ) (13)
The cluster's prior may be used to filter out clusters whose prior is below a given threshold, it is computed as:

P * m = ∑ i∈m P i ( 14
)
Its mean value,

µ * m = 1 P * m ∑ i∈m P i µ i ( 15
)
And its covariance,

S * m = ∑ i∈m P i P * m (x i -x * m ) 2 (x i -x * m )(y i -y * m ) (x i -x * m )(y i -y * m ) (y i -y * m ) 2 (16)
Alternatively, a cluster may also be viewed as a mixture of gaussians, corresponding to individual nodes in the cluster: 2 Hereafter, cluster parameters will be denoted by a superscript asterisk, in order to distinguish them from node parameters In order to compute the covariance matrices S i , we use the points located halfway between i and its neighbors (fig. 3):

P * (p | m) = ∑ i∈m P i N (p; µ i , S i ) (17)
S i = ∑ j∈neigh(i) P j K    x j +x i 2 2 (x j +x i )(y j +y j ) 4 (x j +x i )(y j +y i ) 4 y j +y i 2 2    (18)
Where K = ∑ j∈neigh(i) P j is a normalization constant.

Fig. 3. Estimating the covariance, represented by the ellipse, from the midpoints between a node (center) and its neighbors.

In cases where the algorithm is required to produce interest regions it is often convenient to produce bounding boxes which are slightly larger than the contained object. We have computed the size of these regions using the difference between the maximum and the minimum mean values of the cluster nodes as they were before learning 3 this may be regarded as finding the area bounded by nodes which have not been adapted.

B. Complexity Analysis

The changes in the search algorithm, denoted by the equations 2 and 3, represent a complexity change from O(N f M) to O(N), where N corresponds to the total number of pixels in the image. Noting that now it takes into account all the pixels in the image, instead of just the ones marked as being part of the foreground, and its complexity is independent of the SON grid size. This allows us to weight the pixels background/foreground contribution in a continuous fashion, instead of using a hard threshold, which makes the approach much more robust. The key idea for complexity independence on the SON size is to exploit the fact that the network is processed in a top-bottom, left-right sequence, and to limit the set of nodes in the SON which need to be compared with each pixel by looking in the neighborhood of the last processed pixel.

Thanks to the existence of efficient algorithms, the cost of labeling is linear with respect to the number of nodes in the SON, moreover, the computation of the cluster representation (i.e. gaussian parameters, mixture of gaussian parameters and bounding boxes) may be performed at the same time as labeling. Thus, the algorithm's overall complexity is O(N).

III. STATISTICAL BACKGROUND/FOREGROUND

CLASSIFIER

The goal of the statistical background classifier is to learn a model of the scene background from a sequence of 3 This idea was suggested by David Raulo images. It represents each pixel as an independent normal distribution. For every pixel, the classifier outputs a dissimilarity measure, indicating how different the pixel is from the modeled background, and consequently, how likely it is that pixel belongs to the foreground.

Let us have a given input pixel whose coordinates in the image frame are given by p = [x, y] and whose intensity value is given by I(p). The corresponding background normal will have a mean and standard deviation denoted as µ p and σ p respectively. The dissimilarity measure D(p = f g) is computed using the Mahalanobis distance between I(p) and N(µ p , σ 2 p ):

D(p = f g) = |µ p -I(p)| σ 2 p ( 19
)
Once the measure D(p = f g) is computed, the background is updated as follows:

µ p = D(p = f g)µ p + (1 -D(p = f g))I(p) ( 20
)
σ p = D(p = f g)σ p + (1 -L(p = f g))|I(p) -µ p | (21) µ p ← αµ p + (1 -α) µ p (22) σ p ← ασ p + (1 -α) σ p ( 23
)
where α controls the learning rate.

IV. EXPERIMENTAL RESULTS

The algorithm was tested with three different inputs: (a) A binary bitmap obtained by thresholding the absolute difference between the intensity level of the current and previous video frames; (b) A grayscale bitmap obtained from the statistical background/foreground classifier described above, where the intensity represents the mahalanobis distance between the actual frame and the gaussian estimated background (see III); and (c) a thresholded version of the frame used in (b). The SON size and learning factors are the same for the three cases: 20x20 nodes, ε w = 0.1 and ε n = 0.01. The tests were conducted using the CAVIAR test case scenarios [START_REF]Datasets and videos of the european project caviar[END_REF], which consist of a number of video sequences of people moving in the INRIA Lab's entry hall. The videos come with data files containing the ground truth of the sequences, which has been obtained by hand-labeling the images. A typical image of our detector running on one of these videos is shown in fig. 4.

For each frame processed the result of the extraction is compared with the ground truth and the following parameters are computed:

The mean values for these parameters, obtained from 1042 effective frames processed from the CAVIAR [START_REF]Datasets and videos of the european project caviar[END_REF] The obtained result shows that the detector can be combined with more sophisticated foreground classification methods to improve the quality of detection. Is clear that the better the background/foreground classification is, the better will be the result from the clustering, but what is interesting in the comparison between the results from the inputs (b) and (c) is that the method used to extract object performs better with the non-thresholded input.

In respect to computational costs, the object extraction from a 384×288 grayscale image, using a 20×20 SON consumes less than 50msec (mean obtained from 1044 frames processed), running on a machine equipped with an Intel Pentium IV processor. To illustrate the algorithm complexity property, we repeated the procedure above using a 100 × 100 SON and the increase in the processing time kept around 1msec.

V. CONCLUSIONS AND FUTURE WORK

In this paper we have discussed object extraction from continuous valued bitmaps emphasizing the advantages, but also the three big problems of cluster based algorithms (ie need to know the number of objects to be detected beforehand, sensibility to initialization and complexity) and extended previous work on a Self Organizing Network based on the Growing Neural Gas algorithm which solves the above mentioned problems and keeps the strong theoretical properties of clustering algorithms. Our extension permits makes the complexity of the algorithm independent of the size of the underlying SON, and eliminates the need of obtaining a binary image through a threshold.

We have explained the details of our algorithm, and shown how it may be used to find clusters and represent them using gaussians, mixtures of gaussians or bounding boxes.

Finally, we have discussed the experimental results we have obtained by comparing our approach to a ground truth consisting of hand-labeled data. Our results seem to confirm that our approach is fast, robust and general. We have actually applied the approach to find regions of interest that are fed to a pedestrian detection framework using fixed cameras, and obtained very satisfactory results.

Future work includes continuing our experimental work, in particular by improving the background/foreground classifier. Other possibilities include taking into account temporal information by updating the state of the SON when a new input image is available, instead of working in a frame per frame fashion. Finally, we would like to explore the use use of our SON to perform data fusion on a multicamera system installed in a parking lot.

Fig. 1 .

 1 Fig. 1. Approach overview. Enlarged views showing the different steps of our algorithm using CAVIAR data.

1 )Fig. 4 . 4 )

 144 Fig. 4. Object detection using a statistical background/foreground classifier. The grayscale foreground is passed as it to the object detector.

This work has been partially supported by the european BACS Project. T.