Qingkai Fan

Jean-Arthur Olive

Mathilde Cannat

Thermo-Mechanical State of Ultraslow-Spreading Ridges with a Transient Magma Supply

New model for the thermal structure of ultraslow ridges couples hydrothermal circulation with episodic sill emplacement • Hydrothermal convection in uppermost 6 km of lithosphere with thick boundary layer can explain cold thermal regime • High-output venting at ultraslow ridges likely requires transient shallow magmatic intrusions.

Introduction

New oceanic lithosphere continuously forms at the axis of mid-ocean ridges (MORs) through a combination of magmatic and tectonic processes. The thermal structure of the ridge axis is a critical control on the mechanical properties of young lithosphere, and thus modulates faulting styles and crust emplacement [START_REF] Tapponnier | Necking of the lithosphere and the mechanics of slowly accreting plate boundaries[END_REF][START_REF] Chen | A nonlinear rheology model for mid-ocean ridge axis topography[END_REF][START_REF] Shaw | Causes and consequences of variations in faulting style at the Mid-Atlantic Ridge[END_REF][START_REF] Behn | Effects of hydrothermal cooling and magma injection on mid-ocean ridge temperature structure, deformation, and axial morphology[END_REF][START_REF] Buck | Modes of faulting at mid-ocean ridges[END_REF]. The ~650˚C isotherm, in particular, roughly coincides with the brittle-ductile transition in mafic and ultramafic lithologies (e.g., [START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF][START_REF] Mackwell | High-temperature deformation of dry diabase with application to tectonics on Venus[END_REF]. Its depth is therefore expected to mark the base of the brittle, seismogenic part of the lithosphere. Likewise, the depth to the ~1000˚C isotherm (solidus of basaltic magmas) likely corresponds to the shallowest depth at which magma can pool and reside over long lime scales [START_REF] Cannat | Emplacement of mantle rocks in the seafloor at mid-ocean ridges[END_REF]. In a broad sense, the thermal structure of a MOR reflects a balance between heat advected upwards from the asthenosphere --through mantle upwelling and magma ascent--and heat transferred through the lithosphere by hydrothermal circulation [START_REF] Lin | A finite amplitude necking model of rifting in brittle lithosphere[END_REF][START_REF] Chen | A nonlinear rheology model for mid-ocean ridge axis topography[END_REF][START_REF] Phipps-Morgan | The genesis of oceanic crust: Magma injection, hydrothermal circulation, and crustal flow[END_REF][START_REF] Montési | Mantle flow and melting underneath oblique and ultraslow mid-ocean ridges[END_REF][START_REF] Liu | Models of hydrothermal heat output from a convecting, crystallizing, replenished magma chamber beneath an oceanic spreading center[END_REF][START_REF] Theissen-Krah | Coupled mechanical and hydrothermal modeling of crustal accretion at intermediate to fast spreading ridges[END_REF][START_REF] Lowell | Characteristics of magma-driven hydrothermal systems at oceanic spreading centers[END_REF].

At fast and intermediate MORs, as well as magmatically robust sections of slow MORs, the input of magmatic heat is to first order steady over geological time. One or several sills are continuously replenished or frequently emplaced within the lower oceanic crust, and feed episodic dike intrusions into the upper crust [START_REF] Detrick | Multi-channel seismic imaging of a crustal magma chamber along the East Pacific Rise[END_REF][START_REF] Sinton | Mid-ocean ridge magma chambers[END_REF][START_REF] Kelemen | Faulting and Magmatism at Mid-Ocean Ridges[END_REF][START_REF] Marjanović | A multi-sill magma plumbing system beneath the axis of the East Pacific Rise[END_REF]. Focused intrusions can in turn fracture the crust, rendering it permeable to fluid flow [START_REF] Wilcock | The role of magma injection in localizing black-smoker activity[END_REF]. Heat and permeability enable vigorous hydrothermal circulation above and along the flanks of the shallowest axial melt lens [START_REF] Tolstoy | Seismic identification of along-axis hydrothermal flow on the East Pacific Rise[END_REF][START_REF] Hasenclever | Hybrid shallow on-axis and deep off-axis hydrothermal circulation at fast-spreading ridges[END_REF][START_REF] Fontaine | Along-axis hydrothermal flow at the axis of slow spreading Mid-Ocean Ridges: Insights from numerical models of the Lucky Strike vent field (MAR)[END_REF], which can transfer magmatic heat out of the lithosphere efficiently in quasi-steady state, as long as the conductive boundary layer separating the magma from the permeable zone is thin (100 m, [START_REF] Liu | Models of hydrothermal heat output from a convecting, crystallizing, replenished magma chamber beneath an oceanic spreading center[END_REF]. Models that balance magmatic input and hydrothermal output are generally successful at explaining the thermal state of magmatically-robust MORs, as they accurately predict the depth to the shallowest axial melt lens at intermediate-to fast-spreading ridges (e.g., Phipps-Morgan and [START_REF] Phipps-Morgan | The genesis of oceanic crust: Magma injection, hydrothermal circulation, and crustal flow[END_REF][START_REF] Theissen-Krah | Coupled mechanical and hydrothermal modeling of crustal accretion at intermediate to fast spreading ridges[END_REF].

At the axis of intermediate-to fast-spreading ridges, the brittle lithosphere is typically thinner than ~5 km and hydrothermal convection likely occurs in close proximity to a steady (or frequently replenished) magmatic heat source [START_REF] Maclennan | Thermal models of oceanic crustal accretion: Linking geophysical, geological and petrological observations[END_REF][START_REF] Gillis | The roof of the axial magma chamber: A hornfelsic heat exchanger[END_REF][START_REF] Liu | Models of hydrothermal heat output from a convecting, crystallizing, replenished magma chamber beneath an oceanic spreading center[END_REF][START_REF] Machetel | Numerical model of crustal accretion and cooling rates of fast-spreading mid-ocean ridges[END_REF][START_REF] Theissen-Krah | Modes of crustal accretion and their implications for hydrothermal circulation[END_REF]. This paradigm however breaks down at certain sections of slow and ultraslow MORs where recent seismological studies have revealed the occurrence of brittle deformation at depths below 10-15 km [START_REF] Schlindwein | Mid-ocean-ridge seismicity reveals extreme types of ocean lithosphere[END_REF][START_REF] Grevemeyer | Constraining the maximum depth of brittle deformation at slowand ultraslow-spreading ridges using microseismicity[END_REF]. For example, the "smooth seafloor" section of the Southwest Indian Ridge (SWIR) at 64˚40'E [START_REF] Sauter | The ultraslow spreading Southwest Indian ridge, in Diversity of Hydrothermal Systems on Slow Spreading Ocean Ridges[END_REF][START_REF] Sauter | Continuous exhumation of mantle-derived rocks at the Southwest Indian Ridge for 11 million years[END_REF], which has an effective (i.e., cross-axis) spreading rate of 14 mm/yr features microseismicity down to at least 15 km below seafloor (Fig. 1b) (Cannat et al., 2019b). Such unusually thick axial lithosphere cannot be a result of the ultraslow spreading rate alone [START_REF] Montési | Mantle flow and melting underneath oblique and ultraslow mid-ocean ridges[END_REF], as the Mid-Cayman Spreading Center (MCSC), for example, spreads orthogonally at a similar rate (~15 mm/yr) yet only hosts microearthquakes within 10 km below seafloor near the Mt. Dent detachment fault (Fig. 1a) [START_REF] Grevemeyer | Constraining the maximum depth of brittle deformation at slowand ultraslow-spreading ridges using microseismicity[END_REF]. The difference may stem from the local magma supply, which can be decorrelated from spreading rate if ascending melts undergo strong along-axis focusing [START_REF] Cannat | Emplacement of mantle rocks in the seafloor at mid-ocean ridges[END_REF][START_REF] Dick | An ultraslow-spreading class of ocean ridge[END_REF]. Indeed, the MCSC presently accretes ~4km thick (seismic) crust (von Avendonk et al., 2017), while the SWIR at 64º40'E features virtually no magmatic crust and exposes large tracts of serpentinized mantle at the seafloor. The magma supply of ultraslow MORs is indeed thought to be highly transient in space and time, and often does not permit the accretion of a continuous crustal section [START_REF] Cannat | Emplacement of mantle rocks in the seafloor at mid-ocean ridges[END_REF][START_REF] Dick | An ultraslow-spreading class of ocean ridge[END_REF][START_REF] Standish | Young off-axis volcanism along the ultraslow-spreading Southwest Indian Ridge[END_REF][START_REF] Sauter | Continuous exhumation of mantle-derived rocks at the Southwest Indian Ridge for 11 million years[END_REF]Cannat et al., 2019a). The thermal state of ultraslow ridges is therefore unlikely to reflect the transfer of heat from a continuously replenished magma plumbing system to a hydrothermal convection system. This said, magma-poor, ultraslow ridge sections are not devoid of hydrothermal circulation, as water column anomalies suggesting active venting have been identified above several magma-poor sections of the SWIR [START_REF] Baker | Hydrothermal venting in magma deserts: The ultraslow-spreading Gakkel and Southwest Indian Ridges[END_REF]. Both low-and high-temperature hydrothermal vents were later discovered along the SWIR axis, often in close proximity to detachment faults. The low-temperature Old City field (64ºE, Fig. 1b), for example, sits on the footwall of an active fault primarily exposing serpentinized mantle (Cannat et al., 2019b), while the high-temperature Longqi field (49ºE, Fig. 1c,d) sits on a detachment which clearly dissects a magmatically robust axis, as evidenced by the largely basaltic, hummocky seafloor that surrounds it [START_REF] Zhao | Three-dimensional seismic structure of the Dragon Flag oceanic core complex at the ultraslow spreading Southwest Indian Ridge (49˚39'E)[END_REF][START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF]. Interestingly, microseismicity data suggests that both Old City and Longqi lie above ~15-km thick brittle lithosphere [START_REF] Yu | Lithospheric Structure and Tectonic Processes Constrained by Microearthquake Activity at the Central Ultraslow-Spreading Southwest Indian Ridge[END_REF][START_REF] Chen | Microseismicity constraints on brittle lithosphere thickness at a nearly amagmatic section of the ultraslow spreading Southwest Indian Ridge[END_REF]. This illustrates the complex relationship between the thermal state of the ridge axis and the vigor of hydrothermal convection. [START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF] proposed that the great thermal output of Longqi reflects vigorous hydrothermal convection all the way down to a magmatic heat source ~15 km below seafloor, with fluid upwelling primarily occurring within the detachment fault (Fig. 1d). This is potentially problematic as confining pressure is known to reduce pore space and hinder the formation of well-connected fracture networks [START_REF] Fisher | Permeability within basaltic oceanic crust[END_REF]. In fact, extrapolating permeability measurements from the uppermost 1.5 km of the oceanic crust down to depths >5 km yields values well below 10 -18 m 2 [START_REF] Kuang | An integrated permeability-depth model for Earth's crust[END_REF], which amount to subcritical Rayleigh numbers incompatible with deep hydrothermal convection [START_REF] Lowell | Hydrothermal processes at mid-ocean ridges: Results from scale analysis and single-pass models[END_REF]. While greater permeability is plausible within the damage zone of detachment faults (e.g., [START_REF] Scibek | Multidisciplinary database of permeability of fault zones and surrounding protolith rocks at world-wide sites[END_REF], an average permeability of ~510 -16 m 2 in the bulk of the lithosphere down to 15 km depth --as required for off-fault hydrothermal recharge in the model of [START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF]--exceeds predictions of the [START_REF] Kuang | An integrated permeability-depth model for Earth's crust[END_REF] model by over 4 orders of magnitude. Furthermore, seismic evidence for hydration of the oceanic lithosphere (outside of subduction zone settings) is typically restricted to a few km below seafloor. For example, the velocity of P waves reaches 8 km/s at ~5 km depth beneath the detachment-dominated MARK area of the slow-spreading Mid-Atlantic Ridge, suggesting the underlying peridotites are unaltered. The same is true in 125 Ma oceanic crust south of Bermuda, which formed at the Mid-Atlantic Ridge during an episode of ultraslow spreading [START_REF] Lizarralde | Spreading-rate dependence of melt extraction at mid-ocean ridges from mantle seismic refraction data[END_REF]. Interestingly, recent seismic imaging of the West Iberia amagmatic passive margin suggests that faults act as pathways for fluids to reach the upper mantle beneath thinned continental crust, but that mantle hydration is confined to within ~6 km below seafloor [START_REF] Bayrakci | Fault-controlled hydration of the upper mantle during continental rifting[END_REF]. This maximum depth is consistent with the model put forward by deMartin et al. ( 2004), in which the extent of water percolation in olivine-dominated lithosphere is controlled by its ability to grow interconnected micro-fissures in response to cooling stresses. On the other hand, [START_REF] Schlindwein | Mid-ocean-ridge seismicity reveals extreme types of ocean lithosphere[END_REF] identified extensive regions of young ultraslow-spreading lithosphere deforming aseismically down to ~15 km, and interpreted them as extensively serpentinized, implying deep water circulation. As the lack of seismicity may not be a reliable indicator of serpentinization (e.g., [START_REF] Moore | Frictional strengths of talc-serpentine and talc -quartz mixtures[END_REF], the depth extent of hydrothermal circulation and related hydration very much remains an open question.

With this uncertainty in mind, two conceptual models can be envisaged for the thermal structure of ultraslow ridges. In the first one, permeability --perhaps enhanced by tectonic activity-enables hydrothermal convection over the entire thickness of the brittle lithosphere, which can be as large as ~15 km. This essentially amounts to extending the standard thermal model of intermediate / fast MOR axes down to ~3 times its usual depth range (e.g., Phipps-Morgan and [START_REF] Phipps-Morgan | The genesis of oceanic crust: Magma injection, hydrothermal circulation, and crustal flow[END_REF][START_REF] Liu | Models of hydrothermal heat output from a convecting, crystallizing, replenished magma chamber beneath an oceanic spreading center[END_REF][START_REF] Theissen-Krah | Coupled mechanical and hydrothermal modeling of crustal accretion at intermediate to fast spreading ridges[END_REF]. Because the efficiency of hydrothermal convection at transferring heat (i.e., its Nusselt number) does not depend on the depth of the convective system but solely on the thickness of its conductive boundary layer [START_REF] Lowell | Hydrothermal processes at mid-ocean ridges: Results from scale analysis and single-pass models[END_REF][START_REF] Olive | Smoke without fire: How long can thermal cracking sustain hydrothermal circulation in the absence of magmatic heat?[END_REF], high-temperature venting is possible over a cold, thick lithosphere as long as fluids can circulate in close proximity to a deep magmatic heat source [START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF]. We instead choose to focus on the alternative model in which hydrothermal convection is confined to shallower portions of the lithosphere due to confining pressure preventing the opening of interconnected fractures. If large extents of lithosphere above the magma emplacement zone remain impermeable, convecting fluids can never be in close proximity to the melt, and the resulting thermal structure of the ridge will reflect the balance between a transient magmatic heat supply and the rate of heat conduction through a thick boundary layer, in a manner that has not yet been quantified.

We approach this problem through numerical and analytical geodynamic modeling, by coupling a hydrothermal convection model where the maximum depth of percolation evolves with the temperature field in a self-consistent manner. Convection is fueled by the periodic addition of a small increment of specific and latent heat representing discrete sill intrusions. We run simulations until a stationary state is attained and relate the properties of this steady thermal regime to various seafloor and geophysical observables. We then explore possible transient deviations from this steady state.

2 Numerical model of hydrothermal convection with a transient magmatic heat source

Model setup

We use the finite-volume method of [START_REF] Crone | Modeling the effects of tidal loading on mid-ocean ridge hydrothermal systems[END_REF] and [START_REF] Olive | Smoke without fire: How long can thermal cracking sustain hydrothermal circulation in the absence of magmatic heat?[END_REF] to solve for conservation of mass during flow of an incompressible (Boussinesq) Darcian fluid in a saturated, isotropic and non-deforming porous medium representing young oceanic lithosphere (Fig. 2):

𝛁 • (- 𝑘 𝜇 𝑓 (𝛁𝑷 -𝜌 𝑓 𝒈)) = 0, ( 1
)
where k is permeability and f and f denote the dynamic viscosity and density of the hydrothermal fluid, respectively. The quantity inside the divergence operator corresponds to the velocity of the fluid vf. Bold quantities indicate vectors, e.g., the gradient of pressure P, and the acceleration of gravity g. Equation ( 1) is solved jointly with the heat equation to obtain the evolution of the temperature field T:

𝜌 𝑒𝑓𝑓 𝐶 𝑒𝑓𝑓 𝜕𝑇 𝜕𝑡 + 𝜌 𝑓 𝐶 𝑓 𝒗 𝒇 • 𝛁𝑻 = 𝜆 𝑒𝑓𝑓 ∇ 2 𝑇 + 𝐻 𝑀 . (2)
In equation (2), effective quantities (subscript "eff") indicate the properties of water-saturated rocks, and are obtained by averaging the properties of the fluid and of the dry matrix using porosity  and 1- as weights. C and  denote specific heat and thermal conductivity, respectively. HM is a magmatic heat source term (see Section 2.2).

The thermodynamic properties of hydrothermal fluids are approximated as those of a single-phase 3.2 wt% NaCl-H2O solution (Anderko and Pitzer, 1983;[START_REF] Pitzer | Thermodynamic properties of aqueous sodium chloride solutions[END_REF] with the temperature-dependent viscosity of pure water [START_REF] Holzbecher | Modeling Density-Driven Flow in Porous Media[END_REF]. Our model domain is 10-km wide, 30-km deep (Fig. 2), and uniformly gridded with a cell size of 100 m. The top boundary (seafloor, at a constant pressure of 20 MPa) is open to fluid inflow at a temperature of 0˚C, and can discharge hotter fluids (e.g., [START_REF] Fontaine | Two-dimensional numerical models of open-top hydrothermal convection at high Rayleigh and Nusselt numbers: Implications for mid-ocean ridge hydrothermal circulation[END_REF]. The sides and bottom of the domain are thermally insulated and closed to fluid flow such that the only heat advected out of the system at steady state is provided by repeated magma intrusions as described in section 2.2.

At t = 0, temperature is set to increase linearly from 0˚C to 1200˚C (basaltic liquidus) between the seafloor and 5 km depth, and remains at 1200˚C below (Fig. 2c).

Magmatic input

A magmatic heat source is implemented in the domain as a hAML = 500 m thick and wAML = 4 km wide sill instantaneously emplaced every time interval M in the center of the domain, immediately below the 1000˚C isotherm (basaltic solidus delineating the roof of the lens, Figs. 2, 3). Emplacement of sills at lower temperatures will be considered in section 5.3. Temperature inside the melt lens at the time of emplacement is set to the basaltic liquidus temperature of 1200˚C, which supplies specific heat to the system. Latent heat of crystallization is also accounted for during melt lens cooling following the parameterization of [START_REF] Fontaine | Hydrothermally-induced melt lens cooling and segmentation along the axis of fast-and intermediate-spreading centers[END_REF].

This amounts to expressing the HM term in equation (2) as:

𝐻 𝑀 = -𝐿𝜌 𝑀 𝑇 𝑀 -𝑇 𝑆 𝜕𝑇 𝜕𝑡 ⌈𝑥 𝐿 ⌉ , (3)
where L and M are the latent heat and density of magma, TM and TS the liquidus and solidus temperatures, and xL is the local melt fraction. The ⌈𝑥 𝐿 ⌉ term in equation (3) ensures that latent heat is accounted for only where melt is present, as xL is tracked and updated in the model as a linear function of temperature equal to 0 at T = TS, and 1 at T = TL. To avoid unrealistic temperature decrease with increasing depth below the melt lens, a narrow, 1-element wide dike accompanies each sill and connects it vertically to the base of the domain. The dike has the same properties as the melt lens, except for its latent and specific heat, which are scaled by a geometric factor to represent a 1-m wide magma fracture as a 1-element (100 m) wide column.

The periodicity of magmatic intrusions M is the main parameter that we vary in order to represent different magma budgets. While there is no restriction on the choice of M in our numerical model, we note that magma emplacement events cannot be so frequent that the corresponding volume flux of magma exceeds the accretion rate of new brittle lithosphere, which depends on the spreading rate and sill dimensions. This limit will be further discussed (and quantified) in Section 4. In the meantime, we explore a wide range of M values between 5 and 100 kyrs.

Evolving permeability field

We adopt the simplified model of [START_REF] Olive | Smoke without fire: How long can thermal cracking sustain hydrothermal circulation in the absence of magmatic heat?[END_REF] where permeability is set to a uniform value k in a shallow permeable zone overlying an area of negligible permeability. The permeable zone represents fractured lithosphere where fluids can easily circulate, and is delineated by a cracking front that can migrate over time. As postulated by deMartin et al.

(2004), the position of the cracking front corresponds to the temperature (TH) and pressure (or depth zH) at which tensile micro-cracks can stably propagate from defects present between olivine grains, in response to cooling stresses. This idealized model, detailed in Appendix A and illustrated in Fig. 2b, assumes square grains of size L with flaws of size a undergoing cooling below a visco-elastic transition temperature Tve. Thermo-elastic stresses of magnitude proportional to (Tve -TH) build up as the grains cool, which drives the opening of fractures.

Increasing confining pressure (P0 + gzH, where P0 denotes pressure at the seafloor and  is the density of the lithosphere) however counters the opening of cracks. Stable propagation is possible when the stress intensity factor at the flaws' tips, which balances these two opposing forces, reaches the critical fracture toughness KIC. This can be expressed as

𝐴 𝑇 (𝑇 𝑣𝑒 -𝑇 𝐻 ) -𝐴 𝑃 (𝑃 0 + 𝜌𝑔𝑧 𝐻 ) = 𝐾 𝐼𝐶 , (4)
(Fig. 2b) where AT and AP are coefficients that depend on grain and flaw size. Equation ( 4) is used to update the position of the cracking front at every time step of our simulations. This implies that as soon as stable cracking can occur at the grain scale, a network of interconnected fissures permeable to fluid flow instantaneously develops. This is clearly an oversimplification, which could be remedied by introducing a characteristic time lag between the onset of thermal cracking and the increase in permeability. We however choose to focus on the thermal steady state that develops over time scales likely greater than this time lag, thereby justifying our assumption of instantaneous permeability development.

Parameter values are summarized in Table S1. We assume constant grain and flaw sizes of 7.6 and 0.04 mm, respectively, and discuss the effect of these parameters in Section 4.3. This corresponds to AT = 3734 Pa.m 1/2 .K -1 and AP = 0.0112 m 1/2 (See Appendix A; and section 2.2. of [START_REF] Olive | Smoke without fire: How long can thermal cracking sustain hydrothermal circulation in the absence of magmatic heat?[END_REF]. We consider permeabilities of 510 -16 m 2 and 510 -15 m 2 in the shallow fractured region. The above values constitute reasonable, spatially-averaged estimates for poorly constrained parameters, and the greater value matches that used by [START_REF] Theissen-Krah | Coupled mechanical and hydrothermal modeling of crustal accretion at intermediate to fast spreading ridges[END_REF] to reproduce the equilibrium depth of melt lenses at magmatically robust slow-to fast-spreading MORs. While we do not test a wide range of values in the numerical simulations, we assess their influence on model behavior using an analytical approach in Section 4. Assuming a homogeneous and isotropic permeabity in the shallow domain is clearly an oversimplification [START_REF] Fisher | Permeability within basaltic oceanic crust[END_REF] whose implications have been previously assessed (e.g., [START_REF] Rosenberg | The relationship between flow and permeability field in seafloor hydrothermal systems[END_REF].

One particular consequence of ignoring lateral gradients in permeability is that our simulations cannot focus hydrothermal discharge within zones of high permeability potentially sealed off from lower-permeability recharge zones (e.g., [START_REF] Coumou | The structure and dynamics of mid-ocean ridge hydrothermal systems[END_REF][START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF]. This is likely to underestimate venting temperatures in our simulations. Finally, the 2-D nature of our models, while precluding the formation of realistic pipe-like upflow zones [START_REF] Coumou | The structure and dynamics of mid-ocean ridge hydrothermal systems[END_REF], capture the essential physics of heat transfer through porous convection that set the thermal regime of oceanic lithosphere [START_REF] Lowell | Hydrothermal processes at mid-ocean ridges: Results from scale analysis and single-pass models[END_REF].

Numerical results

We ran a first group of 10 numerical simulations with varying permeability and intrusion frequency (M = 5, 25, 50, 75 and 100 kyrs), using the thermal cracking parameters specified previously. All simulation inputs and results are summarized in Table S2. Every time a sill is intruded, it partially cools and crystallizes until the next episode of magma emplacement. The periodic release of specific and latent heat fuels hydrothermal convection, which typically consists of two plumes ascending along the domain sides, and a broad downwelling area in the middle of the permeable domain (Fig. 3a,d). Hydrothermal cooling progressively depresses the 1000˚C isotherm, and deepens the cracking front. The base of the permeable zone reaches depths of 5-7 km within a few Myrs. The area between the cracking front and the sill acts as a conductive boundary layer and thickens until the depth of subsequent sill intrusions stabilizes between 6 and 25 km, also within ~5 Myrs (Fig. 3b,e). Each simulation is carried out over 40

Myrs to make sure the thermal regime reached after a few Myrs is stationary.

The stationary thermal regime is characterized by a thicker permeable zone and a deeper melt lens, as well as colder seafloor venting temperatures when the recurrence time of intrusions is increased (Figs. 3,4). We find that equilibrium venting temperatures and depth to the cracking front are most sensitive to intrusion periodicities shorter than ~30 kyrs. Melt lens depth, however, increases linearly with M over the entire range considered here (Fig. 4b). Permeability does not strongly affect the steady-state depth of the cracking front or melt lens, but greatly influences venting temperatures. For example, while intrusions every 5 kyrs enable venting at 272˚C with a permeability of 510 -16 m 2 , this temperature is reduced to 112˚C when permeability is increased by a factor of 10.

To translate our predictions of stationary thermal regime into rheological properties, we use the laterally-averaged temperature profiles (Fig. 3b,e) to construct a yield strength envelope (Fig. 3c,d) assuming a background strain rate of 10 -15 , 10 -14 , and 10 -13 s -1 (Brace and Kohlstedt, 1980). These strain rates correspond to a full extension rate of 1 cm/yr distributed over a characteristic distance of 300, 30 (the characteristic extent of the deformation zone at the SWIR, [START_REF] Sauter | Continuous exhumation of mantle-derived rocks at the Southwest Indian Ridge for 11 million years[END_REF], and 3 km, respectively. Our envelopes also assume a friction coefficient of 0.6, 10 MPa of lithosphere cohesion and hydrostatic fluid pressure (fluid pressure ratio = 0.33) in the brittle domain, as well as wet olivine dislocation creep properties for the ductile domain [START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF]. In the case shown in Fig. 3a, which corresponds to a permeability of 510 - 16 m 2 and an intrusion periodicity of 50 kyrs, the brittle-ductile transition (for a reference strain rate of 10 -14 s -1 ) stabilizes at ~11 km (Fig. 3c). Doubling the intrusion periodicity lowers it to ~15 km (Fig. 3f). More generally, the depth to the brittle-ductile transition, which plausibly marks the base of the seismogenic layer, increases linearly with M (Fig. 4b).

In a second group of simulations we ascribe a permeabily of 510 -16 m 2 to the convecting region, but allow hydrostatic fluid pressure to relieve a portion of the confining stress that acts on the cracking front. This amounts to altering the criterion for cracking (equation 4) as:

𝐴 𝑇 (𝑇 𝑣𝑒 -𝑇 𝐻 ) -𝐴 𝑃 ((𝜌 -𝜌 𝑓 )𝑔𝑧 𝐻 ) = 𝐾 𝐼𝐶 , (5)
In equation ( 6), M is the density of the magma, and cPM, LM and TS are specific heat, latent heat of crystallization and solidus temperature, respectively. The above expression accounts for specific and latent heat release and assumes that melt fraction scales linearly with temperature between TS and TM. In the following we will assume TC = 1100˚C (unless otherwise specified), which amounts to 50% of crystallization following each intrusion event. Repeated intrusions every M account for a net vertical heat flux (in W/m 2 )

𝑞 𝑀 = Δ𝑄 𝑀 𝑤 𝐴𝑀𝐿 𝜏 𝑀 . (7)
Assuming steady-state heat transfers, this magmatic input must be matched by the hydrothermal heat flow qH. On the basis of Nusselt-Rayleigh number scalings for porous convection [START_REF] Lowell | Hydrothermal processes at mid-ocean ridges: Results from scale analysis and single-pass models[END_REF], [START_REF] Olive | Smoke without fire: How long can thermal cracking sustain hydrothermal circulation in the absence of magmatic heat?[END_REF] expressed this heat flow as:

𝑞 𝐻 = kΩ 𝑓 (𝑇 𝐻 -𝑇 0 ) 2 , ( 8
)
where T0 is the temperature of the deep ocean (0˚C) and f is a parameter combining several properties of rock and hydrothermal fluid:

Ω 𝑓 = 𝛼 𝑓 𝜆 𝑒𝑓𝑓 𝜁𝑔 𝜅𝜐 𝑓 , (9)
In equation ( 9),  denotes the thermal diffusivity of fluid-saturated rocks, f and f the thermal expansion and kinematic viscosity of hydrothermal fluids.  is the proportionality coefficient between Rayleigh and Nusselt numbers (here set to 0.2: [START_REF] Olive | Smoke without fire: How long can thermal cracking sustain hydrothermal circulation in the absence of magmatic heat?[END_REF]. Equation ( 8)

can thus be used to infer the temperature at the cracking front as a function of the magmatic heat input:

T 𝐻 = T 0 + √ 𝑞 𝑀 𝑘Ω 𝑓 . ( 10
)
This allows us to determine the depth to the cracking front using equation (4):

z 𝐻 = 1 𝜌𝑔 ((𝐴 𝑇 (𝑇 𝑣𝑒 -𝑇 𝐻 ) -𝐾 𝐼𝐶 )/𝐴 𝑃 -𝑃 0 ). (11)
Since the area between the cracking front and the melt lens acts as a conductive boundary layer transferring a heat flow qM, one can add its thickness to zH to get the equilibrium depth of sill intrusion:

z 𝑀 = 𝑧 𝐻 + 𝜆 𝑇 𝑀 -𝑇 𝐻 𝑞 𝑀 . ( 12
)
Finally, if we approximate the brittle-ductile transition as an isotherm TBDT embedded within the conductive boundary layer, its depth is readily obtained as

z 𝐵𝐷𝑇 = 𝑧 𝐻 + 𝜆 𝑇 𝐵𝐷𝑇 -𝑇 𝐻 𝑞 𝑀 . (13)
4.2 Sensitivity of thermal structure to magmatic input and permeability

The combination of equations (6-7) and (10-13) fully characterizes the thermal regime of the ridge axis as a function of permeability and intrusion frequency. Analytical results calculated with the same input parameters as the first group of numerical simulations are plotted in Figs. 4.

The temperature of the cracking front represents a reasonable estimate of simulated venting temperature (Fig. 4a) as vertical gradients in temperature are generally small within fluid upwellings in our numerical models (e.g., Fig. 3d). Our analytical model captures the observed trend of decreasing venting temperatures with increasing intrusion periodicity, as well as with increasing permeability (Fig. 4a). The latter is consistent with previous analytical modeling by [START_REF] Driesner | The interplay of permeability and fluid properties as a first order control of heat transport, venting temperatures and venting salinities at mid-ocean ridge hydrothermal systems[END_REF] and [START_REF] Barreyre | Depth-Dependent Permeability and Heat Output at Basalt-Hosted Hydrothermal Systems Across Mid-Ocean Ridge Spreading Rates[END_REF]. Differences between the numerical simulations and our analytical predictions likely stem from complex thermodynamic properties of the hydrothermal fluid strongly modulating the temperature of upwelling plumes (Jupp and Schulz, 2004). Our predictions of cracking front depth, brittle-ductile transition and melt lens depth (equations 12 and 13, Fig. 4b) are much more accurate, as the cracking rules for the numerical and analytical models are identical, and heat transfer between the permeable zone and the sill is dominantly conductive. This explains why the depth to the brittle-ductile increases linearly with M (at ~70% of the melt lens depth, Fig. 4b). The lack of sensitivity of cracking front depth to intrusion periodicities > 25 kyrs reflects the fact that hydrothermal convection has reached its maximum possible depth of ~6 km given this set of cracking parameters. Increasing permeability from 510 -16 m 2 to 510 -15 m 2 increases the depth of fluid percolation and sill emplacement by ~1 km, primarily reflecting the colder temperature of the cracking front in equations ( 11) and ( 12).

When intrusions are very frequent (< 10 kyrs), fluid convection stabilizes at a much shallower depth (<< 6 km) due to the strong magmatic input qM. We recall, however, that intrusions cannot be so frequent that they account for a crustal flux that exceeds the volume flux of newly created brittle lithosphere. This condition can be expressed as:

𝜏 𝑀 ≥ 𝑤 𝐴𝑀𝐿 ℎ 𝐴𝑀𝐿 𝑧 𝐵𝐷𝑇 𝑈 , (14)
where U denotes the full spreading rate. The corresponding range of implausible intrusion periodicities is shown in gray in Figs. 4 and5. For a spreading rate of 1.5 cm/yr, the minimum M is typically < 20 kyrs.

We next compare our analytical model to the second group of numerical simulations where hydrostatic fluid pressure relieves a portion of the confining stress at the cracking front, and fluids thus reach depths as large as ~10 km (blue curves and symbols in Fig. 5). Again, we find a good agreement between both types of models. This supports the notion that the heat flow that can be transferred through hydrothermal convection is independent of the thickness of the convective system (qH does not depend on zH in equation 8), which stems from the linear scaling between Nusselt and (high) Rayleigh numbers in porous convection [START_REF] Lowell | Hydrothermal processes at mid-ocean ridges: Results from scale analysis and single-pass models[END_REF]. This implies that if permeability was able to develop at even greater depths, our scalings would still provide an accurate description of the depth of the melt lens and brittle-ductile transition relative to the fluid percolation depth.

So far we have only considered 4-km wide, 500-m thick intrusions. However, the characteristic size of sills influences their heat input and competes with M as a control on the equilibrium thermal structure. Specifically, inspection of equations 2 and 3 shows that the heat balance is entirely controlled by the ratio of hAML / M with no influence of wAML. This is because magmatic heat is supplied over a horizontal length wAML and taken up by hydrothermal circulation over a similar length scale. wAML however modulates the net flux of new igneous material, and therefore controls the maximum intrusion frequency (equation 10). Fig. 4c illustrates the tradeoff between hAML and M. As an example, if the repeated intrusion of 500-m thick sills every 45 kyrs is required to maintain the brittle-ductile transition at 10 km, this recurrence time drops to ~10 kyrs for 100-m thick sills. We note that this rapid rate of intrusion does not exceed the rate of lithosphere formation, as the upper bound on M is itself a function of wAML and hAML.

Finally, we investigated the sensitivity of our models to the choice of thermal cracking parameters, by considering a plausible range of grain sizes (1 mm -1 cm) and flaw sizes (Fig. 6). We find that micromechanical parameters can shift the equilibrium depths of percolation, brittle-ductile transition and intrusion by ~4 km, but do not otherwise impact the effect of intrusion frequency and size. We generally find that deeper fluid percolation is favored by bigger grains populated with smaller defects. This said, our percolation model is only based on microfracturing induced by thermal stresses. Considering other drivers of tensile cracking such as tectonic stresses or reaction driven volume changes would likely deepen the reach of hydrothermal circulation.

Application to ultraslow-spreading ridges

Mid-Cayman Spreading Center

Our models provide a framework to interpret the deep seismicity recorded at ultraslowspreading ridges in terms of an intermittent supply of magma. A setting such as the MCSC for example features earthquakes down to ~10 km below seafloor (Figs. 1a; 3c) at the foot of the Mt.

Dent detachment fault [START_REF] Grevemeyer | Constraining the maximum depth of brittle deformation at slowand ultraslow-spreading ridges using microseismicity[END_REF]. This plausibly reflects the depth of the 650˚C isotherm on-axis, which is well explained by a model where a magmatic sill is emplaced every 50 kyrs at ~14 km depth (Fig. 3a,b). This depth is consistent with major element composition of MORB basalts from the MCSC indicating a mean crystallization pressure of ~5 kbars [START_REF] Wanless | Spreading rate-dependent variations in crystallization along the global mid-ocean ridge system[END_REF]. Repeated emplacement of 4-km wide, 500-m thick sills every 50 kyrs amounts to the same magma flux as accreting a 2.7-km thick magmatic crust at a full spreading rate of 15 mm/yr. This is in line with the scattered crustal structure inferred at the axis near Mt. Dent, which likely features a basaltic cover thinner than 3 km [START_REF] Harding | Magmatic-tectonic conditions for hydrothermal venting on an ultraslow-spread oceanic core complex[END_REF]. The characteristic size of sills at ultraslow ridges is however not well known and is likely to vary in space and time, although available geological constraints do point to sizes in the range of those investigated here [START_REF] Dick | Dynamic Accretion Beneath a Slow-Spreading Ridge Segment: IODP Hole 1473A and the Atlantis Bank Oceanic Core Complex[END_REF]. With smaller, 1-km by 100-m sills, each intrusion carries less specific and latent heat, and emplacement events must occur more frequently for the thermal structure to reach the same equilibrium state as in simulations using larger sills. As an example, equilibrating the brittle-ductile transition at ~10 km with 1-km by 100-m sills would require an intrusion periodicity of ~20 kyrs (Figs. 4c,5), and would result in a very tenuous crustal section (~0.3 km equivalent crustal thickness).

Another important aspect of our model is that it predicts hydrothermal convection to be confined within the uppermost ~6 km below seafloor, although depths as large as ~10 km can be envisioned if fluid pressure alleviates some of the confining pressure acting on the cracking front (Fig. 5), or if the lithosphere is made up of large grains with small defects (Fig. 6). While no onaxis seafloor venting has been documented to date at the foot of Mt. Dent where the deepest seismicity occurs, pertinent observations have been made beneath the Von Damm hydrothermal field located on the Mt. Dent detachment footwall, 15 km west of the middle of the axial valley (Fig. 1a). First, [START_REF] Harding | Magmatic-tectonic conditions for hydrothermal venting on an ultraslow-spread oceanic core complex[END_REF] identified a low P-wave velocity anomaly 5 km beneath Von Damm. They proposed that this anomaly could be the manifestation of a fluid-filled cracking front marking the base of hydrothermal convection. In our model --which only considers thermal cracking of idealized olivine aggregates [START_REF] Demartin | Experimental constraints on thermal cracking of peridotite at oceanic spreading centers[END_REF][START_REF] Olive | Smoke without fire: How long can thermal cracking sustain hydrothermal circulation in the absence of magmatic heat?[END_REF]) -fluid pathways can be created at this depth through fracturing of grains that are larger than ~5 mm and contain small intra-granular flaws (Fig. 5). The composition of the Von Damm hydrothermal fluids also suggest that pervasive fracturing may tap heat from a transient magmatic source below the site. Specifically, it was recently proposed that abiotic methane venting at Von Damm could originate from the reduction of CO2 released by magmatic intrusions, subsequently trapped in olivine fractures that formed and healed in overlying rocks.

The resulting CH4 is then tapped by circulating fluids as a new phase of hydrothermal cracking follows an intrusion event [START_REF] Klein | Abiotic methane synthesis and serpentinization in olivine-hosted fluid inclusions[END_REF][START_REF] Grozeva | Chemical and isotopic analyses of hydrocarbon-bearing fluid inclusions in olivine-rich rocks[END_REF].

Off-axis observations at Von Damm suggest that fluid convection can reach depths of ~5 km and be sustained by repeated intrusions. In the framework of our models, the occurrence of micro-earthquakes down to ~10 km on-axis could be explained by a stable boundary layer transferring heat conductively between the region of melt emplacement at ~14 km and the base of the convective system. Such a thick layer however cannot transfer large fluxes of heat, unlike the ~100-m thick boundary layers that exist above stable melt lenses at faster-spreading ridges (e.g., [START_REF] Gillis | The roof of the axial magma chamber: A hornfelsic heat exchanger[END_REF]. This is reflected in the moderate-to-low venting temperatures predicted in our simulations (e.g., 115˚C, Fig. 3a), which in a real system could manifest as diffuse venting that would be difficult to detect. Our model thus plausibly accounts for the thermal state of the MCSC axis near Mt. Dent, but does not apply to the off-axis Von Damm hydrothermal field. Its heat output of ~500 MW and venting temperatures up to 215˚C [START_REF] Hodgkinson | Talc-dominated seafloor deposits reveal a new class of hydrothermal system[END_REF] would indeed require much more frequent magmatic intrusions in the footwall of the Mt. Dent core complex (e.g., M ~ 5 kyrs, Fig. 4). These intrusions would preferentially occur at depths between 5 and 10 km, and could also contribute to the low-velocity anomaly observed beneath Von Damm [START_REF] Harding | Magmatic-tectonic conditions for hydrothermal venting on an ultraslow-spread oceanic core complex[END_REF]. Such a configuration has been inferred based on seismic data beneath the Rainbow site which is hosted on the footwall of a detachment on the Mid-Atlantic Ridge [START_REF] Dunn | Three-Dimensional Seismic Structure of the Mid-Atlantic Ridge: An Investigation of Tectonic, Magmatic, and Hydrothermal Processes in the Rainbow Area[END_REF]. A testable prediction of this model at the MCSC would be the occurrence of shallow microseismicity within ~5 km below the Von Damm site.

Amagmatic corridor, 64˚40'E on the Southwest Indian Ridge

The section of the SWIR located between 64˚30' and 65˚E is a type locality for amagmatic seafloor spreading (Fig. 1b). Its "smooth" seafloor consists of serpentinized mantle exhumed by a succession of large-offset detachment faults of alternating polarity, overlain in places by a scattered basaltic carapace [START_REF] Cannat | Modes of seafloor generation at a melt-poor ultraslow-spreading ridge[END_REF][START_REF] Sauter | Continuous exhumation of mantle-derived rocks at the Southwest Indian Ridge for 11 million years[END_REF]Cannat et al., 2019a). OBS surveys conducted west of this section [START_REF] Schlindwein | Mid-ocean-ridge seismicity reveals extreme types of ocean lithosphere[END_REF] as well as preliminary results from a recent OBS deployment [START_REF] Chen | Microseismicity constraints on brittle lithosphere thickness at a nearly amagmatic section of the ultraslow spreading Southwest Indian Ridge[END_REF] suggest that microearthquakes occur down to at least 15 km below seafloor at the axis. The recent ROVSMOOTH cruise also identified a low-temperature hydrothermal site named "Old City" midway up the scarp of the active axis-bounding detachment fault (Cannat et al., 2019b).

In the framework of our models, the only way for the brittle-ductile transition to equilibrate at a stable depth of ~15 km is to emplace a 4-km by 500-m sill every 100 kyrs at ~22 km (Fig. 3d-f; Fig. 4), or a smaller, for example 1-km by 100-m sill every ~20 kyrs (Fig. 4c).

The latter configuration would only account for the emplacement of a ~200-m thick equivalent igneous crust at a full rate of 15 mm/yr. In reality this crustal section would be extremely scattered and largely reworked by exhumation on large offset detachment faults [START_REF] Sauter | Continuous exhumation of mantle-derived rocks at the Southwest Indian Ridge for 11 million years[END_REF]. It is therefore compatible with the scarcity of the extrusive cover at the seafloor and the prevalence of serpentinized mantle lithologies at 64˚E. Our models are also consistent with the occurrence of low-temperature hydrothermal venting at Old City (Cannat et al., 2019b). Because fluids cannot reach depths greater than ~6 km in our simulations, heat must be transferred through a very thick conductive boundary layer between the depth of repeated magmatic emplacement and the fluid convection system. This strongly limits the heat output of hydrothermal circulation, with venting temperatures that cannot exceed ~100˚C. While direct measurements of output temperatures at Old City are not yet available, the carbonate-brucite chimneys documented there are comparable to observations at the Lost City vent field, Mid-Atlantic Ridge (Früh-Green et al., 2003) and indicative of low-temperature, low-heat flow venting fueled by heat transferred through a thick boundary layer, with a possible contribution from heat released by alteration reactions.

Segment 28 of the Southwest Indian Ridge

Segment 28 of the SWIR located at 49.7˚E spreads orthogonally at a rate very close to that of the MCSC, and also features a prominent detachment fault and focused hydrothermal venting [START_REF] Zhao | Three-dimensional seismic structure of the Dragon Flag oceanic core complex at the ultraslow spreading Southwest Indian Ridge (49˚39'E)[END_REF]. Unlike Von Damm, the Longqi vent field is located at the toe of the Dragon Flag detachment, ~5 km from the neovolcanic axis, and discharges hydrothermal fluids as hot as 379˚C [START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF]. Seismicity data indicates a thicker brittle lithosphere or of 100-m thick sills every ~20 kyr (Fig. 6) at depths of 20-25 km. However, because hydrothermal circulation is confined within the uppermost ~6 km of the lithosphere, the boundary layer separating the convecting system from the heat source must be >15 km thick, implying that deep earthquakes occur in a region devoid of fluid convection. While this steadystate thermal structure is compatible with the weak hydrothermal output of Old City at 64˚E, it cannot explain the high-temperature and high-flux (~250 MW) discharge documented at Longqi [START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF]. Such a high output implies that a heat flow of order 100 W/m 2 is transferred conductively through a < 100-m thick boundary layer [START_REF] Lowell | Characteristics of magma-driven hydrothermal systems at oceanic spreading centers[END_REF]. [START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF] propose that this occurs at a depth of 16 km, which marks the base of the fluid convection system and the locus of sustained magmatic intrusions maintaining a temperature of ~650˚C (Fig. 1d).

This is considerably deeper than the 6 km (~10 km at most) permitted by our cracking model [START_REF] Demartin | Experimental constraints on thermal cracking of peridotite at oceanic spreading centers[END_REF][START_REF] Olive | Smoke without fire: How long can thermal cracking sustain hydrothermal circulation in the absence of magmatic heat?[END_REF], which cannot account for fluid percolation under high confining pressure. Pressure indeed strongly limits permeability past a few km into the oceanic crust, as evidenced by borehole measurements [START_REF] Fisher | Permeability within basaltic oceanic crust[END_REF][START_REF] Kuang | An integrated permeability-depth model for Earth's crust[END_REF]. As described in Section 1, we also lack clear evidence for widespread hydrothermal alteration of the lithosphere at depths of 10-15 km. Nonetheless, deep fracturing could be locally driven by additional sources of stress such as volume changes associated with alteration reactions (e.g., [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], or stresses induced by magmatic intrusions [START_REF] Wilcock | The role of magma injection in localizing black-smoker activity[END_REF], provided they can overcome the inhibiting effect of pressure on crack growth.

An alternative to steady magma emplacement fueling deep fluid circulation would be to invoke infrequent shallow intrusions that temporarily perturb the long-term equilibrium predicted by our numerical models. Episodic emplacement of melt pockets in rocks cooler than the solidus is suggested by both seismic imaging (e.g., [START_REF] Dunn | Three-Dimensional Seismic Structure of the Mid-Atlantic Ridge: An Investigation of Tectonic, Magmatic, and Hydrothermal Processes in the Rainbow Area[END_REF] and petrological studies (e.g., [START_REF] Paquet | Effect of melt/mantle interactions on MORB chemistry at the easternmost Southwest Indian Ridge (61°-67°E)[END_REF]. To assess the thermal impact of an infrequent shallow intrusion, we conducted an additional simulation using the same setup as shown in Fig. 3d-f (same cracking parameters, intrusion size, permeability of 510 -16 m 2 , M = 100 kyrs), this time allowing one every ten intrusions to occur at a depth of 5 km irrespective of temperature, instead of at the 1000˚C isotherm. Results are shown in Figs. 7 and8. We find that the emplacement of a shallow sill once per million year does not bring enough heat to alter the previously described stationary thermal structure in a significant way. The 600˚C and 1000˚C stabilize within the conductive boundary layer at ~15 and 22 km, respectively, as in the initial simulation (Fig. 3d-f), and are not affected by the transient shallow intrusions which occur within the hydrothermal convection domain (Fig. 8a). Flow patterns and temperature structure are however disrupted in the permeable zone for a few hundred thousand years following each shallow intrusion. Specifically, two new upwelling zones develop immediately above a shallow intrusion, and migrate towards the edges of the domain where they merge with the stable upwelling zones, bringing a burst of heat to the seafloor within 50 kyrs of the intrusion. The venting temperature increases from its steady-state value of 99˚C to ~140˚C, and remains at that value for ~10 kyrs. This time span could be extended, and this peak temperature increased if the shallow heat source was maintained (i.e., replenished) for a finite amount of time, as opposed to instantaneously emplaced. Greater permeabilities --on the order of 10 -14 -10 -13 m 2 --would also enable a greater output of heat during these transient events by allowing the formation of a thinner conductive boundary layer above the shallow sill [START_REF] Lowell | Hydrothermal processes at mid-ocean ridges: Results from scale analysis and single-pass models[END_REF]. Overall, our models suggest that a shallow, replenished intrusion could establish a vigorous, but ephemeral hydrothermal circulation system similar to those of faster spreading ridges that could explain high-output venting at Longqi [START_REF] Liu | Models of hydrothermal heat output from a convecting, crystallizing, replenished magma chamber beneath an oceanic spreading center[END_REF][START_REF] Lowell | Characteristics of magma-driven hydrothermal systems at oceanic spreading centers[END_REF], provided this intrusion occurred within the last tens of kyrs. Its effect on the temperature field would however be local and confined to the areas above the shallow sill, as it would not have had the time to conductively warm the underlying, thick steady boundary layer. This means that the shallow sill would cause a local maximum in the geotherm, which would still remain at ~650˚C at 15 km depth, enabling brittle deformation at great depth. This solution is not fully satisfactory at Longqi because the seismic velocity structure of the crust does not suggest a shallow temperature anomaly, or pockets of melt within the uppermost 13 km of the lithosphere [START_REF] Yu | Lithospheric Structure and Tectonic Processes Constrained by Microearthquake Activity at the Central Ultraslow-Spreading Southwest Indian Ridge[END_REF][START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF]. A single, recently intruded sill of width ~100 m may, however, be below the resolution of seismic imaging.

Steady vs. transient thermal regimes at ultraslow ridges

The thermal structure of ultraslow-spreading ridges is unique among MORs as it involves a thick brittle lithosphere which may or may not be entirely permeable to fluid flow, as well as a small, transient input of magmatic heat (Fig. 1). We have designed numerical and analytical models that explore the thermo-mechanical consequences of conductive and convective heat transfers through young ultraslow-spread lithosphere. These models suggest that if magma pools below its liquidus isotherm, and hydrothermal convection is confined within a fixed depth range, the base of the brittle lithosphere and the depth of magma emplacement will equilibrate over time to reflect the thickness of the conductive boundary layer beneath the fluid percolation front.

Less frequent magmatic intrusions provide a lesser heat input, which can be transferred through a thicker boundary layer, making the lithosphere thicker and the magma fractionation zone deeper.

Because hydrothermal heat extraction weakly depends on the thickness of the convective region, the sole effect of creating deeper permeable pathways for fluids is to shift the conductive boundary layer toward greater depths without affecting its thickness.

In this conceptual framework, 10-15 km-deep earthquakes below the quasi-amagmatic seafloor of the SWIR at 64ºE occur in the conductive boundary layer that transfers heat from infrequent, small intrusions at ~22 km to a weak, ~6-km thick hydrothermal convective system fueling the low-temperature Old City vent field. The shallower limit of micro-earthquakes at a more magmatic ultraslow ridge such as the MCSC could simply reflect more frequent intrusions that provide more heat through a thinner boundary layer below the convective zone. This would imply that the SWIR and 64ºE and the MCSC are both in a thermal steady state controlled by the frequency of their magmatic input. Although we did not systematically explore the transient phase that precedes the establishment of this stationary regime, our simulations suggest that isotherm depths and venting temperatures equilibrate within ~1-2 Myr (Fig. 7a). [START_REF] Harding | Magmatic-tectonic conditions for hydrothermal venting on an ultraslow-spread oceanic core complex[END_REF] proposed that a pulse of magmatism occurred at ~2 Ma at the MCSC, manifesting today as a cooled gabbro body exhumed in the footwall of the Mt. Dent detachment. If the supply of magma waned and remained steady after this episode, the present day thermal structure could plausibly reflect a stationary state. At 64ºE on the SWIR, the present-day mode of quasiamagmatic accretion involving flip-flopping detachments has been active for ~11 Myrs [START_REF] Sauter | Continuous exhumation of mantle-derived rocks at the Southwest Indian Ridge for 11 million years[END_REF]. This would have provided ample time for the thermal structure to equilibrate with a steady, subdued magma supply. It is also possible that the time scales of equilibration in our models reflect the average permeability ascribed to the convective domain. [START_REF] Olive | Smoke without fire: How long can thermal cracking sustain hydrothermal circulation in the absence of magmatic heat?[END_REF] indeed found that the cracking front migrates faster for greater permeabilities, in the absence of magmatic heat input. Future work should focus on elucidating how rapidly a ridge can switch from one thermal steady state to another in response to a rapid change in magma supply.

On the other hand, the fact that a fairly magmatic setting such as Segment 28 of the SWIR also hosts 15-km deep earthquakes and vigourous hydrothermal venting challenges those assumptions. One explanation (Fig. 1d) could be that fluid circulation takes place down to ~15 km and convectively cools the lithosphere to such great depths [START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF]. Magma would then get emplaced close to the percolation front, forming a thin boundary layer that transfers large fluxes of heat [START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF]. However, the fact that P waves reach unaltered mantle velocities of 8 km/s at 10 km below the axial volcanic ridge at segment 28 [START_REF] Yu | Lithospheric Structure and Tectonic Processes Constrained by Microearthquake Activity at the Central Ultraslow-Spreading Southwest Indian Ridge[END_REF] appears incompatible with large amounts of igneous crust being emplaced at even greater depths over longer time scales. Our own model, in which fluid circulation is confined to shallower depths, also cannot account for these observations under the assumption of a stationary state, because it would invoke a thermal boundary layer much too thick to allow the heat output measured at Longqi. Instead, we propose that present-day active venting reflects a recent shallow intrusion, possibly too small to be resolved by seismic tomography methods, which occurred within the permeable zone and enabled a thin transient boundary layer and the extraction of a burst of magmatic heat (Fig. 8). In this interpretation, high-temperature venting in thick brittle lithosphere would be an inherently transient process.

Conclusions

We have designed a new model for the thermal structure of ultraslow ridge axes in which episodic magma intrusions provide heat to a hydrothermal convection system confined to depths < 6 km. Our model successfully explains the depth to the brittle-ductile transition as well as the depth of MORB crystallization at the MCSC by invoking repeated emplacement of sills every 20-50 kyrs. It is also consistent with the seismically-inferred depth range of fluid-filled cracks.

Increasing this intrusion periodicity and / or reducing the intrusion size would straightforwardly explain the thicker (~15 km) seismogenic layer as well as the low-temperature venting observed at the amagmatic section of the SWIR at 64ºE. However, unless mechanisms we have not accounted for enable deeper fluid percolation, our theoretical framework requires that hightemperature venting at ultraslow ridges be fueled by episodic magma intrusions into the longterm convective region, which could fuel high-output circulation for a few tens of kyrs.

Figures and captions

  where f = 900 kg/m 3 is the average density of the fluid. Coefficients AT and AP are left unchanged. Consequently, the cracking front must migrate deeper to experience the same effective pressure PH as in the first group of simulations (equation A4). This enables a deeper convective region, down to ~10 km. We run six simulations with this modified cracking criterion, using M = 5, 15, 25, 50, 75 and 100 kyrs. The results are plotted in Fig.5alongside the results of previous simulations with shallower convection ( 6 km) and k = 510 -16 m 2 . The overall model behavior is similar to that of the first set of runs: the effect of allowing deeper fluid convection is simply to shift the steady-state brittle-ductile transition and depth of magma emplacement consistently by ~4 km. 4 Key controls on ultraslow ridge thermal structure 4.1 Analytical heat balance at a magma-poor ridge axis To gain further insight into the behavior of our numerical simulations, we develop a simple analytical model of the heat balance between hydrothermal cooling and transient magmatic input. The amount of heat per unit ridge length QM supplied by a sill of width wAML and thickness hAML cooling from its emplacement (liquidus) temperature TM down to a temperature TC writes: Δ𝑄 𝑀 = 𝜌 𝑀 𝑤 𝐴𝑀𝐿 ℎ 𝐴𝑀𝐿 (𝑇 𝑀 -𝑇 𝐶 ) (

  compared to the MCSC:[START_REF] Yu | Lithospheric Structure and Tectonic Processes Constrained by Microearthquake Activity at the Central Ultraslow-Spreading Southwest Indian Ridge[END_REF] and[START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF] reported earthquakes and microearthquakes down to ~15 km below seafloor throughout Segment 28 and within the detachment footwall. As detailed above, our model requires the repeated emplacement of 500-m thick sills every ~100 kyr to account for such a deep brittle-ductile transition (Fig.3d-f; Fig.4),

1 .

 1 Figure 1. Schematic representation of 3 ultraslow ridge sections with axial seismicity within or exceeding 10 km below seafloor, with inferred heat source and geometry of heat transfer system. DF: Detachment fault. CBL: Conductive Boundary Layer. HT/LT: High/Low-Temperature. a.Mid-Cayman Spreading Center (MCSC) near Mt. Dent[START_REF] Harding | Magmatic-tectonic conditions for hydrothermal venting on an ultraslow-spread oceanic core complex[END_REF]. b. Amagmatic corridor at 64˚40'E on the Southwest Indian Ridge (SWIR)[START_REF] Sauter | Mantle exhumation at the Southwest Indian Ridge; preliminary results of the SMOOTHSEAFLOOR cruise[END_REF][START_REF] Chen | Microseismicity constraints on brittle lithosphere thickness at a nearly amagmatic section of the ultraslow spreading Southwest Indian Ridge[END_REF]. c. Segment 28 of the SWIR near 49˚E. d. Same location, with model for heat transfer from[START_REF] Tao | Deep high-temperature hydrothermal circulation in a detachment faulting system on the ultra-slow spreading ridge[END_REF].

Figure 2 .

 2 Figure 2. Numerical model setup: (a) domain geometry and boundary conditions; (b) details of the cracking front model; (c) Initial temperature profile.

Figure 3 .

 3 Figure 3. Thermo-mechanical state of an ultraslow ridge axis from numerical simulations that have reached stationary state with a shallow permeability of 510 -16 m 2 , subjected to a 400-m  1-km magmatic intrusion every 50 kyrs (a-c) and 100 kyrs (d-f). (a, d) Temperature snapshot showing the axial melt lens (AML) and the base of the permeable zone (dashed white line). (b, d) Laterally-averaged temperature profile. CBL: conductive boundary layer. (c, f) Corresponding yield strength envelope (differential stress, tensile) calculated with the wet olivine rheology of Karato et al. (1986) and a strain rate of 10 -14 s -1 . BDT: Brittle-ductile transition. Green bars are histograms of micro-earthquake depths at the Mid-Cayman Spreading Center (MCSC) and the Southwest Indian Ridge (SWIR), 49˚E, from[START_REF] Yu | Lithospheric Structure and Tectonic Processes Constrained by Microearthquake Activity at the Central Ultraslow-Spreading Southwest Indian Ridge[END_REF].

Figure 4 .

 4 Figure 4. Key model outputs plotted against the periodicity of 500-m thick and 4-km wide magmatic intrusions. Red and blue colors correspond to two different permeabilities. (a) Diamonds: maximum venting temperatures predicted by the numerical simulations in stationary state. Error bars reflect temporal fluctuations when present. Lines: temperature at the bottom of the hydrothermal system predicted from the analytical model. (b) Thickness of the permeable zone predicted by the numerical simulations in stationary state (crosses) and the analytical model (dashed lines, see Section 4). Depth to the brittle-ductile transition in the numerical models(squares, estimated by constructing yield strength envelopes with a strain rate of 10 -14 s -1 , as shown in Fig.3). Depth to the axial melt lens from the numerical simulations (circles) and the analytical model (lines). Horizontal bands indicate the maximum depth of seismicity (proxy for the brittle-ductile transition: BDT) at the Mid-Cayman Spreading Center (MCSC, orange) and at the Southwest Indian Ridge, 49˚E and 64˚E (SWIR, green). (c) Trade-off between thickness and periodicity of magmatic intrusions: red lines show the combination of hAML and M that predict a

Figure 5 .

 5 Figure 5. Effect of deeper fluid percolation on key model outputs. Red and blue colors correspond to two different hypotheses for pressure at the cracking front (full lithostatic vs. partially relieved by fluid pressure, i.e., equations 4 and 5, respectively), assuming a permeability of 510 -16 m 2 in the convective region. (a) Diamonds: maximum venting temperatures predicted by the numerical simulations in stationary state. Line: temperature at the bottom of the hydrothermal system predicted from the analytical model. (b) Thickness of the permeable zone predicted by the numerical simulations in stationary state (crosses) and the analytical model (dashed lines, see Section 4). Depth to the axial melt lens from the numerical simulations (circles) and the analytical model (lines). Other symbols follow conventions of Fig. 4. Here the analytical model assumes each intrusion cools to TC = 1080ºC before the next one, as this slightly improves the fit to the results from the second group of numerical simulations with deeper fluid percolation.

Figure 6 .

 6 Figure 6. Thermo-mechanical properties of an ultraslow ridge axis plotted against the periodicity of 100-m thick and 1-km wide magmatic intrusions. Analytical model assumes a permeability of 10 -16 m 2 and a range of micro-mechanical parameters: 1-mm grains with 0.1-mm and 0.01-mm flaws (green and red); 1-cm grains with 0.1-cm and 0.01-cm flaws (orange and blue). Horizontal bands indicate the maximum depth of seismicity (proxy for the brittle-ductile transition: BDT) at the Mid-Cayman Spreading Center (MCSC, orange) and at Segment 28 of the Southwest Indian Ridge, 49˚E (SWIR, green). (a) Venting temperature. (b) Depth of fluid percolation (long dashed lines), brittle-ductile transition (BDT, solid lines), and axial melt lens (AML, short dashed lines).

Figure 7 .

 7 Figure 7. a. Maximum venting temperature throughout the simulation shown in Fig. 3d-f (blue curve, intrusions steadily occurring every 100 kyrs at the 1000˚C isotherm). In red, venting temperatures from the same simulation setup, with 1 out of 10 intrusions occurring at 5 km depth. b. Zoom on the time span following a shallow intrusion. Timing of the snapshots from Fig. 8 are shown.

Figure 8 .

 8 Figure 8. Snapshots of the temperature structure from a model similar to that of Fig. 3d-f, where 1 out of 10 intrusions (i.e., once per Myrs) occurs at 5 km depth instead of at the much deeper 1000˚C isotherm. Timing of snapshots is indicated in Fig. 7.

Acknowledgments

This work was supported by China Scholarship Council Fellowship #201806010083 (Q.F.), and NSF grant OCE-1536943 (J.-A.O). The authors wish to thank Ingo Grevemeyer and Lars Ruepke for their helpful reviews, and Tim Crone and Javier Escartín for fruitful discussions. The source code used to produce our numerical model is available at https://figshare.com/s/3c506b1aeafd28dff5ee.

Appendix A: Cracking front model for the depth extent of hydrothermal circulation

We assume that the permeable hydrothermal convection zone is bounded by a cracking front that can migrate downward as the system cools [START_REF] Lister | On the penetration of water into hot rock[END_REF][START_REF] Olive | Smoke without fire: How long can thermal cracking sustain hydrothermal circulation in the absence of magmatic heat?[END_REF]. Specifically, in the framework of [START_REF] Fredrich | Micromechanics of thermally induced crack-ing in three crustal rocks[END_REF][START_REF] Fredrich | Micromechanics of thermally induced crack-ing in three crustal rocks[END_REF][START_REF] Demartin | Experimental constraints on thermal cracking of peridotite at oceanic spreading centers[END_REF], anisotropic thermal contraction of olivine grains drives the build-up of thermo-elastic stresses when rocks cool below temperature Tve. These stresses act orthogonally to grain boundaries (oriented along direction x, Fig. 2) and can be expressed as:

Where E = 197 GPa and = 0.25 are the Young's modulus and Poisson's ratio of the grains, = 3.1 10 -6 K -1 is the difference in thermal expansion coefficient between orthogonal directions in olivine grains, and 2L is the grain size. f(x,L) is a shape function defined as:

. (A2)

These stresses may induce the stable growth of tensile cracks from pre-existing flaws of size a located along grain boundaries. This occurs when the fracture toughness of the system KI reaches a critical value KIC. The former can be expressed as:

,

where P H denotes the confining pressure acting at the cracking front depth (z H ). Condition K I = K IC therefore writes:

, (A4)

with AT and AP defined as: ,

(A6)

We further assume that a network of connected cracks allowing fluid percolation forms instantly as KI reaches KIC. We therefore use equation (A4) as the criterion marking the base of the permeable zone in our simulations.