Description: The genesis of oceanic crust at intermediate to fast spreading ridges occurs by the crystallization of mantle melts accumulated in at least one shallow melt lens situated below the ridge axis. Seismic reflection data suggest that the depth of this melt lens is inversely correlated with spreading rate and thereby magma supply. The heat released in it by crystallization and melt injection is removed by a combination of hydrothermal cooling and diffusion. Due to the different time scales of hydrothermal cooling and crustal accretion, numerical models have so far focused on only one of the two processes. Here we present the results from a coupled model that solves simultaneously for crustal accretion and hydrothermal cooling. Our approach resolves both processes within one 2D finite-element model that self-consistently solves for crustal, mantle, and hydrothermal flow. The formation of new oceanic crust is approximated as a gabbro glacier, in which the entire lower crust crystallizes in one shallow melt lens. We find that the depth of the melt lens and the shape of hot (potentially molten) lower crust are highly dependent on the ridge permeability structure. The predicted depth of the melt lens is primarily controlled by the permeability at the ridge axis, whereas the off-axis permeability determines the width of hot lower crust. A detailed comparison of the modeling results with observed locations of the melt lens at intermediate to fast spreading ridges shows that only a relatively narrow range of crustal permeabilities is consistent with observations. In addition, we find significant deviations between models that resolve or parameterize hydrothermal cooling: the predicted crustal thermal structures show major differences for models that predict the same melt lens location. This illustrates the importance of resolving hydrothermal flow in simulations of crustal accretion.

Description: Hydration of the oceanic lithosphere is an important and ubiquitous process which alters both the chemical and physical properties of the affected lithologies. One of the most important reactions that affect the mantle is serpentinization. The process of serpentinization results in a drastic decrease in the density (up to 40%), seismic velocity and brittle strength as well as water uptake of up to 13 wt.% of the ultramafic rock. In this paper, we use numerical models to study the amount and extent of serpentinization that may occur at mid-ocean ridges and its effects on fluid flow within the lithosphere. The two dimensional, FEM model solves three coupled, time-dependent equations: (i) mass-conserving Darcy flow equation, (ii) energy conserving heat transport equation and (iii) serpentinization rate of olivine with feedbacks to temperature (exothermic reaction), fluid consumption and variations in porosity and permeability (volume changes). The thermal structure of the ridge is strongly influenced by rock permeability in addition to the spreading velocity of the ridge. Increased rock permeability enhances hydrothermal convection and results in efficient heat mining from the lithosphere whereas higher spreading velocities result in a higher thermal gradient. Serpentinization of the oceanic mantle, in turn, depends on the aforementioned, competing processes. However, serpentinization of mantle rocks is itself likely to result in strong variations of rock porosity and permeability. Here we explore the coupled feedbacks. Increasing rates of serpentinization lead to large volume changes and therefore, rock fracturing thereby increasing rock porosity/permeability while as serpentinization reaches completion, the open pore space in the rock is reduced due to the relative dominance of mineral precipitation. Although, variations in the relation between porosity and permeability and serpentinization before the reaction reaches completion do not significantly affect the degree of serpentinization, we find that unreasonably large portions of the mantle would be serpentinized if rock closure does not occur at the final reaction stage. The amount of water trapped as hydrous phases within the mantle shows a strong dependency on the spreading velocity of the ridge with water content ranging from 0.18 × 105 kg/m2 to 2.52 × 105 kg/m2. Additionally, two distinct trends are observed where the water content in the mantle at slow-spreading ridges drops dramatically with an increase in spreading velocity. The amount of water trapped in the mantle at fast-spreading ridges, on the other hand, is lower and does not significantly depend on spreading velocity.

Description: A model for the release of Li, Be and B from progressively dehydrating altered oceanic crust during subduction is presented. Combining clinopyroxene/fluid partition coefficients determined experimentally in an earlier study Brenan et al. [Brenan, J.M., Ryerson, F.J., Shaw, H.F., 1998. The role of aqueous fluids in the slab-to-mantle transfer of boron, beryllium, and lithium during subduction: Experiments and models. Geochim. Cosmochim. Acta 62, 3337–3347] with apparent mineral/clinopyroxene partition coefficients as observed in natural high-pressure metamorphic rocks Marschall et al. [Marschall, H.R., Altherr, R., Ludwig, T., Kalt, A., Gméling, K., Kasztovszky, Zs., 2006a. Partitioning and budget of Li, Be and B in high-pressure metamorphic rocks. Geochim. Cosmochim. Acta 70, 4750–4769] results in a set of mineral/fluid partition coefficients for high-pressure metamorphic minerals. Mineral modes of altered oceanic crust as a function of pressure and temperature along a given subduction path can be derived from thermodynamic calculations using the program PerpleX. Combination of these modes with mineral/fluid partition coefficients results in whole rock/fluid partition coefficients at any stage of the P–T path including information on the amount of fluid released at any depth. Based on these data, the concentrations of Li, Be and B in subducting rocks and released fluids along a given P–T path can be modelled. The derived information on B concentrations in rocks and fluids are combined with the temperature-dependent fractionation of B isotopes in order to model the B isotopic evolution of subducting rocks and released fluids. Model calculations are performed for two slightly different chemical compositions (hydrous MORB without K and with 0.5 wt.% K2O), in order to demonstrate the impact of phengite on the boron budget. Provided the necessary input data are available, the concept of such a model could be employed to quantify the trace element release from the slab from any lithology along any reasonable P–T path.

Description: We present 2D and 3D numerical model calculations that focus on the physics of compositionally buoyant diapirs rising within a mantle wedge corner flow. Compositional buoyancy is assumed to arise from slab dehydration during which water-rich volatiles enter the mantle wedge and form a wet, less dense boundary layer on top of the slab. Slab dehydration is prescribed to occur in the 80–180 km deep slab interval, and the water transport is treated as a diffusion-like process. In this study, the mantle's rheology is modeled as being isoviscous for the benefit of easier-to-interpret feedbacks between water migration and buoyant viscous flow of the mantle. We use a simple subduction geometry that does not change during the numerical calculation. In a large set of 2D calculations we have identified that five different flow regimes can form, in which the position, number, and formation time of the diapirs vary as a function of four parameters: subduction angle, subduction rate, water diffusivity (mobility), and mantle viscosity. Using the same numerical method and numerical resolution we also conducted a suite of 3D calculations for 16 selected parameter combinations. Comparing the 2D and 3D results for the same model parameters reveals that the 2D models can only give limited insights into the inherently 3D problem of mantle wedge diapirism. While often correctly predicting the position and onset time of the first diapir(s), the 2D models fail to capture the dynamics of diapir ascent as well as the formation of secondary diapirs that result from boundary layer perturbations caused by previous diapirs. Of greatest importance for physically correct results is the numerical resolution in the region where diapirs nucleate, which must be high enough to accurately capture the growth of the thin wet boundary layer on top of the slab and, subsequently, the formation, morphology, and ascent of diapirs. Here 2D models can be very useful to quantify the required resolution, which we find for a 1019 Pa · s mantle wedge to be about 1 km node spacing for quadratic-order velocity elements.

Description: This study provides new estimates for the global offshore methane hydrate inventory formed due to microbial CH4 production under Quaternary and Holocene boundary conditions. A multi-1D model for particular organic carbon (POC) degradation, gas hydrate formation and dissolution is presented. The novel reaction-transport model contains an open three-phase system of two solid compounds (organic carbon, gas hydrates), three dissolved species (methane, sulfates, inorganic carbon) and one gaseous phase (free methane). The model computes time-resolved concentration profiles for all compounds by accounting for chemical reactions as well as diffusive and advective transport processes. The reaction module builds upon a new kinetic model of POC degradation which considers a down-core decrease in reactivity of organic matter. Various chemical reactions such as organic carbon decay, anaerobic oxidation of methane, methanogenesis, and sulfate reduction are resolved using appropriate kinetic rate laws and constants. Gas hydrates and free gas form if the concentration of dissolved methane exceeds the pressure, temperature, and salinity-dependent solubility limits of hydrates and/or free gas, with a rate given by kinetic parameters. Global input grids have been compiled from a variety of oceanographic, geological and geophysical data sets including a new parameterization of sedimentation rates in terms of water depth. We find prominent gas hydrate provinces offshore Central America where sediments are rich in organic carbon and in the Arctic Ocean where low bottom water temperatures stabilize methane hydrates. The world’s total gas hydrate inventory is estimated at 0.82 x 10sup13 m3 - 2.10 x 10sup15 m3 CH4 (at STP conditions) or, equivalently, 4.18–995 Gt of methane carbon. The first value refers to present day conditions estimated using the relatively low Holocene sedimentation rates; the second value corresponds to a scenario of higher Quaternary sedimentation rates along continental margins. Our results clearly show that in-situ POC degradation is at present not an efficient hydrate forming process. Significant hydrate deposits in marine settings are more likely to have formed at times of higher sedimentation during the Quaternary or as a consequence of upward fluid transport at continental margins.

Description: The surface of the solid Earth is effectively stress free in its subaerial portions, and hydrostatic beneath the oceans. Unfortunately, this type of boundary condition is difficult to treat computationally, and for computational convenience, numerical models have often used simpler approximations that do not involve a normal stress-loaded, shear-stress free top surface that is free to move. Viscous flow models with a computational free surface typically confront stability problems when the time step is bigger than the viscous relaxation time. The small time step required for stability (〈2. Kyr) makes this type of model computationally intensive, so there remains a need to develop strategies that mitigate the stability problem by making larger (at least ~10 Kyr) time steps stable and accurate. Here we present a new free-surface stabilization algorithm for finite element codes which solves the stability problem by adding to the Stokes formulation an intrinsic penalization term equivalent to a portion of the future load at the surface nodes. Our algorithm is straightforward to implement and can be used with both Eulerian or Lagrangian grids. It includes α and β parameters to respectively control both the vertical and the horizontal slope-dependent penalization terms, and uses Uzawa-like iterations to solve the resulting system at a cost comparable to a non-stress free surface formulation. Four tests were carried out in order to study the accuracy and the stability of the algorithm: (1) a decaying first-order sinusoidal topography test, (2) a decaying high-order sinusoidal topography test, (3) a Rayleigh-Taylor instability test, and (4) a steep-slope test. For these tests, we investigate which α and β parameters give the best results in terms of both accuracy and stability. We also compare the accuracy and the stability of our algorithm with a similar implicit approach recently developed by Kaus et al. (2010). We find that our algorithm is slightly more accurate and stable for steep slopes, and also conclude that, for longer time steps, the optimal α controlling factor for both approaches is ~2/3, instead of the 1/2 Crank-Nicolson parameter inferred from a linearized accuracy analysis. This more-implicit value coincides with the velocity factor for a Galerkin time discretization applied to our penalization term using linear shape functions in time.

Description: Most simple models for cooling of sheet like intrusions, the Hawaiian lava lakes being one example, neglect the effect of side wall heat loss on the overall thermal evolution. In this paper we extend a conventional one-dimensional (1D) model for cooling of sheet like intrusions to account for lateral heat loss by either prescribing a fixed side wall heat flux, or a heat flux controlled by heat transport in the surrounding wall rock. In the first part of the study we analyze the general interplay between side wall cooling and the thermal evolution of the system; the second part focuses on a comparison between our modeling results including models without and with lateral heat flux, and the Hawaiian lava lake data. This comparison leads to the following three main conclusions: (1) Side wall cooling does have a significant impact on the cooling history of lava lakes. (2) Models assuming a time dependent temperature profile in the wall rock lead to a better fit with the measured temperature data. (3) Due to the sluggish conductive heat transfer in the mush its thermal evolution is significantly decoupled from the temperature evolution of the convecting bulk liquid.

Description: This study explores a chemo-thermo-dynamic subduction zone model that solves for slab dehydration during subduction. We investigate how changes in the incoming plate's hydration and thermal structure may effect the efficiency of sub-arc water release from sediments, crust, and serpentinized mantle. We find that serpentinized lithospheric mantle may not only be an important fluid source to trigger arc melting but is also an efficient ‘transport-lithology’ to recycle chemically bound water into the deeper mantle. In fact, an old slab may remain sufficiently cold during subduction to retain up to 40% of its initial ‘mantle’ water at 8 GPa (∼240-km depth) after serpentine transforms to higher pressure hydrous phase A. Furthermore, deep water recycling at subduction zones is parameterized in terms of slab age and speed. Coupling this parameterization to a parameterized mantle convection evolution model allows us to calculate the mantle-surface geologic water cycle throughout the Earth's history. We find that the present-day Earth mantle may be highly outgassed containing only a small fraction of the Earth's water, which would mostly be recycled water from the exosphere.

Description: Highlights • Low- and medium pressure TTGs are formed continuously. • High-pressure TTGs are mainly formed during episodic lithospheric recycling events. • The ratio of low- to medium-pressure TTGs corresponds to time-dependent crustal thickness. • Formation of high-pressure TTGs requires a delamination or subduction type process. The appearance of the earliest felsic crust can be estimated by dating zircons and rocks of tonalite-trondhjemite-granodiorite (TTG) composition. However, the necessary geodynamic processes that form the basis for metamorphism and differentiation as well as the role of emerging TTG crust on evolving crustal dynamics is still poorly understood. To investigate the formation of felsic crust with TTG composition, we conduct a detailed analysis of a series of previously published 3D high-resolution magmatic-thermomechanical models at elevated mantle temperature corresponding to Archean conditions. In these models we observed two distinct phases during coupled cyclic tectonomagmatic crust-mantle evolution: a long quiet growth phase followed by a short rapid overturn phase. Results of the detailed model analysis presented here suggest that (1) Low- and medium-pressure TTGs are formed at the base of the crust during both growth and overturn phase. The formation of low- and medium-pressure TTGs is linked with Moho depth. The ratio of low- to medium-pressure TTGs changes with crustal growth or thinning and gives an approximation for crustal thickness. (2) To form high-pressure TTGs an entirely different mechanism is required, as hydrated basaltic rocks need to be buried below the crust. Direct partial melting of cold eclogitic drips can be excluded as a valid mechanism due to their low temperatures and rapid sinking into the deep mantle. Rather we suggest delamination (peeling-off) or subduction as the main process for some high-pressure TTG production.

Description: Highlights • Sedimentation-driven gas hydrate recycling is cyclic in nature with time scales set by reactive multi-phase transport. • Each cycle can be divided into three distinct phases: 1) gas accumulation phase, 2) gas breakthrough phase and 3) uninhibited hydrate build-up phase. • In the presence of sufficient accumulated gas, convex deposition of hydrate acts like a mechanical nozzle for the ascending gas flow. Gas hydrate recycling is an important process in natural hydrate systems worldwide and frequently leads to the high gas hydrate saturations found close to the base of the gas hydrate stability zone (GHSZ). However, to date it remains enigmatic how, and under which conditions, free gas invades back into the GHSZ. Here we use a 1D compositional multi-phase flow model that accounts for sedimentation to investigate the dominant mechanisms that control free gas flow into the GHSZ using a wide-range of parameters i.e. hydrate formation kinetics, sediment permeability, and capillary pressure. In the first part of this study, we investigate free gas invasion into the GHSZ without any sedimentation, and analyse the dynamics of hydrate formation in the vicinity of the base of GHSZ. This helps establish plausible initial conditions for the main part of the study, namely, hydrate recycling due to rapid and continuous sedimentation. For the case study, we apply our numerical model to the Green Canyon Site 955 in the Gulf of Mexico, where the reported high hydrate saturations are likely a result of hydrate recycling driven by rapid sedimentation. In the model, an initial hydrate layer forms due to the invasion of a specified volume of rising free gas. This hydrate layer is consistent with the local pressure, temperature and salinity state. This hydrate layer is then thermally de-stabilised by sedimentation resulting in free gas formation and hydrate recycling. A key finding of our study is that gas hydrate recycling is a cyclic process which can be divided into three phases of 1) gas hydrate melting and free gas nozzling through the hydrate layer, 2) formation of a new gas hydrate layer as the old layer vanishes, and 3) fast uninhibited grow of a new hydrate layer. High hydrate saturations of about 80% can be attained purely through physical, burial-driven recycling of gas hydrates, without any additional gas input from other sources. Hydrate recycling is, therefore, highly dynamic with its own inherent cyclicity rather than a gradual process paced by the rate of sediment deposition.