2D reactive transport simulations of mid-ocean ridge hydrothermal systems

Water-rock interactions in mid-ocean ridge hydrothermal systems are a critical part of Earth system evolution. Extensive insights have been developed from vent fluid chemistry and laboratory experiments, but these leave unanswered many questions about the temporal evolution and spatial structure of the hydrothermal systems that can only be addressed with reactive transport simulations. Other issues are the effects of changing spreading rates and seawater chemistry through Earth history. We are addressing this problem using the Toughreact code, starting with 2D static (no seafloor spreading) simulations of the near-axis region where most of the interaction occurs. The simulations use a dual-permeability grid to represent fractured rocks, and also have a formulation for Sr isotope exchange. Vent fluid Ca, Mg, SO4, and Na concentrations and Sr isotopes can be used as a guide to fluid chemical evolution. Initial simulations reproduce modern vent fluid chemistry even with maximum temperature only at 380°C, and suggest that fluids need not be in equilibrium with the rocks at any point in the system. Model fluids continue to evolve chemically even in the upflow zone prior to venting. The effects of different seawater chemical composition, as proposed for the Cretaceous, for example, can be captured with charge-balance models.


Scientific background
Mid-ocean ridge (MOR) hydrothermal systems are a major component of Earth's geochemical cycles. The chemical exchange that occurs as seawater circulates through basalt and gabbro of varying temperature affects the balance of ions in seawater, and influences the carbon cycle through the coupling of Ca and Mg and carbonate sedimentation. Understanding the geologic history of the oceans, and the relationships between seawater chemistry and climate, requires that the chemical exchange between seawater and oceanic crust be understood not only in its present form, but as it might have been at different stages of Earth history extending back to earliest Precambrian times. For example, differences in the chemical composition of seawater in the geologic past could have had a significant effect on the seawater-basalt exchange in MOR hydrothermal systems, especially at times when seawater Mg, pH, and SO4 were lower, and Ca, Sr, and temperature were higher than at present [1]. There is also a need to formulate better models for the effects of changing seafloor spreading rates and the relative magnitudes and effects of high-temperature and low-temperature seawater-basalt exchange.

Model details
Our simulations use the parallel reactive-transport code TOUGHREACT V3.3-OMP [2] which includes kinetically controlled mineral-fluid interaction, mineral solid-solutions, aqueous kinetics, and Sr isotopic exchange, plus fluid flow and heat transfer. A dualpermeability grid is used to simulate flow through fractures and fracture-matrix chemical and heat exchange. The thermodynamic database is adapted from that of the University of Oregon, and kinetic formulations are first order for both dissolution and precipitation. The partitioning of Sr between fluids and secondary minerals (anhydrite, tremolite, epidote) is done by using Sr-bearing solid solution endmembers and adjusting the ∆Gf (or log Kf) so that an appropriate KSr is obtained (these are still poorly constrained by experimental data).
Mineral dissolution kinetics, which largely control the evolution of the system, are controlled through the reactive surface area of the primary minerals and are adjusted to yield overall reaction rates in the range of 10 -10 to 10 -9 mol/kg fluid/sec at 300-350C. This range is based on available constraints from natural systems [3] and also is limited by fracture-matrix exchange and the ability of the code to reach chemical convergence.
The grid representing ocean floor rocks is 1.2 km high by 1.6 km wide, each grid block is 20m x 20m, and set to represent a fractured medium with 2 meter fracture spacing everywhere, and uniform fracture permeability of 10 -14 m 2 and matrix permeability of 10 -18 m 2 . The matrix has 5% porosity and the fractures have 50% porosity. We experimented with other fracture spacings and permeabilities. The permeable rock part of the grid is overlain by 200 meters of ocean water that starts with a temperature of 4°C for simulations of modern systems, and 15°C for Cretaceous systems and with chemical composition, pH, and 87 Sr/ 86 Sr appropriate to each time period.
Flow is driven by heating from below with a specified heat flux, adjusted so that the maximum temperature does not exceed about 380°C because of the limitations of the thermodynamic database. For the permeability we use, the heat flux is 8 W/m 2 at the ridge axis, dropping to 2 W/m 2 at 1.6 km from the ridge with a Gaussian profile. The maximum heat flux is equivalent to having the top of a magma chamber at 1200°C situated about 200 meters below the bottom of our grid near the ridge axis.
The characteristics of the system present a challenge for the code. Temperature varies from 4°C to 380°C, so the rates of reactions vary by about 6 orders of magnitude. The system takes 5000+ years to reach steady state temperature. To achieve the results we have thus far, the system was started with a lower heating rate (maximum 2 W/m 2 ) and with reactions set to be extremely slow by adjusting the mineral reactive surface areas. This combination of parameters was run for about 6000 model years until thermal steady state was achieved. Heating rate was then increased to 6 W/m 2 , for 500 model years and then 8 W/m 2 for 500 model years. Reactive surface areas were then increased, gradually, first running the system for 500 years for each step increase in reactive surface area, and then 100 years each step until realistic mineral reaction rates were reached.

Fluid flow and temperature field
In our first simulations we used a fracture spacing of 10 meters, a maximum heat flux of 12 W/m 2 near the ridge axis, and fracture permeability 10 -13 m 2 in the vertical direction (parallel to the ridge) and 10 -14 m 2 horizontal. This produced reasonable fluid fluxes in that the upflow fluxes were about 2-4 x 10 -4 kg/m 2 /sec [4]. However, only a tiny region of the grid was heated to above 300°C, and the fluid spent so little time in the high temperature part of the system that there was little chemical evolution. Model vent fluid Mg was >45 mM and nowhere in the system was it much lower. Decreasing the vertical permeability to 5 x10 -14 m 2 did not change the result much. With a uniform fracture permeability of 10 -14 m 2 , the fluid chemistry evolves substantially and model vent fluids have zero Mg 2+ and SO4 2and elevated Ca 2+ . Because we wanted to evaluate controls on the fluid evolution, we used the uniform permeability model for further simulations. The resulting model temperature and fluid flow pattern (Figure 1) are similar in many respects to other flow models [4]. The upflow zone at the ridge axis has a maximum fluid flux of about 10 -4 kg/m 2 /sec, which is the right order of magnitude but a little lower than in other simulations of fluid flow [4]. The fraction of the grid where steady state temperature is above 300°C is very small (about 20 grid blocks out of 4800). The axial upflow zone is narrow (50-100 meters half-width) and there is a distinct transition to downflow just 250 m from the axis. A substantial fraction of the total downflow is in the region between 250 and 450 m from the axis. This near-axis downflow is enhanced if vertical permeability is increased, and its effect, because the cold downwelling fluid is entrained into the hot upflow, is to dilute and cool the upflowing evolved fluids. There is some recycling of reacted fluid from the upflow zone into the near-axis downflow region (Figure 2).
With the permeability we are using, fluid moves through the >300°C zone quickly (50-100 years), but long enough to react sufficiently with the rocks to strip Mg and SO4 from the fluid. The permeability we are using is difficult to compare with estimates based on core measurements [5], but overall is consistent with observations. If fracture spacing is 2 meters, the key fractures won't be represented in measurements of core samples. Our matrix porosity of 5% may be somewhat high, but this helps to involve the matrix in the chemical exchange more. If matrix porosity were lowered to 2-3%, it would need to be offset with smaller fracture spacing to get the fluid chemistry to evolve enough.
Another key insight is that the effective Damkohler number in the system doesn't vary quite as much as might be expected because of the correlation between fluid velocity and temperature. In the low temperature regions (T = 60-80°C), fluid fluxes are small (ca. 10 -6 kg/m 2 /sec). In the high T regions (T= 250-350°C), fluid fluxes are 100x greater. The contrast in reaction rate of plagioclase is a factor of 5000, so the Damkohler number is only about 50x higher at 300°C than at 75°C. Fluid in the upflow zone traverses the 1.2 km vertical distance in 200 years.  An interesting feature is that Mg increases slightly before it decreases, which might be a kinetic effect resulting from olivine dissolution being too fast for secondary mineral precipitation to keep up. Vent fluid has 28 mM Ca, similar to the modern vent fluid average. Seawater Ca is almost completely stripped from the fluid before it enters the higher temperature reactive zone where it then increases. The vent fluid Ca concentration anti-correlates with Na + concentration, and both depend on the efficiency of albite precipitation relative to amphibole. In this simulation, albite precipitation was enhanced relative to the rate expressions available in the literature.

Mineralogy and fluid chemistry
An important question is where and how much anhydrite precipitates. Figure 3 shows the fraction of fracture volume filled with anhydrite at the end of the simulation for modern seawater. Another point is that nowhere in the system does plagioclase ever reach equilibrium; it remains undersaturated by at least 10,000x everywhere even though vent fluid chemistry is similar to modern vents. Epidote forms only at the highest temperatures and in the upflow zone where fluid velocity is highest.
In changing the seawater chemistry to one characteristic of mid-Cretaceous times (35 mM Ca, 30 mM Mg, 10 mM SO4), we obtained results compatible with simple charge balance models [1]. Using Cretaceous seawater the vent fluid Ca is much higher (56 mM) and the amount of anhydrite formed is about half, relative to the modern system. Vent fluid Mg and SO4 are still zero. Major differences in Sr isotopes are expected but we have not yet completed those simulations.

Conclusions
Our initial models are simple but we can begin to understand what controls the results. The models will need other features to get a more robust idea of how MOR systems work. One significant issue is seafloor spreading, which can be reasonably fast in comparison to the chemical evolution of the system. At 3 cm/yr half-spreading rate, the grid would move 30 meters (1.5 grid blocks) every 1000 years. Since it takes significantly longer than 1000 years for the flow and thermal structure to reach steady state, the thermal effects of spreading can be significant. Spreading may produce a different thermal and flow structure, and of course this will depend on spreading rate. Dikes emplaced at the axis may affect heating. Also, because all of the rocks in the spreading system start at the ridge axis where temperature and fluid flow is highest, the entire grid probably should be composed of altered rocks. We believe we can account for most of these effects in future models, and also evaluate the effects of different permeability structure.