Numerical model of crustal accretion and cooling rates of fast-spreading mid-ocean ridges

We designed a thermo-mechanical numerical model for fast-spreading mid-ocean ridge with variable viscosity, hydrothermal cooling, latent heat release, sheeted dyke layer, and variable melt intrusion possibilities. The model allows for modulating several accretion possibilities such as the “gabbro glacier” (G), the “sheeted sills” (S) or the “mixed shallow and MTZ lenses” (M). These three crustal accretion modes have been explored assuming viscosity contrasts of 2 to 3 orders of magnitude between strong and weak phases and various hydrothermal cooling conditions depending on the cracking temperatures value. Mass conservation (stream-function), momentum (vorticity) and temperature equations are solved in 2-D cartesian geometry using 2-D, alternate direction, implicit and semi-implicit finite-difference scheme. In a first step, an Eulerian approach is used solving iteratively the motion and temperature equations until reaching steady states. With this procedure, the temperature patterns and motions that are obtained for the various crustal intrusion modes and hydrothermal cooling hypotheses display significant differences near the mid-ocean ridge axis. In a second step, a Lagrangian approach is used, recording the thermal histories and cooling rates of tracers travelling from the ridge axis to their final emplacements in the crust far from the mid-ocean ridge axis. The results show that the tracer’s thermal histories are depending on the temperature patterns and the crustal accretion modes near the mid-ocean ridge axis. The instantaneous cooling rates obtained from these thermal histories betray these discrepancies and might therefore be used to characterize the crustal accretion mode at the ridge axis. These deciphering effects are even more pronounced if we consider the average cooling rates occurring over a prescribed temperature range. Two situations were tested at 1275–1125 C and 1050–850 C. The first temperature range covers mainly the crystallization range that is characteristic of the high temperature areas in the model (i.e. the near-mid-oceanic-ridge axis). The second temperature range corresponds to areas in the model where the motion is mainly laminar and the vertical temperature profiles are closer to conductive. Thus, this situation results in less discriminating efficiency among the crustal accretion modes since the thermal and dynamic properties that are described are common to all the crustal accretion modes far from the ridge axis. The results show that numerical modeling of thermo-mechanical properties of the lower crusts may bring useful information to characterize the ridge accretion structure, hydrothermal cooling and thermal state at the fastspreading ridges and may open discussions with petrological cooling rate results.


Introduction
There are remaining uncertainties in how the oceanic crust is accreted at fast mid-ocean ridges, both in terms of accretion geometry but also in terms of roles and efficiencies of the cooling processes like hydrothermal convective circulation.During the last decades, three main families of structures have been proposed to take into account the local thermal, seismic or geophysics properties of the mid-ocean ridges.Thus, Norman Sleep (1975) proposed a gabbro glacier (G structure) mechanism where crystallization occurs below the sheeted dykes at the floor level of a shallow melt lens.Gabbros flow downward and outward to build the entire lower oceanic crust below the sheeted dykes.This ridge structure P. Machetel and C. J. Garrido: Numerical model of crustal accretion is compatible with the geophysical observations collected at the East Pacific Rise (EPR) and also with the structural studies of the Oman ophiolite (Nicolas et al., 1988;Kent et al., 1990;Sinton and Detrick, 1992, Henstock et al., 1993, Nicolas et al., 1993;Quick and Denlinger, 1993;Phipps Morgan and Chen, 1993).However, it seems that other geophysical measurements at the EPR (Crawford and Webb, 2002;Dunn et al., 2000;Nedimovic et al., 2005) and field observations from the Oman ophiolite (Kelemen et al., 1997;Korenaga and Kelemen, 1997) also suggest mixed accretion mechanisms (M structures) that involve melt lenses at both shallow depth and Moho Transition Zone (MTZ) (Boudier and Nicolas, 1995;Schouten and Denham, 1995;Boudier et al., 1996;Chenevez et al., 1998;Chen, 2001).Furthermore, several authors (Bédard and Hebert, 1996;Kelemen et al., 1997;Korenaga and Kelemen, 1997;Kelemen and Aharonov, 1998;MacLeod and Yaouancq, 2000;Garrido et al., 2001) argue in favor of melt intrusions at various depths through superimposed sills at the ridge axis between the Moho and the upper lens (S structure).
A lively scientific debate opposes these three possibilities according to their effects on the thermal structure near the ridge.Indeed Chen (2001) argued against the M and S propositions observing that the latent heat release during crystallization would melt the lower crust if a significant quantity of gabbro were generated deep in the crust without efficient extraction of heat from the hot/ductile crust by efficient hydrothermal cooling.However, the scientific debate about the depth and the temperature for which hydrothermal cooling remains efficient is still pending with opinions that it may be as deep as Moho, consistently with the seismic observations at EPR (Dunn et al., 2000).The depth where hydrothermal convective processes cool the lower crust is related to the thermal cracking temperature of gabbros that depends on cooling rates, grain sizes, confining pressures and viscoelastic transition temperatures (DeMartin et al., 2004).A review of this question has been proposed by Theissen-Krah et al. (2011), with a cracking temperature ranging from 400 to 1000 • C. Hence, analyzing high temperature hydrothermal veins (900-1000 • C), Bosch et al. (2004) found hydrous alterations of gabbro active above 975 • C, requiring hydrothermal circulations until the Moho.Their results are in agreement with those of Koepke et al. (2005), which have proposed hydrothermal activity of the deep oceanic crust at very high temperature (900-1000 • C) and with those of Boudier et al. (2000), which proposed a temperature cracking higher than the gabbro solidus.Conversely, Coogan et al. (2006), bounds the hydrothermal flows in the near-axis plutonic complex of Oman ophiolite at a temperature of 800 • C, in agreement with the model of Cherkaoui et al. (2003).This very exciting scientific debate motivated us to improve a tool that may be useful exploring the effects of the melt accretion structures and hydrothermal cooling hypotheses on the nearridge thermal and dynamic patterns.
In this work we present a series of cases calculated from a numerical code written to explore the sensitivity of the thermal and dynamic patterns of the mid-oceanic ridge versus the hydrothermal cooling efficiency and the crustal accretion mode.However, in spite of their common features, each ridge is also a particular case, with its own local properties and configurations.This diversity, combined with the still pending debates among the scientific community about the ridge structures, the location and the strength of the hydrothermal cooling, the value of the cracking temperature values, the viscosity contrast in the upper crust and the feedbacks between these uncertainties, requires broad parameter studies that are sometimes difficult to present synthetically.Our code has been designed to allow broad, easy (and quite cheap in terms of computer time -from a few hours to a few days on a modern PC) explorations of these assumptions.This paper presents a study of the effects of the three most common melt intrusion geometries, with modulation of viscosity and hydrothermal cooling according to depth and cracking temperature.The Fortran sources (and data files), which are easy to customize for whoever may be interested, may be modified to focus on other effects as the assessments of the feedback interactions between these different processes that may help to understand the physical behaviors occurring at mid-oceanic ridge and may help to interpret the petrologic or geophysical observations.
We have chosen to simulate hydrothermal cooling thanks to an enhanced thermal conductivity triggered by threshold cracking temperatures.Following Chenevez et al. (1998) and Machetel and Garrido (2009), the numerical approach consists of iteratively solving the temperature and motion equations.They are linked by the non-linear advection terms of temperature equation and by variable viscosity.In their work, Chenevez et al. (1998) did not open the possibility of melt intrusion other than the M structure.Later thermal models by Maclennan et al. (2004 and2005) assumed melt intrusions at several positions rather than at the bottom of the ridge axis, but without explicit coupling between motion and temperature.The present work follows the one of Machetel and Garrido (2009) that opened the possibility of modulating the melt intrusion structure with a consistent coupled solving of motion and temperature equations.However, this work didn't take into account the sheeted dyke layer that strongly modifies the thermal structure of the lower crust at shallow depth and was inducing, in certain cases, unrealistic flow patterns at the ridge axis.This flaw is corrected in the present work that also describes the physical coupling of hydrothermal effects with temperature structure, crystallization, viscosity and crustal accretion mode.

Theoretical and numerical backgrounds
A global iterative process couples the temperature and motion equations in a Eulerian framework until reaching Geosci.Model Dev., 6, 1659-1672, 2013 www.geosci-model-dev.net/6/1659/2013/steady-state solutions.The latter are used, in a second step in a Lagrangian approach, to compute the thermal histories of tracers along their cooling pathway in the lower crust.If the global principles, the basic equations and the numerical methods described by Machetel and Garrido (2009) seem similar, the new internal and boundary conditions used to take into account the sheeted dyke structure induced deep modifications in the program.Thus, the index of the computation grid have been subjected to several row and column displacements in order to take into account the dyke injection and the horizontal melt injection at the upper lens level.Such calculations were not possible with the code given by Machetel and Garrido (2009).In the present one, the 2-D computation area depends on two indexes (ranging from 1 to N x , horizontal direction; and from 1 to N y , vertical direction).The 2-D Laplacian operators for motion (within a stream function -vorticity framework) and temperature are solved with finite-difference alternate scheme that leads to inverse tri-diagonal matrix.Each row or column of the computation grid can be split up into to five horizontal or vertical resolution segments of at least three computation nodes and which extremities may potentially be used to prescribe external boundary conditions (if the node is located on the external boundary of the grid) or, if the point is inside the grid, internal conditions.These local properties are defined in the subroutine "computation_grid" that can be carefully modified to further developments of the code if some modifications as finite lengths of sills or melt lens geometries had to be taken into account.Only three crustal accretion modes have been described in details in the paper: (1) the Gabbro Glacier, (2) the mixed MTZ and shallow lenses, and, (3) the superimposed sill configuration.They can be considered as benchmark cases to initiate numerous situations where internal conditions are applied to explore the effects of injection geometries including asymmetric melt intrusion at the ridge or asymmetric expansion of plates.
In this paper, we simulate the thermal and dynamical behaviors of the crustal parts of two symmetrically diverging lithospheric plates V p = 50 mm yr −1 , and symmetric crustal accretion.Three conservative principles have been applied to ensure the momentum, mass and energy conservations of the fluid including the effects of melt intrusions, latent heat releases, viscosity variations and hydrothermal cooling.Within this framework the Boussinesq, mass conservation equation can be written as Eq.(1).
The mathematical properties of this zero-divergence velocity field allow introducing a stream-function, ψ, from which it is possible to derive the velocity components (Eq.2).
The stream-function approach ensures a mathematical checking of the zero divergence condition for the velocity field.It also allows accurate, easy to operate local prescriptions of discharges for melt injection.Indeed, thanks to the mathematical and physical meanings of the stream-function, the difference of values between two points measures the flux of matter that flows in the model through that section.Furthermore, the stream-function contour maps reveal the tracer trajectories that offer direct visualizations of the flow analogous to virtual smoke flow visualization (e.g.Von Funck et al., 2008).
Beyond the equation of continuity, the fluid motion must also verify the conservation of momentum through an equation that links the body forces, the pressure and the stress tensor with the physical properties of the lower crust.From the physical values given in Table 1, the Prandtl number, Pr = (c p η)/k that characterizes the ratio of the fluid viscosity to the thermal conductivity, ranges between 2.5 × 10 15 and 2 × 10 18 .Such high values allow the neglect of the inertial terms of the motion equation, which finally reduces to Eq. (3), where the stress tensor, τ , can be written as Within this context, we also assume that the crustal flow remains two-dimensional in a vertical plane parallel to the spreading direction.Then, the vorticity vector ω = curl v reduces to only one non-zero component equivalent to a scalar value ω.After some mathematical transformations, it is possible to rewrite the continuity and momentum equations (Eqs. 1 and 3) in a system of two coupled Laplace equations for stream-function and vorticity (Eqs.5 and 6) This formulation allows keeping all the terms that appear with the effects of a non-constant viscosity.
Finally, the vorticity and stream-function equations need to be completed with an equation ensuring the conservation energy.The left-hand side of the temperature equation (Eq.7) includes advection of heat through the material derivative while its right-hand side takes into account the hydrothermal cooling (through an enhanced thermal conductivity), the heat diffusion, the latent heat (according to the variations of the crystallization function c described later) and the heat produced by viscous heating.
The used for non-linear terms of partial derivative equations.This introduces constraints on time stepping, which, for temperature equation, follows the Courant criterion while it is overrelaxed for the stream-function and vorticity elliptic operators.The formalism described above offers two main advantages.The first is the numerical stability of the coupled ellip-tic, Laplace operators for ω and ψ.The second is the physical meaning of the stream function since its amplitude difference between two points of the grid measures the discharge (m 3 s −1 m −1 ) of matter flowing between.Furthermore, as a direct consequence of its mathematical relationship to velocity (Eq.2), the velocity vectors are tangent to the lines of the contour map of ψ.
To set the geometry of the melt intrusion at the ridge axis we consider the mass conservation properties that stipulate the balance of inflows (at the bottom ridge) and outflows (through the right and left lateral boundaries) (Fig. 1).Since, in this paper, the left and right spreading rates are equal, the amplitude of the stream-function jump at the bottom of the ridge axis must be ψ c = 2V p H .The stream-function is mathematically defined as an arbitrary constant that allows setting, once in the grid, an arbitrary value.We have chosen to impose ψ lb = 0.5ψ c at the bottom left boundary of the grid.Then, starting from this value, and turning all around the box, we can determine the external boundary conditions that respect the flows escaping or going in the computation box through these boundaries (Fig. 1).At the left lateral boundary the stream-function decreases linearly: (ψ l = ψ lb −y•V p ) to reach ψ t = 0 at the top.Since no fluid can escape through the upper boundary, the stream-function condition remains zero from going from the left to the right upper corners.Conversely, starting from this left upper corner, the value of the stream function on the lateral right boundary decreases linearly from zero to ψ r = V p (y − H ) at depth y, and to ψ lb = −0.5ψc at the lower right corner.At the ridge, the crustal accretion mode will be set thanks to internal conditions applied to the internal values of the stream function on the two central columns of the computation grid (Fig. 1).This method does not affect the Alternate Direction Implicit (ADI) solving of the temperature and motion equations since it is easy to keep the tri-diagonal forms of the inversion matrix by splitting the horizontal or vertical resolution segments into shorter ones surrounding the location where the internal condition has to be applied.Then both ends of the segments are used by the algorithm as internal or external boundary conditions.The amplitudes of stream-function jumps, located on the two central columns of the computation grid, from the Moho to the upper lens level, now simulate the discharges of sills and lenses defining the hypothesized melt intrusion pattern (Fig. 1).Free-slip boundary conditions (ω = 0) have been applied to the vorticity equation.
However, it is also necessary to add bottom, lateral and internal conditions to solve the temperature equation.We have considered that the oceanic crust is embedded in an halfspace cooling lithosphere (Eq.8) with a ridge temperature equal to the melt intrusion temperature This choice of using half-space cooling thermal boundary conditions has been done to limit as much as possible the arbitrary prescriptions on thermal crustal surrounding, although no direct measurements of temperature or dynamic states are available at MTZ.This approach follows the global geodynamic agreement that considers the half-cooling law as a consistent first-order approximation of temperature in oceanic lithospheres.It is also justified by the fact that, in the present study, we focus our attention on the thermal and dynamic properties near the ridge axis, far from the lateral boundaries, which effects are expected to be significantly weaker than those of accretion structures and hydrothermal cooling.Indeed, while, near the ridge, steep vertical thermal gradients ensue from hydrothermally enhanced heat extraction and motion is deeply influenced by the mass conservation and the crustal accretion mode, far from the ridge, the crustal evolution is mainly driven by thermal conduction and laminar motions.It is clear that this assumption generates side effects near the lateral boundaries of the computation area.However, such situations probably exist in nature where hydrothermal cooling induces near the ridges axis enhanced vertical heat extraction that may exceed the heat transport by conductive processes alone, resulting in areas cooler than the rest of the crustal part of the lithosphere.
The injection of the sheeted dyke layer at the ridge axis is simulated, at the roof of the upper lens, by prescribing internal conditions on the stream function values.Their amplitude differences correspond to the fluxes of matter required to build the left and right sheeted dyke layers (Fig. 1).Above the lens, linear stream function increases are also prescribed from the roof of the upper lens to the surface on the two central columns of the computation box.Then, starting from these two central points and going respectively to the right and left lateral boundary, the stream function is maintained constant all along the horizontal rows corresponding to the sheeted dyke layer.This ensures, in the sheeted dyke layer, horizontal velocities equal to the plate spreading and zero vertical velocities.Consistently with the solving of temperature in the remaining of the solution, the thermal behavior in the sheeted dyke layer is obtained from Eq. ( 7), in which the vertical advection term is cancelled by the zero value of the radial velocity that ensures a conductive vertical heat transfer through the sheeted dyke layer while it is conductive and advective in the horizontal direction.
Second-order accuracy and ADI finite-difference schemes have been applied to solve the Laplace operators for streamfunction, vorticity and temperature (e.g.Douglas and Rachford, 1956) with a 100 × 600 node grid corresponding to a 6 × 40 km oceanic crust area.The computational process iteratively solves the vorticity, stream-function and temperature equations until the maximum relative evolutions of the temperature between two time steps falls below 10 −7 at each node of the computational grid.

Introducing the crust physical properties
We also need to describe the links between viscosity, thermal conductivity, hydrothermal cooling and temperature, continuity and motion equations.In order to minimize accuracy losses that could result from numerical differentiation, the   local variations of the physical parameters have been written as hyperbolic, tangent-like step functions (Eq.9).Such functions, their derivatives and potencies are continuous and display accurate analytical expressions, evolving from 0 to 1.With this formalism, 88 % of the transition occurs over a 2δ range, centered on a threshold value, d T .The quantity d may stand either for distances, temperatures or crystallization.Table 1 recalls the characteristic values used for the various physical parameter of this study.Even if most of them are subjects to still pending scientific debates, they are representative of oceanic crust properties.In any cases, they allow for the illustrating of the sensitivity of the numerical model to the assumptions and uncertainties about the physical values and processes that affect the ocean ridge dynamics.
Hence, following the work of Kelemen and Aharonov (1998) on the sharpness of the crystallization, we will consider that the melt fraction varies rapidly around a threshold temperature, T c = 1230 • C, with a transition width δT c = 60 • C (Eq. 10 and Fig. 2-top).In our study, we have chosen to link the viscosity to crystallization through Eq. ( 11) because several authors have emphasized such steepness for the viscosity variations versus crystallization (e.g.Pinkerton and Stevenson, 1992;Marsh, 1996Marsh, , 1998;;Ishibashi and Sato, 2007).In any case, the hypotheses done on the viscosity can be easily modified in the subroutine called "viscosity" of the numerical code where small changes in the description of physical dependencies of viscosity can be taken into account without changing the motion resolution scheme.The global iterative process used in our approach couples temperature, vorticity and stream function by successive solving.Therefore viscosity, which explicitly depends on temperature through Eqs. ( 10) and ( 11), evolves with it all along the computing process.
The present work assumes that viscosity ranges between a strong, cold phase for 0 % melt fraction (η s = 5 × 10 15 Pa s) and a weak, hot phase for 100 % melt fraction (η w = 5 × 10 12 Pa s or η w = 5 × 10 13 Pa s).These viscosity values may seem low but they belong to the range used by Chenevez et al. (1998)   effects of its variations.The two curves in the lower panel of Fig. 2 displays the resulting viscosity-temperature relationship according to the contrast of viscosity assumed between the strong and weak phases.Cases involving the two viscosity contrasts, 5 × 10 12 versus 5 × 10 15 and 5 × 10 13 versus 5 × 10 15 (Pa s) have been calculated but we will not show all the results (except later in Fig. 6) because they are so close that they are impossible to distinguish with the naked eye in the figures.
Cryst (x, y) = Hyperbolic, tangent-like step functions have also been used to simulate hydrothermal cooling by linking the enhancement of thermal conductivity to depth and temperature.First we consider, through Equation 12, that k decreases with depth from a high value, k h near the surface, to a low value k l at the Moho (numerical values are given in Table 1).Then we also assume, through Eq. ( 13), that conductivity depends on a cracking temperature T crack , for which intermediate (700 • C) and high (1000 • C) values have been tested.The resulting thermal conductivity that is used by the numerical model solving the temperature equation is obtained through Eq. ( 14), combining Eqs. ( 12) and ( 13).
The dashed and solid lines of Fig. 2 (middle panel) displays the variation of thermal conductivity obtained at the surface level (y = H ) according to the values assumed for the cracking temperature.The two series of cases corresponding to both cracking temperature assumptions will be shown in the following.

Effects of the accretion on the thermal and dynamic states of the ridge
Three series of cases have been computed to illustrate the potential effects of the crustal accretion mode.The first is a gabbro glacier structure (G); the second is a mixed structure with two lenses below the sheeted dyke and above the Moho (M); and the third is a sheeted sill structure (S) with superimposed sills delivering melt at the ridge axis.The G hypothesis consists of a melt intrusion through a shallow lens located just below the sheeted dyke (4.5 km above the MTZ).The M structure assumes two shallow and deep lenses respectively located just below the sheeted dyke and a few hundreds of meter above the MTZ (0.3 km and 4.5 km above the MTZ).Finally, the melt delivered for the S structure comes through nine sills, evenly stacked at the ridge axis, every 0.45 km, above the MTZ.
Figure 3 displays the temperature patterns (color palette) and the stream functions isocontours (black lines) obtained with an intermediate cracking temperature (T crack = 700 • C) for the G, M and S crustal accretion modes respectively in the top, middle and bottom panels.The contours of the streamfunctions reveal the trajectories resulting from the different accretion scenarios.As expected, the G structure (Fig. 3, top) induces, near the ridge, predominantly descending gabbro motion, while, the M motion is partly descending from  the upper lens and partly rising from the MTZ lens (Fig. 3, middle); and, for the S structure the numerous superimposed intrusive sills induce nearly horizontal motion (Fig. 3, bottom).In all the cases, far from the ridge, the stream functions converge toward laminar behaviors with cancelling of the vertical velocity.The top panel of Fig. 3 also depicts the temperature field obtained for the G crustal accretion mode (color scale).As they are sensitive to the heat advection and viscosity through crystallization, the temperature patterns are also dependent on the crustal accretion mode near the ridge.These behaviors result in variations that may locally reach several tens of degrees, but may remain sufficiently low to be difficult to be evaluated by the naked eye in most of the solution.To foil this inconvenience, the temperature fields obtained for the M and S structures have been represented as differences with one of the G hypothesis (two lower panels of Fig. 3).With the G accretion structure context, all the melt necessary to build the entire upper crust crosses the upper lens, carrying a lot of heat at a shallow level, where it is efficiently extracted by the cold thermal shallow gradient and enhanced hydrothermal cooling.This is no more the case for the M and S crustal accretion modes.Less heat is injected at shallow lens since melt is shared with the MTZ lens or with superimposed sills at the ridge axis.Compared to the thermal pattern obtained for the G structure, cold anomalies appear below the upper lens where, compared to the G structure temperature field, there is now a deficit of heat, and hot anomalies appear near the ridge in the deeper part of the solution where more heat is now advected.Far from the ridge, the influence of the melt intrusion patterns on the temperature damp rapidly converge toward conductive temperature profiles.
A second series of cases has been calculated with a higher gabbro cracking temperature (1000 • C) in order to illustrate the impact of the increase of the depth of hydrothermal cooling penetration in the crust.With this higher cracking temperature, the efficiency of hydrothermal cooling increases and penetrates deeper that induces colder temperatures than in the previous cases (Fig. 4, top panel).However, the evolution remains small for the cases presented I this study.They are more decipherable on the medium and lower panel of Fig. 4, where the lower temperature are betrayed by slightly stronger stream-function isocontours shapes, characteristic of a stronger viscosity due to the slightly colder temperature.From a dynamic point of view, the results, displayed on Fig. 4 remain strongly dependent on the G, M or S crustal accretion modes and similar to the ones previously obtained with an intermediate cracking temperature.The near-ridge gabbro motion is predominantly descending for the G structure (Fig. 4, top), is partly descending and partly rising for the M structure (Fig. 4, middle), and nearly horizontal for the S structure (Fig. 4, bottom).
As a common feature observed with both values of the cracking temperature, it appears that the global temperature pattern implies steeper thermal gradient for the G accretion mode than when the heat is shared with the MTZ lens or with superimposed sills at the ridge axis.This result is in qualitative agreement with the observations of Chen (2001), recalled in the introduction, since, in our case, at similar hydrothermal cooling, hot anomalies of approximately 100 • C appears in the lower crust with the M and S crustal accretion modes compared to the G accretion mode.
Figures 3 and 4 give x-y, Eulerian representation of steady state temperature patterns and motions reached at the end of the computing processes.From these solutions it is possible to calculate Lagrangian representations of the thermal   histories of tracers, following the T -t-x-y (temperature, time, offset, depth) trajectories of cooling gabbros in the lower crust.The panels in the left column of Fig. 5 give a representation of the thermal history of tracers during their travel from the intrusion at ridge level to their final emplacement in the cooled lower crust.During their transfers, the depth of tracers evolve with offset (and therefore time) according to the values of the vertical velocity, which near the ridge axis is strongly dependent on the crustal accretion mode.Therefore the vertical axis of Fig. 5 (left column) does not represent the depth of tracers, but the final height that are reached above MTZ at the lateral boundary of the computation grid.The left panels of Fig. 5 therefore displays detailed cooling histories of gabbros as they could be recorded in cooled crustal sections far from the ridge axis.The panels corresponding to the G and S crustal accretion modes (left column of Fig. 5) reveal monotonic increases with depth of the times spent at high temperature by tracers.This is no longer the case with the M crustal accretion mode, for which different marked thermal histories are obtained at depths located between the melt intrusions of the upper and lower lenses.The convergence of advected heat from the upper and MTZ lenses along with the slowness of the tracers near the ridge increases the times spent at high temperature for the tracers crossing at that crustal level.The slowness of the tracers mo-tions are portrayed in Figs. 3 and 4 (middle panels) by the smoothness of the stream-functions that denote low spatial derivatives and therefore, according to Eq. ( 2), slow velocity components.This result explains the particular shape of the thermal histories (Fig. 5, left and middle columns) for the M crustal accretion mode.Obviously, these thermal and dynamic effects due to merging of sills cannot exist for the G crustal accretion mode but are also present, for the S accretion structure, at the locations where the streams of neighbor sills merge (Fig. 5, left bottom panel).
The cooling histories of tracers are also shown in the middle and right columns of Fig. 5 that depict the T -t evolutions and the instantaneous cooling rates of tracers for each of the G, M and S melt intrusion geometries.As the depth of tracers varies during their trajectory, the tracers are reported at different final depth and final height that is reached above MTZ at the lateral boundary of the computation grid.Hence, they portray the cooling and cooling rate evolution of gabbros since their intrusion at the ridge axis until its final emplacement in cooled oceanic crustal section.For the G structure, the relative positions of the thermal history curves results are consistent with the monotonic evolution of thermal history described in the global thermal history (Fig. 5, left column, top panel).During the time evolution, the curves are regularly superimposed in such an order than the curves corresponding to the closest locations near MTZ at their final emplacement in the cooled crust display the hottest temperature all along the tracer trajectories.This T -t evolution of tracers is rather similar for the S structure, but different to that of the M structure, wherein inversions of the above curves order occur and (tracers at MTZ height 1818 m, red solid line; and at 2424 m, yellow solid line) several significant times are spent by tracers of shallower depth at higher temperatures than for deeper tracers (tracers at MTZ height equal to 606 m, solid blue line and 1212 m, green solid line).The instantaneous cooling rates (panels in the right column of Fig. 5) vary along the T -t flow paths for all accretion models.As longer time corresponds to farther distance from the on-axis ridge intrusion, they show that, from a Eulerian point of view, instantaneous cooling rates vary not only according to the intrusion mode but also as a function of the distance from the ridge axis and temperature.For most cases, the instantaneous cooling rates decreases with temperature and distance from the ridge axis.The slowest cooling rates take place generally at temperatures above or near the solidus and near the axis (i.e., shorter times).From these results, it is expected that natural proxies of igneous cooling rates (e.g.crystal size; Garrido et al., 2001) will differ in extent and variability from those based on intracrystalline diffusion (e.g., different minerals and diffusing species with diverse diffusion velocities).A comparison with these natural proxies of cooling rate would require, however, simulation of intracrystalline diffusion and crystallization along the T -t-x-y trajectories provided by our thermomechanical model.To better illustrate for the variability of gabbro cooling rates along the trajectories of different tracers as a function of the final distance above the MTZ, we have calculated averaged cooling rates following Eq.( 15) between two temperatures.In Eq. ( 15), dt is the time interval during which temperature ranges between T h and T l , the high and low temperatures encountered during the tracer travels.

CR =
Two arbitrary low and high temperature ranges have been chosen to illustrate the variability of cooling rates as a function of temperature interval and cracking temperature.The high temperature range starts at T h = 1275 • C, just below our melt intrusion temperature and ends with a low temperature T l = 1125 • C, thus covering most of our temperature crystallization range of the lower oceanic crust.The second corresponds to subsolidus temperature ranges, starting with T h = 1050, and ending with T l = 850 • C. With these temperature range choices, the average cooling rates record the cooling properties in different places of the model (those for which the temperature is actually ranging between, the high and low boundaries).According to the main locations of the 1275-1125 and 1050-850 • C isotherms in Figs. 3 and 4, the first will be more sensitive to the thermal structure near the ridge axis while the second will mainly record the thermal structures a few kilometers off-axis.
The left panel of Fig. 6 presents the average cooling rate profiles obtained for the three series of cases within the crystallization temperature range (while the right panel shows the result for the sub-solidus temperature range).Red, green and blue curves correspond respectively to the G, M and S crustal accretion modes.The results obtained with viscosity contrasts of two orders of magnitude and cracking temperatures T crac = 1000 • C have been drawn using heavy lines; the viscosity contrast of three orders of magnitude with cross symbols; and the viscosity contrast of two orders of magnitude but T crac = 700 • C with dashed lines.The comparisons of the relative location of the solid and cross-shaped symbols   in the various curves of Fig. 6 confirm that increases of two to three orders of magnitude viscosity contrasts have almost no effect on the results.In all the cases, the cross symbols superimpose almost perfectly with the corresponding heavy curves.The comparisons of solid line curves with dashed lines provide visualizations of the effects of the cracking temperature level.In the previous section, comparing Figs. 3 to 4, we emphasized that the final effects of the cracking temperature on the global thermal patterns of solutions was remaining of a few tens of degrees.Nevertheless, the slight deviations of trajectories and distortions of temperature due to these changes modify the average cooling rates of Fig. 6.The slightly colder environments, induced by the enhancement of the deep, near ridge cooling with the high cracking temperature contexts expose the gabbro to lower temperatures that results in higher cooling rates.However, these effects are weak and the slopes of the average cooling rates curves with depth remain almost unchanged.
For the average cooling rate profiles obtained with the lower temperature range, the effects of variable cracking temperature are similar but less pronounced (Fig. 6, right panel).In these cases the sampled areas of the model are far from the ridge axis, where the motion become laminar and the vertical temperature evolution tend toward conductive profile.It results decreases of the cooling rate differences versus accretion mode and gathering of the curves depicting the average cooling rates obtained for the different crustal accretion modes.The differences among the G, M and S accretion models are better discriminated by the average cooling rate obtained with the higher temperature range (Fig. 6, left).
As result, which it might been expected from instantaneous cooling rate evolution (Fig. 5), the profiles of integrated cooling rates with distance from the MTZ are monotonic for the G structure (red curves), display bi-modal shapes with marked minimum values at the levels where the flows from the upper and lower lenses merge for the M structure, and present saw tooth-like shapes for the S structures where the sills are merged.The same analysis is more difficult to apply to the average cooling rate profiles obtained with lower temperature range (Fig. 6, right).Indeed, the shift of the sampled areas far from the ridge and the lower value of the low temperature prevent the calculation of the average cooling rates in the lower part of the lower crust.It is clear from these results that the ability of average cooling rate to differentiate the various accretion scenarios will be better for the higher temperature ranges that actually correspond in the models to locations where the temperature and motions are the most influenced by the crustal accretion mode.

Summary and discussion
Our thermo-mechanical model offers a tool to explore the effects of deep, near off-axis hydrothermal cooling, viscosity contrast variable and crustal accretion mode on the thermal and dynamic patterns of flow near fast-spreading mid-ocean ridges.The series of cases presented in this paper simulate gabbro glacier "G", mixed MTZ and shallow lenses "M" or superimposed sills "S" crustal accretion modes, with various viscosity contrasts and, through the effects of a cracking temperature, various depths of hydrothermal cooling.Other accretion modes may be explored with our thermo-mechanical model, however.The differences of thermal structures obtained for the G, M and S hypotheses induce minor differences in temperature with depth and distance off-axis, which make it difficult to use temperature (or geophysical proxies of temperature) to discriminate among different crustal accretion scenarios.All cases investigated in this paper are consistent with the temperature structure at the ridge axis derived from geophysical studies at the East Pacific Rise (Dunn et al., 2000;Singh et al., 2006) that show a 8-12 km wide magma chamber (T < 1150 • C) with steep isotherms near the ridge axis.However, our results indicate that combinations of near-ridge flow patterns with local temperature differences both depending on the crustal accretion mode induce significant differences in the cooling histories of lower crustal gabbros.These differences are portrayed variations of the instantaneous cooling rates with time (Fig. 5) and average cooling rate with distance from the MTZ (Fig. 6) of the lower crust, which can both be useful to discriminate among different crustal accretion scenarios.Depending on the temperature interval used to average the instantaneous cooling rate, the profiles of the cooling rate with distance from the MTZ are however more or less able to discriminate among different crustal accretion modes.Cooling rates obtained from P. Machetel and C. J. Garrido: Numerical model of crustal accretion averaging a higher temperature interval sample area of the mode that are closer to the ridge axis, where the differences of average cooling rates are more discriminant of the crustal accretion mode.Conversely, average cooling rates obtained with a lower temperature interval sample areas of the accretion models where the thermal state and dynamic of the lower crust converge, respectively, toward a vertical conductive temperature profile and laminar motions.Hence, they are less discriminant of the accretion model.
Cooling rates obtained from petrographic and/or mineral compositional data in crustal samples of ophiolite or active mid-ocean ridges are integrated of cooling rates over T -t interval, which values are intrinsic to the methodologies used to derive the cooling rates.Absolute quantitative cooling rates of the plutonic crust have been determined by thermochronology (John et al., 2004), Crystal Size Distribution (CSD) of plagioclase in plutonic rocks (Marsh, 1988;Marsh, 1998, Garrido et al., 2001), an elemental diffusion in minerals from the plutonic crust based on geospeedometry (Coogan et al., 2002;Coogan et al., 2007;VanTongeren et al., 2008).These different proxies of magmatic cooling rates record the cooling of oceanic gabbros in T -t interval between the liquidus and the solidus temperature (i.e., the crystallization time; CSD of Garrido et al., 2001) or elemental diffusion in minerals from the plutonic crust, as those based on geospeedometry (Coogan et al., 2002;Coogan et al., 2007;VanTongeren et al., 2008), record the cooling rate in the T -t interval where a characteristic exchange diffusion is effective.Garrido et al. (2001) measured CSD from plagioclase in the Khafifah section of the Wadi Tayin massif and found evidence of a transition from conduction-dominated cooling in the lower gabbros (below 1500 m above Moho) to hydrothermally dominated cooling in the upper gabbros (above 2500 m).They concluded their data were consistent with the S model of accretion.However, their cooling profiles did not show the same kind of evolution with depth than the average cooling rate presented in this numerical study.They were displaying an upper crust value 1.5 to 2 times faster than lower crust values.In light of the present numerical results, this could be compatible with the three accretion structure hypotheses.Coogan et al. (2002), using the Ca diffusion in olivine from a Wadi Abyad massif crustal section, reported that cooling rates decrease rapidly with depth by several orders of magnitude between the top and bottom of the lower crust.They also mention that the cooling depth profile matches that of conductive models.These authors concluded in favor of a crystallization occurring inside of the magmatic chamber and hence a G ridge structure.VanTongeren et al. (2008) extended the work of Coogan et al. (2002) in the Wadi Tayin massif of the Oman ophiolite.The two studies differ both in amplitudes and shapes of cooling rates profiles (see the comparison in Fig. 7 of VanTongeren et al., 2008 andFig. 3 of MacLennan et al., 2005).VanTongeren et al. (2008) argued that these differences reflect distinct ther-mal histories due to differences in crustal thickness and/or the geodynamic setting.However, the cooling rates recalculated by VanTongeren et al. ( 2008) using Coogan's data (Coogan et al., 2002) remain several orders of magnitude faster than the ones calculated by the former.Such differences between the results of VanTongeren (2008) and those of Cogan et al. (2002) suggest, as recalled by Coogan et al. (2002) and MacLennan et al. (2005), that large uncertainties in petrology may probably come from a deficit of constraints on the values of diffusion parameters.
The present thermomechanical model provides much more detail of cooling rate history of the oceanic crust through the instantaneous cooling rates than that obtained from natural proxies.A strict comparison with petrological derived cooling rates would require simulation of crystallization and chemical diffusion along the T -t-x-y trajectories employing numerical models of net-transfer and exchange reactions in combination with estimates of intracrystalline diffusion.Such simulation is beyond the scope of the numerical model, which, however, lays the foundation for this development.The results of our model indicate, however, that some assumptions often made using petrological-derived cooling rates to discriminate between accretion models are simplistic.For instance, our results show that cooling rates at supersolidus conditions are generally slower than those subsolidus conditions.It is hence likely that natural proxies of cooling rates at super solidus conditions will provide different values than those using proxies based on subsolidus intracrystalline diffusion.The present models show that monotonic variations of the cooling rates with depth are not necessarily symptomatic of conductive cooling or a G crustal accretion structure (Coogan et al. 2002), as this variation may be also produced by an S accretion geometry.
In spite and because of these uncertainties, our results suggest that numerical modeling of crustal accretion modes and their consequences in terms of instantaneous and average cooling rates may provide efficient tools to try to discriminate between different crustal accretion modes at fastspreading mid-ocean ridges.Further use of the present thermomechanical model to discriminate between crustal models would require benchmarking the results with geophysical observables.

Figure 1 :
Figure 1: Sketch of the stream-function (top) and temperature (bottom) internal and boundary conditions.The computation grid is composed of Nx vertical columns (from 1 to Nx, x = -L to x =L) and Ny rows (from 1 to Ny, y = 0 to y = H).The ridge axis is located between the points Nx/2 and (Nx/2)+1 where the amplitude of the stream-function jump is equal to c = 2 Vp H, the total flux of crust that leaves the computation box through the left and right lateral boundaries.At the top of the crust and at the Moho level, is constant (impervious boundary), except at the ridge axis.Its bottom left value has been set to bl=Vp H. Internal stream function jumps,  on the two central columns, drive the melt intrusion through sills and/or lenses.Initial and internal boundary conditions are also applied to the temperature field according to a half-space lithospheric cooling law(Eq.14).

Fig. 1 .
Fig. 1.Sketch of the stream-function (top) and temperature (bottom) internal and boundary conditions.The computation grid is composed of N x vertical columns (from 1 to N x , x = −L to x = L) and N y rows (from 1 to N y , y = 0 to y = H ). The ridge axis is located between the points N x /2 and (N x /2) + 1 where the amplitude of the stream-function jump is equal to c = 2V p H , the total flux of crust that leaves the computation box through the left and right lateral boundaries.At the top of the crust and at the Moho level, is constant (impervious boundary), except at the ridge axis.Its bottom left value has been set to bl = V p H . Internal stream function jumps, δ , on the two central columns, drive the melt intrusion through sills and/or lenses.Initial and internal boundary conditions are also applied to the temperature field according to a half-space lithospheric cooling law (Eq.14).

Figure 3 :
Figure 3: Temperature and stream function obtained for the G (top), M (middle) and S (bottom) crustal accretion modes.To increase the readability of the figure the temperature field obtained for the G structure has been represented in the top panel.Values of temperature correspond to the top color scale.For the M and S structures (medium and lower panels, the temperatures are represented as difference with the gabbro glacier (G) structure above.The direction of velocity is directly given by the streamfunction contours (black lines).White lines display the locations of the 1125, 1050 and 850°C isotherms that will be used for the computation of the average cooling rates (see text and Fig.6).The results have been obtained with an intermediate hydrothermal cracking temperature Tcrac = 700°C and a viscosity contrast of two orders of magnitudes.

Fig. 3 .
Fig. 3. Temperature and stream function obtained for the G (top), M (middle) and S (bottom) crustal accretion modes.To increase the readability of the figures the temperature field obtained for the G structure has been represented in the top panel.Values of temperature correspond to the top color scale.For the M and S structures (medium and lower panels, the temperatures are represented as difference with the gabbro glacier (G) structure above.The direction of velocity is directly given by the stream-function contours (black lines).White lines display the locations of the 1125, 1050 and 850 • C isotherms that will be used for the computation of the average cooling rates (see text and Fig.6).The results have been obtained with an intermediate hydrothermal cracking temperature T crac = 700 • C and a viscosity contrast of two orders of magnitudes. 31

Fig. 5 .
Fig. 5. Thermal history of the gabbro versus their final height above MTZ in the cooled lower crust and the G (top panels), M (middle) and S (bottom) crustal accretion modes.Left panels display the T (t) evolution of Lagrangian tracers along their trajectories in the lower crust; middle panels give this temperature evolution for different depths (given in height above MTZ in panels); evolutions of the instantaneous cooling rates are given in the right panels.The results have been obtained with T crac = 700 • C. 33

Figure 6 :
Figure 6: Average cooling rates (ACR) calculated from Eq. 15 using two temperature intervals.The first, from 1275 to 1125°C, covers most of our crystallization range (left panel) the second, from 1050 to 850°C involves mainly sub-solidus temperatures (right panel).The vertical coordinate corresponds to the final depths reached by the tracers at the end of computation.ACR curves are drawn for the G (red curves), M (green curves) and S (blue curves) crustal accretion modes.In each panel, heavy solid lines correspond to cases with two orders of magnitude viscosity contrasts and high cracking temperature (1000 °C).The cases obtained by changing the viscosity contrasts superimpose perfectly (X symbols).Finally, the results obtained with an intermediate hydrothermal cracking temperature of 700 °C are displayed with dashed line.

Fig. 6 .
Fig. 6.Average cooling rates (ACR) calculated from Eq. 15 using two temperature intervals.The first, from 1275 to 1125 • C, covers most of our crystallization range (left panel) the second, from 1050 to 850 • C involves mainly sub-solidus temperatures (right panel).The vertical coordinate corresponds to the final depths reached by the tracers at the end of computation.ACR curves are drawn for the G (red curves), M (green curves) and S (blue curves) crustal accretion modes.In each panel, heavy solid lines correspond to cases with two orders of magnitude viscosity contrasts and high cracking temperature (1000 • C).The cases obtained by changing the viscosity contrasts superimpose perfectly (X symbols).Finally, the results obtained with an intermediate hydrothermal cracking temperature of 700 • C are displayed with dashed line.