Retrieving timescales of oceanic crustal evolution at Oceanic Core Complexes: Insights from diffusion modelling of geochemical profiles in olivine

The building of oceanic crust at Oceanic Core Complexes (OCC) has been described as a complex process involv- ing multiple intrusions of magma over a protracted period of time. The migration of primitive magmas (i.e., Mid Ocean Ridge Basalts, MORBs) can lead to melt-rock interactions during reactive porous ﬂ ow processes through the lithosphere. The timescales of these reactive processes and the subsequent cooling of the modi ﬁ ed crystal matrixremain unconstrained.Diffusion modelling hasbeen widely used to retrieve timescales ofmagmaticpro- cesses. In this study, we use diffusion models to constrain (i) the minimum timescales of melt-rock interactions and (ii) the cooling ratesof the gabbroic sequence forming the oceanic crust atan OCC. We chose samples of the mostprimitiveolivine-richtroctolitesfromthegabbroicsequencesampledinIODPHoleU1309D(AtlantisMassif OCC, Mid-Atlantic Ridge 30°N). Olivine-rich troctolites were interpreted as marking local partial assimilation of mantle intervals into the oceanic crust, and thus allowed to understand the dynamics of mantle assimilation and the formation of slow-spreading oceanic crust at the Atlantis Massif OCC. Olivines in olivine-rich troctolites represent relicts of pre-existing mantle olivine, while clinopyroxenes and plagioclases are crystallized during reactive percolation. Olivinechemicalcompositionsshow thatolivine-rich troctolitesinherit chemical heteroge- neity fromthemantleprecursor. Flat geochemicalpro ﬁ les inolivine indicate complete chemicalre-equilibration of olivine crystals with the locally modi ﬁ ed percolating melt. Exception is made for most Ca pro ﬁ les that show lower Ca contents at the olivine rim compared to the relative crystal core, as the result of subsolidus cooling. Three-dimensional (3D) diffusion models at magmatic conditions ( T = 1210 – 1300 °C and P = 2 kbar) reveal that complete chemical re-equilibration of 4 mm-size mantle-derived olivine with percolating MORB-type melts can beattained within durationsofless than 300yr. The Ca-in-olivine geospeedometer revealthatcooling rates from ~1200 °C to ~1050 °C are constant downhole and on average 0.004 °C/yr; they are comparable with lower temperature cooling rates(850 °C – 250 °C) estimatedatthe AtlantisMassif OCC. The minimum timescales from 3D models point to rather fast re-equilibration of olivine. The downhole chemical heterogeneity inherited from the precursor mantle, coupled with the timescales of diffusive re-equilibration suggest that the partial as- similation of the upwelling mantle and its incorporation into the oceanic crust occurred in the time-frame of a singlemeltinput.Ourcooling data areconsistentwithbuilding oftheoceanic crustatOCCscontrolledbycontin- uous uplift, inturn governed bylong-lived detachment faults. Thelattercontributetotherapidcooling of theas-similated mantle intervals and the magma bodies. Diffusion models of geochemical pro ﬁ les in olivine from a single crustal section allow to reconstitute the early magmatic processes leading to mantle assimilation and early crystallization of gabbros, and the cooling history of the oceanic crust at OCC from magmatic conditions to hydrothermalism.


Introduction
Melts generated by adiabatic melting of the upwelling mantle beneath mid-ocean ridges are supplied to the magmatic plumbing system, where they undergo complex processes of magma crystallization, mixing and/or assimilation of a pre-existing crystal matrix prior to eruption on the seafloor (e.g., Brandl et al., 2016;Coogan and O'Hara, 2015;O'Hara, 1977;O'Neill and Jenner, 2012). These processes are recorded in the lower oceanic crust beneath oceanic spreading centers. Recent studies have demonstrated the importance of reactive porous flow and melt-rock interactions in the construction of the gabbroic lower crust, in particular along slow-spreading ridges (e.g., Dick et al., 2010;Drouin et al., 2009Drouin et al., , 2010Lissenberg and Dick, 2008;Lissenberg and MacLeod, 2016;Lissenberg et al., 2019;Sanfilippo et al., 2015) and in their ophiolite analogues (e.g., Basch et al., 2019aBasch et al., , 2019bBasch et al., , 2018Bédard, 1993Bédard, , 2001Bédard and Hébert, 1996;Rampone et al., 2020;Sanfilippo et al., 2014). These reactive processes are thought to occur in cooling yet still partially molten magmatic systems, possibly at the transition between lithospheric mantle and the magmatic lower oceanic crust (e.g., Collier and Kelemen, 2010;Ferrando et al., 2018;Rampone et al., 2020;Sanfilippo et al., 2015). To better understand the mechanisms driving the formation of slow-spreading oceanic crust, it is imperative to constrain the timescales of the processes involved in its evolution from early magmatic differentiation of melts (i.e., crystallization and dissolution-precipitation reactions) to cooling of the crustal sequence (i.e., crystal-dominated system).
Timescales of magmatic emplacement beneath slow-spreading ridges have been determined for crustal sections exposed along the Mid-Atlantic Ridge (MAR; e.g., Grimes et al., 2008), while duration of cooling have been concurrently estimated by other authors for the lower oceanic crust beneath both the MAR and the Southwest Indian Ridge (SWIR) (e.g., Grimes et al., 2011;John et al., 2004;Schoolmeesters et al., 2012). John et al. (2004) documented the thermal structure of oceanic lithosphere beneath the ultraslow-spreading SWIR. They applied thermochronometric techniques on minerals in felsic veins (crystallization at~850°C), thus lacking direct constraints on the earlier magma emplacement history at higher temperatures (from 1250°C to~850°C). A broad reconstruction of timescales of oceanic crustal evolution from high temperature magma emplacement to cooling of the gabbroic sequence has not yet been established for a single suite of samples. In this contribution, we use the geochemical composition of olivines from the widely investigated primitive olivine-rich troctolites recovered at Atlantis Massif (IODP Hole U1309D; Blackman et al., 2006) to provide timing of magmatic processes occurring at temperatures of about 1230°C and of cooling of the lower oceanic crust down to about 1000°C.
Olivine-rich troctolites from the oceanic crustal section at the Atlantis Massif Oceanic Core Complex (OCC; IODP Hole U1309D; Blackman et al., 2006) are one of the best studied olivine-rich lithologies (i.e., with >70% modal olivine) worldwide. They were interpreted to result from partial mantle assimilation during reactive porous flow of a MORB-type melt (Drouin et al., , 2009(Drouin et al., , 2010Ferrando et al., 2018;Suhr et al., 2008). Olivine-rich troctolites from Hole U1309D were selected for this study because they are unique samples from active spreading ridges that record the oceanic crustal evolution from early magmatic processes of melt-rock interaction and crystallization of a primitive MORB, to exhumation and cooling during OCC formation. We apply diffusion modelling to provide a comprehensive dataset of timescales related to these magmatic processes and to reconstruct the geodynamic evolution of the oceanic crust from the Atlantis Massif OCC.
Diffusion modelling of chemical species in rock-forming minerals is a tool widely used to quantify the timescales of magmatic processes such as magma assimilation (e.g., Bindeman et al., 2006;Costa and Dungan, 2005), magma mixing (e.g., Chamberlain et al., 2014;Lynn et al., 2017a), magma residence time before eruption (e.g., Lynn et al., 2018;Nakamura, 1995), and subsolidus cooling rates of magmatic bodies (e.g., Coogan et al., 2007Coogan et al., , 2002Faak et al., 2013;Faak and Gillis, 2016;Sun and Lissenberg, 2018). Many theoretical and experimental studies have investigated the crystallographic and chemical parameters controlling diffusion of major, minor and trace elements in mafic minerals, most particularly olivine, thus significantly improving the available diffusion coefficient database (e.g., Coogan et al., 2005;Petry et al., 2004;Spandler and O'Neill, 2010). Concurrently, progresses in analytical techniques now allow measurement in situ of a broad range of chemical elements down to very low concentrations (<<ppm) such as Rare Earth Elements (REE) in olivine (e.g., Basch et al., 2018;D'Errico et al., 2016;Drouin et al., 2009;Ferrando et al., 2018;Rampone et al., 2016;Sanfilippo et al., 2014). These developments allow the collection of chemical profiles for many elements that have different rates of diffusion, thus permitting detailed investigations on the timescales of complex processes.
In this study, we used geochemical profiles (Ferrando et al., 2018 and this study) in olivine because it is the only phase present from the mantle protolith to the reactive formation and cooling of olivine-rich troctolites. Geochemical profiles record two processes: (i) flat profiles of olivine forsterite content (Fo) and Ni, Mn, Co and Zn compositions indicate complete re-equilibration of olivine with the percolating melt during melt-rock interactions; and (ii) convex profiles of Ca, Y and Yb document subsolidus re-equilibration during cooling. However, timescales of the high-temperature (hereafter referred to as 'magmatic temperatures') reactive process and cooling remain to be estimated.
(i) We reproduce the measured homogeneous element profiles by numerically modelling their chemical diffusive re-equilibration and retrieve minimum timescale estimates of high temperature (> 1200°C) re-equilibration. To account for diffusion anisotropy in olivine (D [001] > D [010] ≈ D [100] for Fe\ \Mg, Ni and Mn, D [001] >> D [100] > D [010] for Co, while Ca, Y and Lu are nearly isotropic; Table 1 and references therein), and to consider the influence of spatial dimensions, crystal morphology, and sectioning, we used a three-dimensional (3D) numerical model (Jollands and Müntener, 2019;Lynn et al., 2017b;Shea et al., 2015). Simultaneous diffusive re-equilibration of eight diffusing species comprising major, minor and trace elements (Fe\ \Mg, Ni, Mn, Ca, Co, Zn, Y, Yb) in olivine were modelled. Our results provide the first constraints on the timing of melt-rock interactions largely described at oceanic spreading centers.
(ii) Subsolidus cooling rates have been previously investigated in the oceanic crust at the Atlantis Massif OCC, although most constraints document the latest stages of cooling (below 780°; Grimes et al., 2011;Schoolmeesters et al., 2012) and cooling rates in the temperature interval from magmatic conditions (~1200°C) down to 900-800°C remain unconstrained. In this study, calcium zoning in olivine is modelled using the Ca-in-olivine geospeedometer (Coogan et al., 2007(Coogan et al., , 2002 to determine~1200°C to~1000°C subsolidus cooling rates of Hole U1309D olivine-rich troctolites.
The results of this study characterize the timing of mantle incorporation into the oceanic crust, which are used to reconstruct the thermal history and subsolidus cooling of the gabbroic sequence formed at an OCC, from high magmatic temperatures down to~250°C.

The oceanic crustal sequence at Atlantis Massif
Atlantis Massif is a domal structure located at 30°N on the western flank of the slow-spreading MAR, at the intersection with the Atlantis transform fault (Fig. 1a). The corrugated core of the dome is a~2 Ma old OCC where oceanic crust is exposed via a long-lived and low-angle detachment fault (e.g., Blackman et al., 2002;Ildefonse et al., 2007).
On the central dome of Atlantis Massif, IODP Hole U1309D (Fig. 1a) penetrated 1415.5 m below seafloor (mbsf) through a complex gabbroic sequence (IODP Expedition 304/305; Blackman et al., 2011). The recovered section is highly heterogeneous and comprises~85% of gabbro, gabbronorite and olivine gabbro, with lesser oxide gabbro (7%) and minor MORB-type diabase intrusions (3%). The downhole lithostratigraphy is characterized by numerous inter-fingered intrusive bodies varying in thickness, with more evolved lithologies generally intrusive into less evolved lithologies (e.g., Blackman et al., 2006). Bulk Mg-values (Mg# (cationic ratio) = 100 × Mg/(Mg + Fe total )) vary from 60 to  Blackman et al., 2006). (b) Examples of recovered intervals of (left) Ol-T1 (Core 248R-3, 82-104 cm) composed of wehrlitic (WEHRL), troctolitic (TROCT) and minor dunitic (DUN) domains, and (right) Ol-T2 (Core 248R-2, 5-26 cm) mainly composed of plagioclase-dunitic domains and cut by a gabbroic vein. (c) Downhole composition of IODP Hole U1309D (from left to right): 20 m running average of rock type recovered (white indicates no recovery), variations in whole-rock Mg # (Mg# (cationic ratio) = 100 × Mg/(Mg + Fetotal); modified after Blackman et al., 2006), and Zircon Pb/U Ages where numbers indicate averages of measurements (Grimes et al., 2008). (d) Ni and Li average composition of olivine from single samples of Ol-Ts selected for this study in the interval 1100-1200 mbsf. 90 mol%, with local exceptions of lower values characterizing the most evolved lithologies (Fig. 1c). No systematic compositional trends are observed downhole (i.e., trend of fractional crystallization; Godard et al., 2009). Such complex structural and textural relationships, together with chemical evolution confined in discrete lenses and the variable downhole zircon ages indicate that melt injections occurred at random depths during the construction of the gabbroic sequence (Blackman et al., 2006;Godard et al., 2009;Grimes et al., 2008). The total duration of magmatic accretion was inferred to be~200 ka by radiometric dating of zircons (Grimes et al., 2008;Fig. 1c) and occurred over two main periods, with an older intrusive event forming the deepest interval from 600 mbsf to the bottom of the hole (~1.24 Ma average) and a younger event occurring throughout the core (~1.17 Ma average). Grimes et al. (2008) inferred that emplacement of the oldest and deepest 635 mthick lower crustal interval in Hole U1309D occurred within~150 kyrs at a continuous growth rate of 1.6 cm/yr. These estimations, together with the recognition of over 250 intrusive igneous contacts throughout the Hole, allowed assessing an average thickness of single magma sills of 10 m and their emplacement every 630 yrs (Grimes et al., 2008).

Hole U1309D olivine-rich troctolites
Ol-T is the dominant lithology in the interval between 1100 and 1300 mbsf (Fig. 1) and is locally very fresh with <1% serpentinization. Drouin et al., 2009, 2010and Suhr et al. (2008 interpreted Ol-Ts from the Atlantis Massif OCC as resulting from melt-rock interaction in a reactive porous flow process. Ferrando et al. (2018) demonstrated that the multi-stage reactive process is triggered by melt infiltration into mantle harzburgite. Reactive percolation of MORB-type melts leads to mantle olivine dissolution, which locally modifies the composition of melts. In turn, the tholeiitic crystallization suite is modified: concomitant crystallization of interstitial clinopyroxene and plagioclase occurs whilst olivine is re-equilibrating with the migrating melt (Drouin et al., 2009;Ferrando et al., 2018). Ferrando et al. (2018) posit that, as temperature decreases, locally modified melts partially crystallize in cross-cutting gabbroic veins and finally form part of the gabbroic intrusions building the crustal sequence at the Atlantis Massif OCC.
At the grain scale, during the open-system process of mantle-melt interaction, olivine is eroded (i.e., partial dissolution) and partially precipitated at crystal rims . The large variations and non-systematic downhole correlations of Ni and Li contents (Ni = 1800-2820 ppm; Li = 1.5-2.9 ppm) in olivine from Ol-Ts (Fig. 1d) are locally inherited from the precursor heterogeneous mantle (Fig. 10 in Ferrando et al., 2018). Ferrando et al. (2018) defined two endmember Ol-Ts on the basis of their structural and textural characteristics, and of their olivine composition. Ol-T1 (<77 vol% modal olivine; Table 2) is the most reacted end-member (i.e. more mantle olivine dissolution) characterized by olivine grains showing Mg#~85 mol% and Ni contents between 1870 and 2820 ppm, embedded in large oikocrysts of plagioclase and clinopyroxene ( Fig. 4e in Ferrando et al., 2018). Ol-T2 (>77 vol% modal olivine; Table 2) is the least reacted end-member with olivine displaying more evolved composition compared with Ol-T1 (Mg#~84 mol% and Ni = 1790-2130 ppm) and forming aggregates with interstitial plagioclase and minor clinopyroxene ( Fig. 4f in Ferrando et al., 2018). In both Ol-T1 and Ol-T2 olivine compositions are in equilibrium with adjacent clinopyroxene and plagioclase (Figs. 8-9 in Ferrando et al., 2018).
Geochemical profiles in olivine (Fig. 2, Supplementary Figs. S1a and S1b) were performed by Ferrando et al. (2018) and in this study along preferred crystallographic directions that were determined by Electron Backscatter Diffraction analyses (EBSD). These profiles record two distinct processes that involved chemical diffusive re-equilibration, first at magmatic temperatures (~1250°C) during melt-rock interactions and subsequently at decreasing temperature under subsolidus conditions (down to~1000°C) during cooling.
At magmatic temperatures, the fast diffusive transport of elements in olivine (e.g., , which overall shows less than three orders of magnitude difference in diffusivity (Table 1) . S1a and S1b). Flat profiles in olivine are interpreted to be the result of complete diffusive re-equilibration with the reacted and modified percolating melt during melt-rock interactions. These profiles can be used to quantify the minimum timescales over which chemical reequilibration likely occurred. This process of olivine re-equilibration is herein referred to as 'diffusive re-equilibration at magmatic temperatures (MT diffusive re-equilibration)'.
Although Ca in olivine was likely also completely re-equilibrated during melt-rock interactions, Ca profiles show lower Ca contents at the olivine crystal rim compared to the relative crystal core (Fig. 2). Similar convex profiles are documented in some olivine crystals also for Y and Heavy-REE (HREE, e.g., Yb; Fig. 2 and Supplementary Figs. S1a and S1b). Being preferentially hosted in clinopyroxene, Ca, Y and HREE can diffuse from olivine into clinopyroxene at subsolidus conditions due to the dependence of element partitioning on temperature decrease (e.g., Coogan et al., 2002;Witt-Eickschen and O'Neill, 2005). In turn, the rate of diffusion decreases at decreasing temperature (e.g. Coogan et al., 2005) hampering complete re-equilibration of olivine core with the relative olivine rim at subsolidus conditions. We thus infer that convex profiles do not result from the previous complete MT diffusive re-equilibration during melt-rock interactions (i.e., leading to flat profiles). They rather record subsolidus cooling after magma emplacement (see discussion below).

Geochemical profiles in olivine
The mineral modal contents, textures (grain size and phase relationships) and microstructures analysed by EBSD are described in Ferrando et al. (2018). We focused this study on core-rim geochemical variations in 42 olivine grains from samples of Ol-T1 and Ol-T2, and we report complementary analyses of in situ major and trace element compositions ( Fig. 2 and Supplementary Material Table S1). A subset of 12 most representative rim-core-rim profiles is shown in Supplementary Figs. S1a and S1b. Absolute concentrations of Yb are the highest of the Rare Earth Elements (REE) in olivine (Supplementary Material  Table S1), therefore we report Yb in Fig. 2 as representative of olivine REE concentrations.
Because element diffusion is overall anisotropic, geochemical profiles were collected along preferred directions selected parallel to (at least) one of the three principal crystallographic axis of olivine having axes plunge (measured from the sample surface) <5° . All three crystallographic directions of olivine were investigated in each sample. Core-rim and rim-to-rim traverses were measured with 19-80 μm spacing for major and minor elements, and 70-500 μm for trace elements (Supplementary Material Table S1).
Major elements and trace elements were determined using analytical instruments from the Microsonde Sud (Géosciences Montpellier, University of Montpellier). Major elements were measured by Electron Probe Micro Analyser (EPMA) using a CAMECA SX100 equipped with five wavelength-dispersive X-ray spectrometers (WDS); accelerating potential was set to 20 kV, beam current at 10 nA, and counting times were 30 s for all elements. Natural minerals and synthetic oxides are used as standards. In situ trace elements were measured after EPMA analyses to minimize loss of sample. Surface cleaning was performed on the selected EPMA spot by pre-ablation in 5-6 pulse per second. We used a Thermo Scientific Element 2 XR (eXtended Range) high resolution -Inductively Coupled Plasma Mass Spectrometry (ICPMS). The ICP-MS is coupled with laser ablation (LA) system, a Microlas (Geolas Q+) automated platform with a 193 nm Excimer Compex 102 laser from LambdaPhysik. An in-house modified 30 cm 3 ablation cell with a helium atmosphere was used to enhance sensitivity and reduce interelement fractionation (Günther and Heinrich, 1999). The laser energy density was set to 12-15 J cm 2 and repetition rate at 8-10 Hz. The laser spot size was 102-77 μm. Data were collected in time resolved acquisition mode with the background signal collected for 2 min followed by 1 min of sample ablation, or with 2.33 min for the blank and 40 s sampling. Concentrations were calibrated against the NIST 612 rhyolitic glass using the values given in Pearce et al. (1997). Data were reduced with the GLITTER software package (Van Achterberg et al., 2001) using the linear fit to ratio method, and 29 Si was used for internal standardization relative to EPMA data. Signals were carefully monitored and data were filtered for spikes on an element by element basis. Detection limits were between 0.07 and 1 ppm for Ni, <55 ppb for Mn, Cu, Zn and <22 ppb for Co; they were <5 ppb for most incompatible elements (i.e., Y and HREE). Reference basalt BIR-1G was used as internal standard to monitor accuracy as well as reproducibility within single series and between runs. This resulted in reproducibility better than 5% for Co and REE, and it is <12% for all other elements (see Ferrando et al., 2018 for details and Tables).

Numerical diffusion modelling
Diffusion in solids is described using Fick's law (J i = − D i ∂C i /∂x), which states that the concentration evolution of a chemical species i with time (flux J i ) is governed by its rate of diffusion (diffusion coefficient, D i , m 2 /s) and depends on the concentration gradient of i along a given direction x (∂C i /∂x), where C i is the concentration of i. Most numerical diffusion models are solutions of the Fick's second law (Crank, 1975). Several solutions of this equation have been published in the literature to evaluate the timescales of diffusion in minerals in various although simplified magmatic systems.
Flat profiles of Fo, Ni, Mn, Co and Zn were modelled to quantify the minimum timescales of MT diffusive re-equilibration ( Fig. 2, Supplementary Figs. S1a and S1b). To account for diffusion anisotropy in olivine, and to consider the influence of crystal size, morphology, and sectioning, MT diffusive re-equilibration was simulated using threedimensional (3D) diffusion models modified after Shea et al. (2015). Additionally, cooling rates were estimated using one-dimensional models (Ca-in-olivine geospeedometry; Dodson, 1973;Coogan et al., 2002Coogan et al., , 2007 to fit CaO convex profiles measured in olivine (Fig. 2, Supplementary Figs. S1a and S1b). The choice of using a 1D model for Ca diffusion at subsolidus conditions is supported by the negligible dependence of Ca diffusion rates on the crystallographic orientation of olivine (Coogan et al., 2005); the CaO profiles are convex in most olivine crystals and in all measured directions, thus indicating that crystal shape and diffusion anisotropy (3D effects) did not affect Ca diffusion at subsolidus conditions. Among CaO convex profiles, those displaying flat plateaus at the crystal core were preferentially selected for models of Ca-in-olivine geospeedometry in this study. We avoided asymmetric CaO convex profiles and those showing dipping plateau, as they are common in sections that are off centre and/or oblique to principal crystallographic axes and can contribute to the uncertainty in retrieved cooling rates (e.g., Shea et al., 2015).
Data in this study are from olivines that span the range of observed grain sizes from 1 to 4 mm. In the following sections we illustrate the details of MT diffusive re-equilibration models and Ca-in-olivine geospeedometry.

MT diffusive re-equilibration models
Chemical re-equilibration in olivine was simulated for eight elements that were chosen on the basis of their compatibility in olivine, from compatible (Fe\ \Mg [Fo], Ni) to moderately incompatible (Mn, Co, Zn, Ca) and highly incompatible elements (Y, Lu). Trace elements were investigated as common tracers of melt-rock interaction processes, during which their concentrations are strongly modified (e.g., Basch et al., 2018;Ferrando et al., 2018;Rampone et al., 2016;Sanfilippo et al., 2014). Moreover, the modelled elements cover a wide range of diffusion rates in olivine ( Table 1) allowing characterization of a broad range of minimum re-equilibration timescales.
Simulations were performed using three solutions of the 3D form of Fick's second law (Crank, 1975). Solutions for anisotropic diffusing species were applied for Fe\ \Mg, Ni, and Mn using the concentrationdependent equation (Eq. (A1)), and Ca using the non-concentration dependent form (Eq. (A2)); trace elements were modelled using the non-concentration dependent form for isotropic diffusion (Eq. (A3)) (see Appendix A for details on the equations). Because anisotropic diffusion of Co was documented exclusively at a single temperature condition (i.e., 1300°C; Spandler and O'Neill, 2010) and the Arrhenius relationship for D Co is only available for Co diffusion along the crystallographic axis [001] (D Co[001] , i.e., fast direction; Ito et al., 1999), we modelled isotropic diffusion of Co using D Co[001] in Eq. (A3). We are aware that this simplified approach possibly underestimates the timescales computed from Co profiles (i.e., use of the fastest D Co ), but we emphasize that the overall D Co is comparable to the D i of the other trace elements (see Table 1) thus leading to timescales in the same order of magnitude as Y and REE.
Although simple crystal geometries and analytical solutions are often used to calculate timescales of element diffusion, Shea et al. (2015) demonstrate that using spherical crystals lead to systematic overestimations of timescales. Following these findings, olivine is here modelled with polyhedral morphology (Shea et al., 2015) allowing to reproduce a realistic olivine shape (Fig. 3). The best preserved olivine crystals in the studied Ol-Ts have dimensions typical of euhedral olivine crystals with c-axis >> b-axis > a-axis and c-axis ≈ 2 * a-axis (Supplementary Fig. S2). Comparable olivine crystal shapes have been previously documented by Welsch et al. (2014) in samples from the same Hole U1309D. Therefore, we assumed a euhedral crystal similar to . Models were run twice, with 1 and 3 mm c-axis lengths. (b) Initial olivine (C 0 ) and melt (C i = C1) compositions used in the models. An initially homogeneous crystal with a sharp compositional boundary to the surrounding melt was then allowed to diffuse at the conditions specified for each model and element. olivines often described in erupted lavas (e.g., Lynn et al., 2017aLynn et al., , 2017b. The 3D olivine + melt models have dimensions of 221x221x221 voxels (a voxel being a pixel in three dimensions). The modelled olivine had dimensions of 201 voxels along the c-axis, 121 voxels along the b-axis, and 95 voxels along the a-axis (e.g., Shea et al., 2015;Fig. 3). Two crystal sizes reflecting the minimum and maximum olivine sizes observed in the thin sections were modelled with different voxel resolutions. For small grains c-axis was set~1 mm long (resolution of 4x4x4 μm per voxel; olivine crystal of dimensions x = 804 μm, y = 484 μm, z = 380 μm) and for coarser grains~4 mm long (resolution of 20x20x20 μm per voxel; olivine crystal of dimensions x = 4020 μm, y = 2420 μm, z = 1900 μm).
The advantage of the 3D model approach is that mineral core compositions can be tracked through time. This provides insights into how resilient original mantle olivine compositions are to being overprinted by secondary processes such as melt-rock interactions. A simplified estimate of diffusion distance using x = ffiffiffiffiffiffiffiffiffiffi D i t ð Þ p has been applied by numerous previous studies to infer the time (t) required to generate a diffusion profile of a given length (x) for an element of interest with known D i . However, crystals are 3D objects subject to diffusive fluxes in all dimensions. While this simplified 1D equation might broadly characterize the time required for diffusive re-equilibration to reach a crystal's core, our 3D models allow us to determine how long it will take to completely re-equilibrate that core composition. Thus, this modelling approach is essential to understand how quickly minerals re-equilibrate with their surroundings and cannot be estimated using the simplified 1D equation.
The continuous change of melt composition at the interface between crystal and melt during dissolution-precipitation processes would be best modelled using moving compositional boundary conditions (e.g., Chakraborty, 2008) combined with changes in crystal shape (Chakraborty, 2018). However, experimental investigations of reactive percolation of a primitive melt (MORB-type melt) through a dunitic (Borghini et al., 2018) or troctolitic matrix (e.g., Yang et al., 2019) documented that at magmatic conditions (T~1250°C; similar to U1309D Ol-Ts formations) dissolution-reprecipitation reactions of olivine are rapid, while Liang (2003) in an experimental and theoretical study of the kinetics of melt-rock reaction demonstrated that diffusion of multiple elements in basaltic melts (of the order~10 −11 m 2 /s; e.g., Kress and Ghiorso, 1995;Liang, 2003) is faster than diffusive re-equilibration in minerals (see Table 1). Also, we have no constraints on (i) the variation of melt volume due to dissolution-precipitation, (ii) the composition of olivine and melt at a given stage of the reactive process, and (iii) the effect of changes in crystal shape on diffusion rates. Therefore, for simplicity, we modelled MT diffusive re-equilibration assuming static crystal shape and static boundary conditions (e.g., constant melt composition). This assumption implies instantaneous change in melt composition (i.e., fast element diffusion in melt; e.g., Liang, 2003) produced by melt-rock interactions and melt-olivine equilibrium at the crystal rim.
Diffusion simulations were run until complete re-equilibration with the boundary condition was achieved. 20 models were saved at equal intervals throughout the diffusion simulation to track the evolution of the crystal's composition with time. These models were then sectioned through the olivine crystal's core perpendicular to the c-axis (Fig. 3). The composition of the central pixel in the 2D section (Fig. 3) was sampled at regular intervals throughout the model diffusion time to track the core compositional evolution. We also sampled 1D numerical traverses parallel to the a-axis in the ideal 2D section to illustrate the evolution of compositional profiles with time. The extent of MT diffusion is expressed in terms of % re-equilibration (noted % req; see Lynn et al., 2017b), which allows a direct comparison among diffusing chemical species with different absolute concentrations (e.g., ppm, wt%, mol%; Fig. 4a,c): where C initial is the composition of mantle olivine before the onset of diffusion, C measured is the composition of the olivine core after a given time (expressed in years, examples in Fig. 4b), and C equilibrium is the composition of olivine in chemical equilibrium with the surrounding melt. The magnitude of zoning in a 1D traverse parallel to the a-axis across olivine diminishes as time elapses and the profile becomes difficult to resolve analytically; distinctive inter-element zoning becomes less likely and the crystal nears "complete re-equilibration" beyond the kinetic window. We introduce the "effective % req" to take into account typical analytical uncertainties on each element (i.e., σ i = average analytical error of the modelled element i). The effective % req adjusts the reported timescales with the analytical profiles that could realistically be resolved in natural samples as they near 100% equilibrium with the boundary conditions. Using the effective % req, we consider complete re-equilibration when olivine core composition in the 3D model reaches C equilibrium + σ i (Table 3).
We ran the MT diffusive re-equilibration models at pressure of 2 kbar and oxygen fugacity at QFM (for Fe\ \Mg, Ni, Mn, and Ca), which are appropriate for the Atlantis Massif based on the estimated depths at which the gabbroic sequence was formed (~7 km; Grimes et al., 2008). Temperatures were set at 1230°C for models of Fo, Ni, Co, Mn, Ca (Fig. 4a), which correspond to the calculated crystallization temperature of plagioclase (Drouin et al., 2009). For the listed major and minor elements, models were also ran at 1210°C and 1250°C to account for uncertainties in the calculated timescales related to the +/− 23°C uncertainty of temperature estimate from Drouin et al. (2009). From these models at three different temperatures we calculated Δt T , which is the uncertainty in modelled timescale related to the uncertainty in the temperature estimates. The uncertainty at the high end of the T range (+Δt T ) is the timescale at 1210°Ctimescale at 1230°C, and the low end of the T range (−Δt T ) is the timescale at 1250°Ctimescale at 1230°C. As no temperature-dependent equation is available for calculating D Zn , D Lu and D Y (single value at 1300°C, Spandler and O'Neill, 2010; Table B1 in Appendix B), we also ran MT diffusive reequilibration models of all elements at 1300°C and fO 2 = 10 -8.3 bars (QFM-1; corresponding to the conditions of experiments in Spandler and O'Neill, 2010;Fig. 4c).
3.2.1.1. Diffusion coefficients. For the MT diffusive re-equilibration models, D i were calculated at appropriate magmatic temperatures, pressure and oxygen fugacity conditions using the Arrhenius relationships reported in Table B1. Uncertainties in D i are related to uncertainties in activation energies (~10%; e.g., Costa et al., 2008), in turn leading to uncertainties in computed timescales of diffusion. Di calculated at temperatures higher or lower than the temperature at which Fig. 4. Results of MT diffusive re-equilibration model to simulate timescales of melt-rock interactions at 1230 ± 20°C (a-b) and 1300°C (c) using C 0 a as starting olivine composition (Table B2). (a-c) Evolution of the modelled % re-equilibration (% req) of olivine core as function of time of re-equilibration (years). Symbols represent olivine grain size of 1 mm and lines are for grain size of 4 mm. Colors distinguish each modelled element and reported in the legend. The time of re-equilibration is compared with the time of emplacement of a single 10 m thick sill (~630 years, light grey) after Grimes et al. (2008). In (a) shaded areas in green, Co, and black, Ni, represent the ±Δ in time of re-equilibration (±Δt T ) at given % req related to uncertainties in crystallization temperature (Drouin et al., 2009). They are reported for Co and Ni in 1 mm models only to provide examples; ±Δt T of all other elements can be found in Table 3. (b) 1D rim-to-rim profile evolution along the a-axis across 1 mm sized olivine (see section '3.2.1 MT diffusive re-equilibration models' for procedure of 1D profile sampling) are reported for Fo, NiO (wt%) and MnO (wt%) from MT diffusive re-equilibration model at 1230°C.. Colored lines (elements are represented by colors as reported in the legend) highlight olivine composition after 50% of diffusive re-equilibration and olivine composition at effective % req (marked with * in the legend). For details on the MT diffusive re-equilibration models see the text and Appendix A. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) D i was experimentally constrained can have uncertainties of up to 3-4 orders of magnitude. In this study, D i were calculated at temperatures comprised in the temperature interval at which their Arrhenius relationships were experimentally calibrated. This minimizes the uncertainties in timescales, which are likely encompassed by Δt T and uncertainties for olivine initial composition (see later in this section). Also, greater uncertainties in timescales in these natural rocks are related to the variable grain size of olivine crystals that are here treated modelling different crystal sizes.
To account for Fo variations, D Fe-Mg(Fo) and D Mn were calculated using the diffusion coefficient equations from , and D Ni from Petry et al. (2004) (Table B1). As no such concentration-dependent equations are available for the Ca and Co, we calculated D Ca and D Co using the Arrhenius relationships (Coogan et al., 2005 for Ca andIto et al., 1999 for Co; Table B1) calibrated for olivine Fo~90.
Diffusion of REE is currently matter of debate and vast discrepancies exist between trace element diffusivities in olivine with rates that span over several orders of magnitude (10 −19 -10 −20 m 2 /s by Cherniak, 2010;~10 −15 m 2 /s by Spandler and O'Neill, 2010). This range may be related to the different melt-powder compositions used in the experiments and, specifically, on the activity of SiO 2 and Ti contents in the melt (Supplementary Fig. S3; see discussions in Jollands et al., 2016 andBurgess andCooper, 2013). The source of diffusant used by Spandler and O'Neill (2010) is a trace-element-doped MORB-type melt containing 12 wt% of Al 2 O 3 and 5 wt% of TiO 2 , in contrast to the synthesized mixture of REE (La, Dy, or Yb) aluminate and synthetic forsterite powders used by Cherniak (2010) that does not reproduce a silicate melt (TiO 2 = 0 wt%). On the other hand, the experiments of Spandler and O'Neill (2010) (i) also use a San Carlos olivine with doped trace elements at levels that do not reflect natural olivine similarly to Cherniak (2010), and (ii) were conducted only at single temperature, pressure and oxygen fugacity conditions, which prohibit the extrapolation to temperatures other than 1300°C. Moreover, REE diffusion rates have been suggested to also depend on the mechanism of diffusion and concentration of diffusing species in olivine (Chakraborty, 2018). Similar to findings of Li diffusion (Dohmen et al., 2010), in most natural systems REE may diffuse via two mechanisms that differ in diffusion rate (Chakraborty, 2018); though, dependence of the active diffusion mechanism on the concentration of the diffusing species has never been defined, nor diffusion coefficients for these different mechanisms have ever been quantified.
Because only the experiments of Spandler and O'Neill (2010) include the buffering of all chemical activities (e.g., aSiO 2 ) and because the computed melt composition in equilibrium with interstitial clinopyroxene (see Ferrando et al., 2018) contains~12-14 wt% of Al 2 O 3 we posit that the D REE,Y (for D REE only D Lu is available) values from Spandler and O'Neill (2010) are the most appropriate for modelling minimum timescales of re-equilibration after melt-rock reactions at Atlantis Massif. For this reason, we selected D Zn , D Lu and D Y single values from Spandler and O'Neill (2010). The use of fast diffusivities may lead to underestimation of timescales, and our model results represent minimum timescales of re-equilibration.
3.2.1.2. Initial and boundary conditions. Following the outcome of previous studies (see Section 2.2) on the reactive origin of Ol-Ts from assimilated mantle, we assumed that the initial composition of olivine is that of olivines in mantle harzburgites (C initial in Eq. (1)). Mantle harzburgites sampled at modern ridges are compositionally heterogeneous due to different extents of mantle melting and impregnation processes (e.g., D'Errico et al., 2016;Regelous et al., 2016;Tamura et al., 2008). Also, among the little data available of in situ trace element composition of olivines in abyssal harzburgite (e.g., D'Errico et al., 2016 ;Regelous et al., 2016) none is from samples collected at the Atlantis Massif OCC. The Gakkel Ridge is a slow-to ultra-slow avolcanic spreading ridge, similar to the region of MAR where the Atlantis Massif OCC is located. For these reasons, we explored the variation in timescales of MT diffusive re-equilibration using two different initial compositions (Table B2): (i) C 0 a is the average composition of olivine in harzburgites from Gakkel Ridge (D'Errico et al., 2016;Regelous et al., 2016); (ii) C 0 b is the composition of a single harzburgite olivine corresponding to that used by Ferrando et al. (2018) (values of Zn and Co are averages of the most primitive olivines from the global database of abyssal peridotite compositions in Regelous et al., 2016). The resulting variation in modelled timescale related to initial compositions is expressed as Δt x = │t(C 0 a) -t(C 0 b)│, where t(C 0 a) and t(C 0 b) are timescales calculated for C 0 a and C 0 b, respectively.
Boundary conditions for the models are set to the measured olivine compositions, which represent the last reacted melt composition that is in equilibrium with adjacent clinopyroxene and plagioclase present in Ol-T (C equilibrium in Eq. (1)). Boundary melt compositions are constrained by the averages of olivine composition in the U1309D Ol-T1 summarized in Table B2), as representative of the most reacted Ol-Ts.

Ca-in-olivine cooling speedometry
After all melt has crystallized, the temperature dependence of partition coefficients (K d ) between minerals (e.g., Coogan et al., 2002;Witt-Eickschen and O'Neill, 2005) controls the re-distribution of elements during subsolidus cooling at the rims of adjacent phases, and a corerim compositional gradient is produced within single crystals. The down-temperature dependence of Ca, Y and REE partitioning between olivine and clinopyroxene leads to the re-distribution of these elements from olivine rims to adjacent clinopyroxenes. Diffusion rates in olivine also decrease with decreasing temperature (e.g., Coogan et al., 2002;, resulting in incomplete chemical reequilibration of the crystal, in turn generating core-to-rim variations in slow-diffusing elements, such as Ca content. Calcium is often used to retrieve timescales of subsolidus cooling (e.g., Coogan et al., 2007;Faak and Gillis, 2016). In this study, we determine the cooling rates and closure temperatures for 7 selected olivine grains that record a progressive decrease of Ca, Y and REE contents toward the crystal rim, applying the Ca-in-olivine geospeedometry (Eq. (A4) in Appendix A) following the methods of Coogan et al. (2007). Here, the term 'closure temperature' describes the low temperature at which no appreciable Ca diffusion occurs in the studied olivine (measured using the geothermometer by Köhler and Brey, 1990; see Coogan et al., 2007). We used the D Ca from Coogan et al. (2005) (Eq. (A5) in Appendix A) as it is consistent with the most recent measurements by Bloch et al. (2019), who combined LA-ICP-MS, SIMS and atom probe (LEAP) analytical techniques. The Ca-in-olivine method assumes clinopyroxene is an infinite reservoir for Ca, thus we selected only olivines in contact with clinopyroxene. Because clinopyroxene has much higher Ca, Y and REE contents (here on average, CaO = 21 wt%) compared to those of olivine (here on average, CaO = 0.08 wt%), no core-to-rim variations are generally observed in clinopyroxene (i.e., infinite reservoir boundary condition).
Initial conditions for the diffusion model were homogeneous and set to the Ca content measured in the core of the olivine crystal ( Fig. 5 and Supplementary Fig. S4). Initial model temperatures were determined using the partition coefficient for Ca between olivine and clinopyroxene (Köhler and Brey, 1990) and the Ca content preserved in the olivine core plateau. This initial temperature was then used to calculate the initial D Ca used in the diffusion model (D Ca from Eq. (A5) in Appendix A). After one model iteration, the time elapsed (∂t) was used to determine the change in temperature (ΔT) as a function of an assumed cooling rate. We recalculated the partition coefficient, temperature, and diffusion coefficient for Ca after each time-step throughout the models to extract linear cooling rates; different cooling rates were imposed on the models until a best fit to the measured profile was obtained ( Fig. 5 and Supplementary Fig. S4). The total time elapsed represents only the diffusion time recorded by the zoning generated between the initial and closure temperatures. Thus, they are minimum estimates of the duration of high-T cooling histories.

Timescales of mantle assimilation into the crust
Timescales of diffusive re-equilibration after mantle-melt interactions are reported in Table 3 and in Fig. 4 at given % req. Because there is no chemical zoning in olivine (except Ca; Fig. 2

and Supplementary
Figs. S1a and S1b), we have no constraints on the time elapsed after complete re-equilibration (past the kinetic window) and therefore, the determined timescales are first order minimum durations. Timescales can change as a function of temperature, initial conditions and grain size. The related variability of MT diffusive re-equilibration timescales is reported in the following.
The 1230°C 3D modelling of 4 mm olivine grains using C 0 a as starting olivine composition show complete MT diffusive reequilibration (i.e., effective~100% req) within 220 years for all major and minor elements (Fe\ \Mg, Ni, Mn, Co, Ca; Table 3 and Fig. 4a). Using C 0 b as starting composition for the same set temperature and grain size provides timescales lower than 250 yr, leading to Δt X < 100 yr (Table 3). Models ran at 1210°C and 1250°C allowed us to test the variability in timescales related to the uncertainty of the geothermometric estimates (T = 1230 +/− 23°C). The 1210°C models yield the highest values of minimum durations of diffusive re-equilibration among all models. The variability in timescales (Δt T ) between models at 1230°C (taken as reference values) and at 1210°C and 1250°C ranges from +/− 5 yr to +/-40 yr between all investigated elements (Table 3). At any given T, initial composition and % req, the time of re-equilibration of an element decreases with decreasing grain size (from 4 mm to 1 mm in Fig. 4a,c), as equilibrium can be achieved more rapidly in smaller crystals.
The models of MT diffusive re-equilibration at 1300°C allow comparison between all investigated elements including Zn, Y and Lu (see 3.1 MT diffusive re-equilibration models). They yield shorter timescales compared to those simulated at 1230°C (Fig. 4) due to the Tdependence of D i (faster diffusivities at higher T). For olivine grains of 4 mm size and C 0 a as starting olivine composition, the 1300°C 3D models show that Zn displays the shortest re-equilibration time of  Table 3 Timescales (yr = years) of chemical re-equilibration, after mantle-melt interactions; % req is the effective % req (see text for details). Millimeters (mm) refer to the diameter of modelled olivine. Δt T is the difference between timescales calculated at 1210°C (t 1210 ), 1230°C (t 1230 ), 1250°C (t 1250 ) (+Δt T = t 1210 -t 1230 ; −Δt T = t 1250 -t 1230 ). Δt X (│Δt X │) is the difference between timescales calculated using two different starting olivine composition C 0 a and C 0 b.  (Table 3 and Fig. 4c). The construction of 2D sections (Fig. 3) also permits sampling of 1D rim-corerim profiles across the olivine models. With these data we can track the evolution of the composition of olivine zoning patterns at regular intervals throughout the 1 mm 3D models at 1230°C for representative major elements (Fo), and compatible (Ni) and moderately incompatible (Mn) minor elements (Fig. 4b). Elements are re-equilibrated at relatively constant rate up to 80% req (Fig. 4) of the olivine core. Fo content decreases and Ni and Mn increase until 80% relative to the starting composition (C 0 a) within 15 years (Fig. 4b).
To attain complete reequilibration other~10 years are necessary, showing that the reequilibration process progressively slows down as crystal nears complete re-equilibration (last~20% req.), due to decreasing compositional gradient (Fig. 4).
Regardless of the grain size, set temperature and initial composition, all elements spanning from Fe\ \Mg to incompatible REE (i. e., Lu) and Y show complete re-equilibration in olivine within 300 years (Table 3 and Fig. 4). We recall that in MT diffusive re-equilibration models from this study we assumed a polyhedral morphology of olivine (i.e., euhedral crystal). On the other hand, models of crystals with more equant aspect ratios, and thus ultimately a smaller volume, would result in even shorter re-equilibration timescales (e.g., Shea et al., 2015).
Overall, the computed minimum durations of re-equilibration are lower than the 630 yrs (light grey band in Fig. 4) estimated for the frequency of emplacement of 10 m-thick sills, and overall lower than the total~150 kyrs emplacement of the oldest and deepest 635 m-thick lower crustal interval (from 600 to 1235 mbsf) in the Atlantis Massif OCC (Grimes et al., 2008).

Timescales of subsolidus re-equilibration
Ca, Y and REE are preferentially hosted in clinopyroxene. Ca contents (in wt% from EPMA analyses) are systematically lower at olivine rim compared to the crystal core ( Fig. 2 Table S1 and Supplementary Figs. S1a, S1b). Y and Yb concentrations (in ppm from LA-ICP-MS analyses) are either constant in single olivines, or slightly lower at the crystal rim in comparison with the corresponding core (Fig. 2 Table S1 and Supplementary Figs. S1a, S1b). Nonetheless, because of serpentinization along olivine edges, most geochemical analyses were performed at a minimum distance of 50 μm from olivine grain boundary leading to a gap in data at the actual mineral/mineral interface. Moreover, the low concentration of trace elements in olivine force setting large spot sizes (102-77 μm in this study), which possibly prevent the analyses of compositional variations at the rim of the olivine crystal. As a consequence, we cannot preclude possible local chemical variations at olivine rims for some elements, which can explain the flat Y and Yb profiles measured in most discarded olivine grains and some selected olivines (Supplementary Figs. S1a and S1b).

and Supplementary Material
Core-to-rim chemical variations in olivine suggest that Ca, and to a lesser extent Y and Yb, diffused from olivine into the adjacent phases (i.e., plagioclase and clinopyroxene), as the effect of subsolidus diffusive re-equilibration driven by a temperature sensitive K d during cooling of the gabbroic sequence at the Atlantis Massif.
Partitioning of other elements, such as Mg\ \Fe and Ni, is also temperature-sensitive (e.g., Faak et al., 2013;Roeder and Emslie, 1970;Witt-Eickschen and O'Neill, 2005). Olivine acts as an infinite reservoir for compatible elements, nominally Ni and Mg, at decreasing temperature similarly to Ca in clinopyroxene. Cooling of gabbroic rocks can induce chemical re-distribution of Ni and Mg into olivine from clinopyroxene and plagioclase, respectively (Faak et al., 2013;Witt-Eickschen and O'Neill, 2005). If this was the case for olivines in Hole U1309D Ol-Ts, diffusion of Ni and Mg# from clinopyroxene and plagioclase into olivine would have resulted in convex profiles of Mg# in plagioclase (e.g., Faak et al., 2013;Sun and Lissenberg, 2018) and Ni in clinopyroxene during cooling. However, no zoning of Mg# or Ni contents is observed in clinopyroxene and plagioclase ( Supplementary  Fig. S5). Because element diffusion in clinopyroxene and plagioclase is slower than diffusion in olivine (D Cpx~1 0 −19 -10 −21 ; e.g., Van Orman et al., 2001;Zhang et al., 2010;D MgPlag =~10 −16 ;Faak et al., 2013), Ni and Mg# contents decrease at the very edge of plagioclase and clinopyroxene crystals without further concentration change (i.e., lack of element diffusion) toward the crystal cores. Such concentration variations (likely at the order of few μm) cannot be detected at the resolution of the laser spot size (77 μm), and therefore it is plausible that subsolidus re-equilibration of plagioclase and clinopyroxene was simply not detected in the studied Ol-Ts. On the olivine side, Mg# and Ni contents would increase of just a negligible amount because they are substantially more concentrated in olivine than in the adjacent phases (i.e., olivine infinite reservoir). At the resolution of EPMA and LA-ICP-MS used in this study, these chemical variations are also not detectable. Moreover, olivine, clinopyroxene and plagioclase should record equilibrium temperatures close to 1000°C if compositions were modified at subsolidus conditions (e.g., Coogan et al., 2002;Sun and Liang, 2014). Mineral chemical compositions indicate that these phases in Hole U1309D Ol-Ts are in equilibrium at temperatures between 1190°C and 1205°C with a melt showing Mg# = 58-65 mol% . The high equilibrium temperatures, together with the (apparent) lack of Ni and Mg zoning and olivine acting as infinite reservoir, indicate that flat profiles of elements preferentially hosted in olivine record magmatic processes, and that those concentrations in the studied olivine were likely not modified, or increased to a negligible extent, during cooling.
Initial temperatures recorded by olivine core Ca contents are on average~1200°C. Closure temperatures recorded by olivine rims range between 1014°C and 1042°C. To limit uncertainties on the computed closure temperatures, we also calculated equilibrium temperatures between olivine and clinopyroxene rims using the lattice strain model for REE and Y distribution among mantle minerals (Sun and Liang, 2014). These equilibrium temperatures range between 1000°C and 1050°C overlapping the computed closure temperatures. The duration required to generate a good fit to the measured Ca profile represents a minimum timescale over which cooling might have occurred. The higher-T and lower-T (< 900°C) history of cooling is outside the kinetic window of Ca, and thus the total duration of cooling is longer than that calculated here via Ca-in-olivine cooling geospeedometry ( Fig. 5 and Supplementary Fig. S4). The linear cooling rates obtained from the models span a range from 0.01 to 0.001°C/yr (Fig. 6), with a mean value of 0.004°C/yr ( Fig. 5 and Supplementary Fig. S4). Fig. 6a displays data from this study compared with downhole cooling rates of the oceanic crustal sequence from Hole U1309D determined by combined U\ \Pb zircon crystallization ages, (U\ \Th)/He zircon thermochronometry and multicomponent magnetic remanence data Schoolmeesters et al., 2012). These cooling rates were determined at lower temperatures, from 780°C tõ 250°C. All data (this study; Grimes et al., 2011;Schoolmeesters et al., 2012) overlap within error, showing that cooling rates are relatively constant with depth and time over a large range of temperatures from magmatic conditions to present (Fig. 6a).

Preserved mantle chemical heterogeneities
Reactive processes through oceanic crustal gabbros modify the composition of migrating melts and pre-existing crystal matrix (e.g., Boulanger et al., 2019;Lissenberg et al., 2013;Lissenberg and Dick, 2008). The spatial extent of these chemical modifications depends on the intensity of melt-rock interactions, which are in turn strictly related to the frequency of magma inputs and the melt/rock ratio (e.g., Basch et al., 2018;Higgie and Tommasi, 2012;Lambart et al., 2019;Fig. 6). Continuous melt infiltration at high melt/rock ratios and reactive melt transport can lead to progressive over-enrichments in incompatible elements (Godard et al., 1995) from the deepest section of lower oceanic crust to shallower levels (e.g., Hess Deep crustal sequence; Lissenberg et al., 2013). High melt/rock ratios enable chemical homogenization and re-equilibration of the protolith minerals with the reacted melt (Fig. 6b) at the scale of few tens of meters. Extensive melt-rock interactions are able to shift the composition of percolating melts toward the composition expected for fractionation of MORB at elevated pressures (Lissenberg and Dick, 2008). On the other hand, episodic inputs of magma and low melt-rock ratios integrated over time lead to changes in chemical composition of percolating melt and crystal matrix at local, centimeter-scale (Basch et al., 2018;Fig. 7c), having a significant effect on the geochemical budget of erupted basalts (Paquet et al., 2016).
The lack of up-hole over-enrichments (Fig. 1c) and the low-pressure magmatic differentiation signature of Hole U1309D gabbroic section suggest that percolating melts were overall not chemically modified by reactive porous flow at Atlantis Massif. Assuming continuous melt infiltration (e.g., beneath fast-spreading ridges), fast MT diffusive reequilibration timescales (i.e., at magmatic conditions) would be responsible of chemical homogeneity of olivines over tens of meters (Fig. 7b). In contrast, the composition of olivines is heterogeneous throughout Hole U1309D, as evidenced by Ni and Li (Fig. 1d), suggesting that the distribution of flow paths was confined into discrete intervals. The strong variations in Ni and Li are inherited from pre-existing chemical heterogeneity of the mantle protolith  and record local re-equilibration with the modified melt after the reactive process (Fig. 7c) at relatively high melt/rock ratios.
The presence of Ol-Ts at different depths throughout Hole U1309D (Fig. 1c), together with their inheritance of mantle heterogeneity and the fast MT diffusive re-equilibration of olivine crystals modelled in this study, point to episodic small-scale magmatism. Although estimated re-equilibration timescales are minimum durations, we document that Ol-Ts formed and cooled during the ongoing uplift of the Atlantis Massif OCC (see 6.3 Role of detachment faults on crustal cooling), suggesting that they remain shortly at depth after formation.

Depths of magma emplacement at Atlantis Massif
Multiple episodes of magmatic infiltration triggered partial assimilation and incorporation of lithospheric mantle slivers into the oceanic crust. Subsequently, the melts locally modified by mantle partial assimilation migrated upwards and crystallized at variable depths to form gabbroic intrusions. It is worth noting that modified melts likely represent only a portion of the parental magmas forming the crustal sequence. Injections and crystallization of magma preserving their primary composition (i.e., not modified by reactive processes) also participated to the building of the lower oceanic crust at the Atlantis Massif OCC. In the following, we use geothermometric estimates and cooling rates recorded in rocks from Hole U1309D to reconstruct the history of formation and cooling of the oceanic crust at the Atlantis Massif OCC.
Estimations of crystallization temperatures of the interstitial phases in Ol-Ts from the Atlantis Massif (Drouin et al., 2009;Ferrando et al., 2018) indicate that melt-rock interactions occurred at temperatures of 1230°C. At such magmatic temperatures, our MT diffusive reequilibration models (Table 3; Fig. 4) demonstrate that pre-existing mantle olivine crystals are completely re-equilibrated with the reactive percolating melt in either the time-frame of a single melt input (i.e. < 630 yrs. assuming the emplacement of a 10 m-thick sill; Grimes et al., 2008), or within the~150 kyrs of emplacement of the 635 m-thick deepest interval in Hole U1309D (i.e., between 600 and 1235 mbsf; Fig. 6. Cooling rates (°C/year) obtained from Ca-in-olivine geospeedometry on Ol-Ts from this study (orange dots, Ol-T1, and red dots, Ol-T2) compared with cooling rates of slowand fast-spreading oceanic crust. (a) Dots are downhole Hole U1309D cooling rates determined for gabbroic rocks by combined U\ \Pb zircon crystallization ages and (U\ \Th)/He zircon thermochronometry  and multicomponent magnetic remanence data (Schoolmeesters et al., 2012). For the latter, we selected initial cooling rates (780°C-250°C) to ignore the effect of late-stage hydrothermal circulation, which buffer temperatures and decreases the cooling rates to~0.0003°C/yr (Schoolmeesters et al., 2012). The dotted blue box indicates the depth interval in (b). (b) Zoom-in of the studied interval (1190-1196 mbsf). Dotted black lines indicate constant downhole conductive cooling during uplift toward the surface modelled assuming temperatures of 0°C at surface and 1300°C at depth of magma emplacement (model from Coogan et al., 2007). Models of conductive cooling are reported for uplift rates of 10 and 20 mm/year. (c) Ranges of cooling rates of the (i) slow-spread oceanic crustal section at Atlantis Bank OCC (black bars) retrieved using Ca-in-olivine geospeedometry (Atlantis Bank OCC and another OCC along the Mid Atlantic Ridge at 23°N (MARK area); Coogan et al., 2007) and thermochronometric data (John et al., 2004), and (ii) shallower section of the fast-spread crustal sequence at Hess Deep Rift (grey bar; Coogan et al., 2007 Faak andGillis, 2016;Sun and Lissenberg, 2018). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Grimes et al., 2008). Using 1230°C as formation temperature of Ol-Ts and the reconstruction of isotherm depths previously documented at segment ends along the MAR (Grimes et al., , 2008Schoolmeesters et al., 2012), we can predict the depth of Ol-Ts formation and mantle assimilation. The depth of the~900°C isotherm has been estimated to lay at 6-6.5 km beneath the ridge axis at Atlantis Massif (Grimes et al., , 2008Schoolmeesters et al., 2012). Considering a geothermal gradient of~300-350°C/km at these depths (e.g., Grimes et al., 2011;Oxburgh and Turcotte, 1968), we can extrapolate the 1230°C isotherm to lay at~7.5-8 km depth during active magmatism beneath the ridge (Fig. 8). This is consistent with the documented depth of clinopyroxene + plagioclase crystallization front (i.e. double saturation at 8 km) beneath the axis of a slow-spreading ridge (spreading rate~20 mm/yr; Hebert and Montési, 2010).
Dissolution-precipitation reactions within the mantle require relatively high temperatures for the reaction to prevail on melt crystallization (Kelemen, 1990); the efficiency of reactive processes thus decreases with decreasing temperature. Combining our estimations of subsolidus cooling rates (Ca-in-olivine: 0.004°C/yr; Fig. 6) and the OCC uplift rate of 15 mm/yr (Grimes et al., 2008), we estimate~700 m depth interval over which melt-rock interactions occurred and Ol-Ts formed and cooled from 1230°C to 1050-1100°C. Accordingly, we conclude that the 1050°C isotherm lies at 6.8-7.3 km depth (Fig. 8).
At depths shallower than ⁓7 km (i.e., T < 1050°C), melts are no longer modified by melt-rock interactions and likely start to crystallize the gabbroic sills (Fig. 8). Nonetheless, Grimes et al. (2008) and Schoolmeesters et al. (2012) estimated that most gabbroic sills are intruded in a rather cold lithosphere (850°C), at 5-7 km depth. Therefore, Ol-Ts must have moved vertically to cool conductively of additional 200°C (from~1050°C to~850°C) before being intruded by the gabbro parental magmas. Assuming that the position of isotherms was stable during magmatic activity at Atlantis Massif, we performed calculations of the vertical displacement using uplift rates (15 mm/yr) associated to cooling duration of the Ol-T, using the cooling rates modelled in this study (0.004°C/yr, Figs. 5 and 6). We estimate a vertical uplift of 750 m from 1050°C to 850°C, for a total of 1.5 km during cooling of Ol-Ts from magmatic temperatures (1230°C, 7.5-8 km depth) to cooler environments (850°C). This indicates that Ol-Ts were located at 6-6.5 km depth when gabbros formed (Fig. 8), consistently with the depth of gabbroic intrusion estimated by Grimes et al. (2008) and Schoolmeesters et al. (2012).
Taken as a whole, our results document an evolution from deep formation of Ol-Ts at 7.5-8 km depth (1230°C) to their subsequent cooling to 1050°C at 6.8-7.3 km during constant uplift. From~7 km depths at 1050°C, reactive processes forming the Ol-Ts cease, and melt modified by melt-rock interactions start crystallizing. The oldest intrusion event (1.24 Myr; Grimes et al., 2008) at Hole U1309D was identified in the lower interval of the 1415 m-deep gabbroic section (below ⁓600 mbsf). This suggests that the oldest intrusion event represents the deepest emplacement and early crystallization of melts that were previously modified by reactive processes. The continuous exhumation of the OCC by the detachment fault then brings the lower oceanic crustal section to 6-6.5 km depth (850°C), where younger and smaller intrusions (<10 m thick) are emplaced throughout the Hole U1309D (Blackman et al., 2006;Grimes et al., 2008) and the Ol-Ts are integrated within the gabbroic sequence.

Role of detachment faults on crustal cooling
Cooling rates retrieved from olivines in Ol-Ts from IODP Hole U1309D (Fig. 6a,b) are on average 0.004°C/yr (from Ca-in-olivine geospeedometry; Fig. 5 and Supplementary Fig. S4). These cooling rates are comparable with cooling rates measured at different depths throughout the oceanic crustal section from the same Hole, and obtained by U\ \Pb zircon crystallization ages combined with (U\ \Th)/ He zircon thermochronometry and multicomponent magnetic remanence data (Schoolmeesters et al., 2012;Fig. 6a). The combination of Ca-in-olivine geospeedometry from this study and previously published thermochronometric investigations document that cooling rates of the oceanic crust exposed at the Atlantis Massif OCC are constant with depth (Fig. 6a). Moreover, cooling rates of the studied oceanic crustal section are not only constant in space (i.e., over depth), but also in time, as evidenced by similar cooling rates from high-temperature (1230-850°C; 0.004°C/yr; this study) to low temperature (850-250°C; 0.002-0.004°C/yr; Schoolmeesters et al., 2012). Comparable downhole constant cooling rates estimated using Ca-in-olivine geospeedometry have been determined for gabbroic rocks from a similar, but compositionally different OCC along the ultraslow-spreading SWIR (Atlantis Bank OCC; Coogan et al., 2007;Fig. 6c). One plausible explanation for constant cooling rates over space and time provided by Coogan et al. (2007) lies in the efficiency of crustal uplift in removing heat conductively. Also, John et al. (2004) determined cooling rates of the oceanic crust at the Atlantis Bank OCC using a different approach based on thermochronometric data. By applying a 2-D plate-cooling model for oceanic lithosphere, they demonstrated that the oceanic crust at an OCC mimics conductive cooling (at Atlantis Bank OCC over the temperature range~900-330°C; John et al., 2004). Those cooling rates are comparable with cooling rates of the gabbroic sequence exposed at the Atlantis Massif OCC (Fig. 6c), which we thus infer to have cooled conductively.
Parental melts of the exposed gabbroic sequences (i.e., Atlantis Massif and Atlantis Bank OCCs) may have crystallized at depth within the lithosphere, and then cooled as the detachment fault transported the lower oceanic crust to shallower depths. In this scenario, conductive cooling would be simply controlled by the uplift rate (Coogan et al., 2007). To test this hypothesis we compare cooling rates of the Atlantis Massif with modelling of uplift-controlled cooling rates (modelling after Coogan et al., 2007;Fig. 6b). Constant cooling rates of the studied oceanic crust can be reproduced by conductive cooling during constant uplift of the oceanic crust at rates of 10 to 20 mm/yr (Fig. 6b). This is consistent with the average vertical transport rate of 15 mm/yr estimated at the Atlantis Massif OCC by zircon thermochronometry (Grimes et al., 2008), and predicted by numerical models of oceanic crust accretion accommodated by detachment faulting at slowspreading ridges (i.e., OCC; Tucholke et al., 2008). This simplified model suggests that denudation by detachment fault is able to cool conductively the oceanic crust.
The constant cooling rates point to continuous vertical transport of Ol-Ts induced by the Atlantis Massif detachment fault (e.g., Canales et al., 2008;Escartín et al., 2008;Ildefonse et al., 2007;Tucholke et al., 2008). The progressive uplift transported Hole U1309D crustal sequence from depths of~6 km (isotherm of 850°C) to~1.5 km where gabbros passed through the~250°C isotherm (zircon thermochronometry estimations; Grimes et al., 2011;Schoolmeesters et al., 2012). Simple exhumation calculations based on the downhole constant cooling rates (0.003-0.004°C/yr) and the constant exhumation rate of 15 mm/yr indicate that the oceanic crustal sequence has indeed been uplifted~4-4.5 km during its cooling from 850°C to 250°C (Fig. 8). The oldest gabbroic intrusions were likely uplifted after their emplacement at depth and their cooling was driven by detachment fault-related exhumation. Uplifting may not represent the sole mechanism of heat removal: magma emplacement in a 'cold' lithosphere (<850°C) could account for fast and downhole constant cooling (e.g., Sleep and Warren, 2014) for the youngest magma injection event.
Lower oceanic crust exposed at other OCCs, such as the Atlantis Bank (SWIR, 57°E; Coogan et al., 2007;John et al., 2004;Fig. 6) and along the Fifteen-Twenty Fracture Zone (MAR, 14°N-16°N; Grimes et al., 2011), also record fast and constant downhole cooling rates (0.002-0.005°C/yr). The latter are comparable to our estimates of cooling rates of the gabbroic sequence at the Atlantis Massif OCC (0.004°C/yr on average). These similarities suggest that denudation of oceanic crust by long-lived detachment faults is able to remove heat from the slow-spreading lower oceanic crust worldwide.
Hydrothermalism is widely described along slow-spreading ridges (e.g., Boschi et al., 2008;Früh-Green et al., 2016) and also affects cooling of the oceanic crust. Heat removal by fluid circulation is well constrained in the upper section of lower oceanic crust at fastspreading ridges (Hess Deep, East Pacific Rise) and analogous ophiolites (Samail Ophiolite, Oman), where very fast cooling rates are observed ( Fig. 6c; Coogan et al., 2007;Faak and Gillis, 2016;Sun and Lissenberg, 2018). They are significantly higher than those observed at OCCs. At OCCs, hydrothermal circulation appear to dominantly affects rocks nearest to the detachment fault where higher permeability allows fluid pathway (Hirose and Hayman, 2008); they undergo slow cooling (~0.0003°C/yr; Schoolmeesters et al., 2012) from 250°C to present day temperature. However, the presence of high-temperature amphiboles downhole Hole U1309D (Blackman et al., 2006) suggests that hydrothermal convection could contribute to heat removal at higher temperatures deeper in the section.

Concluding remarks
Retrieving the timescales of magmatic processes leading to the formation of oceanic crust is crucial for constraining magma supply at ridges and to understand the control of exhumation mechanisms on cooling of the oceanic crust. We have shown that flat geochemical profiles in olivine can be used to calculate minimum duration of complete re-equilibration of mantle olivine matrix with the melt modified during melt-rock interactions. Despite mantle olivine can be completely re-equilibrated in rather short timescales, textural and chemical heterogeneities of the precursor mantle are preserved at the cm-scale throughout Hole U1309D. Such preserved smallscale heterogeneities rule out continuous magma supply and indicate that magmatism at the Atlantis Massif OCC occurs as discrete melt inputs. In this context, our minimum re-equilibration times (<300 years) estimated using multiple elements demonstrate that mantle assimilation and formation of Ol-Ts can occur in the timeframe of a single melt input. The lithospheric mantle is incorporated in the lower oceanic crust within the~150 kyrs of emplacement of the oldest and deepest magmatic interval in Hole U1309D. Subsequent subsolidus re-equilibration during cooling of the gabbroic sequence from magmatic temperatures (1230°C) down to 1050°C was constrained for the first time at the Atlantis Massif OCC in this study; we estimated cooling rates of 0.003-0.004°C/yr. These timescale estimations, combined with cooling rates at lower temperatures (from 850°C to 250°C) allow us to reconstruct the geodynamic evolution of Ol-Ts from their formation at depth to their integration into the oceanic crustal sequence during continuous exhumation. At Atlantis Massif, magma inputs at depth (7.5-8 km) assimilate discrete intervals of lithospheric mantle through reactive porous flow processes, forming Ol-Ts. As detachment faulting initiates, Ol-Ts are transported to shallower depths (at~6-6.5 km) and concomitantly cool down to~850°C. In this rather cold environment, the reacted melts that previously formed the Ol-Ts start crystallizing. Subsequently, the gabbroic sequence is continuously uplifted at a rate of~15 mm/yr and cools conductively at constant cooling rates.
The evolution of oceanic crust at OCC is controlled by the decreasing temperature of the uplifting system. This decrease in temperature governs the rheology of the host rock and rules the ability of the melt to react at depth (T > 1050°C) and to segregate into magmatic intrusions at lower temperatures (T~850°C). The constant downhole cooling rates from magmatic temperatures down to 250°C indicate that the gabbroic sequence continuously formed and progressively cooled over~150 ka, with no apparent change in spreading and exhumation rates.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Simulations of MT diffusive re-equilibration were performed using three solutions of the three-dimensional form of Fick's second law (Crank, 1975). The concentration dependent equation for anisotropic diffusing species was used for Fe\ \Mg, Ni, and Mn where C i is the concentration of element i (ppm or wt%), t is time (s), and D x , D y and D z are the diffusion coefficients (m 2 / s) of element i along the principal crystallographic axes [100], [010] and [001], noted in x, y, and z dimensions, respectively. The MT diffusive reequilibration models of Ni and Mn were run concurrently with Fe\ \Mg models to account for the concentration dependence of these elements on the olivine Fo content. Diffusive re-equilibration of Ca was modelled using the non-concentration dependent form for anisotropic diffusion: and trace elements (Co, Dy, La, Y, Lu, Zn) were modelled using the nonconcentration dependent form for isotropic diffusion: For the Ca-in-olivine geospeedometry we used a non-concentration dependent 1D form of Fick's second law: We utilized the simplified diffusion coefficient for Ca in olivine (D Ca ) from Coogan et al. (2005): where D o is the preexponential factor (m 2 /s), E is the activation energy (kJ/mol), R is the gas constant, T is temperature (K), and fO 2 is the oxygen fugacity.
Element C 0 a (initial condition) C 0 b (initial condition) C1 (boundary condition) -averages of EPMA and LA-ICP-MS measurements in Ol-T 1