2.1 Hydrothermal circulation model

Hydrothermal flow circulating around midocean ridges is simulated in 2D based on conservations of energy, mass and momentum (Steefel and Lasaga, 1994; Cherkaoui et al., 2003; Iyer et al., 2010). Assuming a steady state, the mass conservation of fluid is represented by

where ∇ is the vector differential operator ($\mathrm{\nabla }=(\partial /\partial x,\partial /\partial y)$) and q is the mass flux vector ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). Conservation of momentum is realized by adopting Darcy's law (e.g., Steefel and Lasaga, 1994):

where k is the permeability of oceanic bulk rock (m²), ν is the kinematic viscosity of water (m² yr⁻¹), P is the fluid pressure (Pa), ρ_f is the fluid density (kg m⁻³), and g is the gravity vector given by $\mathit{g}=(\mathrm{0},-g)$ where g represents the acceleration by gravity (m yr⁻²; e.g., Steefel and Lasaga, 1994). The energy conservation is represented by the following:

where t is time (yr), ϕ is the porosity, ρ_m is the density of oceanic rock (kg m⁻³), $c_{\mathrm{p}}^{\mathrm{f}}$ and $c_{\mathrm{p}}^{\mathrm{m}}$ are the specific heat capacity of water and oceanic rock, respectively ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), T is the temperature (K), and κ is the thermal conductivity of oceanic rock ($\mathrm{J}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$). The thermodynamic and transport properties of water ($c_{\mathrm{p}}^{\mathrm{f}}$, ρ_f and ν) are obtained through a FORTRAN90 library STEAM, which is based on the NBS steam table package (Meyer et al., 1983; Haar et al., 1984). As in Iyer et al. (2010), we assume pure water for hydrothermal fluid. The petrophysical parameters, except for the permeability (i.e., κ, ϕ, ρ_m and $c_{\mathrm{p}}^{\mathrm{m}}$), are assumed to be the same as those in Iyer et al. (2010). We assume that log k exponentially decreases from −11.8 to −16.8 with a length scale of 300 m for an e-fold decrease, which is consistent with observations (Fisher, 1998, see the Supplement). Following Cherkaoui et al. (2003), oceanic rocks below the crust–mantle interface (6 km depth from the ocean–crust interface), and with temperatures above the critical temperature for rock cracking (600 ^∘C), are impermeable. As long as the model is consistent with observations, changing the assumptions about permeability does not affect the general conclusions given in the present study (see the Supplement).

The boundary conditions adopted in the present study follow those by Iyer et al. (2010) and/or Cherkaoui et al. (2003). Pressure at the ocean–crust interface is assumed to be constant at 25 MPa, corresponding to the assumption of a 2.5 km water depth. This assumption is modified when we explore the Precambrian hydrothermal circulation in Sect. 3.3. The right and bottom sides of the calculation domain are assumed to be insulating. The ocean above the top boundary is assumed to have a constant temperature of 2 ^∘C. The ridge axis on the left-hand boundary is assumed to supply a boundary heat flux J_b (J m⁻² yr⁻¹):

where w is the spreading rate (m yr⁻¹) and T_m is the temperature of the intrusion (1200 ^∘C). As a standard value, we assume $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{2}}$ m yr⁻¹ for w, and this assumption is again changed in Sect. 3.3 in which we consider the Precambrian hydrothermal circulation. Free flow is allowed at the top boundary, and the other boundaries are assumed to be impermeable (Iyer et al., 2010). The calculation domain is wide enough (see below) so that making the right-hand boundary permeable with respect to heat and water will change the results only negligibly (cf. Supplement; Cherkaoui et al., 2003). See Table 1 for the definitions and values of the parameters used in the present study.

Table 1Symbols and their definitions and values.

Download Print Version | Download XLSX

A finite difference approach is taken to solve Eqs. (1)–(3) for q, P and T using the second-order central differencing scheme for the second-order differential terms and the first-order upwind and forward-differencing schemes for the first-order spatial and temporal differential terms, respectively (e.g., Steefel and Lasaga, 1994). The calculation procedure follows that of Iyer et al. (2010). First, we obtain P by solving Eqs. (1) and (2). Then, q is obtained from the calculated P and Eq. (2). Finally, the calculated q is used in Eq. (3), which is solved to obtain T for the same time step. The calculated P and T are used to update the thermodynamic and transport properties of water for the next calculation step. The above procedure is repeated with a time step of 3×10⁴ years until steady states are reached, which is accomplished within 10³ time steps or 3×10⁷ model years (cf. Cherkaoui et al., 2003). Note that only steady states are considered in this study because of the long timescale over which oceanic oxygen isotopic composition can change ($>\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{7}}$ year; e.g., Holland, 1984; Walker and Lohmann, 1989; Wallmann, 2001). The model grid extends as deep as 12 km from the ocean–crust interface, allowing changes in the location of the crust–mantle interface (though not varied in this study); and it extends to a 30 km distance from the ridge axis, wide enough to simulate major changes in oxygen isotopic composition of oceanic crust which have been observed to be completed within less than 10 Myr from the rock formation at midocean ridges (Muehlenbachs, 1979, see the Supplement for more details). The calculation domain (12×30 km²) is divided into a 320×200 irregular grid, with the grid-cell size increasing horizontally from the ridge axis (1.1 to 330 m) and vertically from the ocean–crust interface (0.17 to 82 m; cf. Cathles, 1983; Cherkaoui et al., 2003).

2.2 Reactive transport model of oxygen isotopes

Oxygen isotopes of solid rock are assumed to be transported by the spreading of oceanic rocks, while those of porewater are transported by water flow and molecular diffusion (cf. Norton and Taylor, 1979; Lécuyer and Allemand, 1999). At the same time, isotope exchange reactions transfer ¹⁸O between the two phases (e.g., Norton and Taylor, 1979; Lécuyer and Allemand, 1999). Then, from the mass conservation of ¹⁸O in the two phases, the time rates of change in the ¹⁸O∕total O mole ratios of solid rock and porewater (F_r and F_p, respectively) can be represented by

where m_s and m_f are the oxygen concentrations of solid rock and porewater, respectively (mol kg⁻¹), ρ_b is the density of bulk rock given by ${\mathit{\rho }}_{\mathrm{b}}=\mathit{\varphi }{\mathit{\rho }}_{\mathrm{f}}+(\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}_{\mathrm{m}}$, k_ex is the rate constant for oxygen isotope exchange (${\mathrm{mol}}^{-\mathrm{1}}\mathrm{kg}{\mathrm{yr}}^{-\mathrm{1}}$), α is the oxygen isotope fractionation factor, and D is the effective diffusion coefficient for ¹⁸O (m² yr⁻¹) determined by molecular diffusion and hydrodynamic dispersion. The first and second terms on the right-hand sides of Eqs. (5) and (6) represent oxygen isotope transport and oxygen isotope exchange, respectively. The kinetic expression of oxygen isotope exchange in Eqs. (5) and (6) is formulated based on Cole et al. (1983, 1987) and the rate constant is given by an Arrhenius equation:

where E is the apparent activation energy (5×10⁴ J mol⁻¹) and R_g is the gas constant (8.314 $\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$). The value of 10^−8.5 ${\mathrm{mol}}^{-\mathrm{1}}\mathrm{kg}{\mathrm{yr}}^{-\mathrm{1}}$ at reference temperature 278 K (5 ^∘C) is obtained by reducing the range 10^−7.2–10^−6.6 ${\mathrm{mol}}^{-\mathrm{1}}\mathrm{kg}{\mathrm{yr}}^{-\mathrm{1}}$ at 5 ^∘C extrapolated from 10^−2.1–10^−1.5 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ at 300 ^∘C (Cole et al., 1987) with $E=\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$ J mol⁻¹ (cf. Cole et al., 1987) and 10³ m² kg⁻¹ specific surface area of marine basalt (cf. Nielsen and Fisk, 2010) by a factor of 10^1.3–10^1.9. The reducing factor (10^1.3–10^1.9) accounts for the fact that the reaction rate in the field is generally slower than in the laboratory by a factor of up to 10³ (e.g., Wallmann et al., 2008, see the Supplement for details). A similar factor (up to 10⁴) has been adopted by Cathles (1983) for the formulation of oxygen isotope exchange kinetics. Note that the general results and conclusions in the present study are not affected by changes in the reference k_ex value within a plausible range (see the Supplement). The oxygen isotope fractionation factor for andesite by Zhao and Zheng (2003) is adopted in the present study as follows:

where β=0.876 (Zhao and Zheng, 2003) because Eq. (8) with β=0.876 yields comparable α values to experimental results for basalt by Cole et al. (1987) at 300–500 ^∘C, while being applicable over the wider range of temperature considered for the present study (2–1200 ^∘C). The effective diffusion coefficient D considers both molecular diffusion and hydrodynamic dispersion; the former is obtained from the modified Stokes–Einstein relation for ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{16}}\mathrm{O}$ diffusion by Krynicki et al. (1978) and the isotope effect from Harris and Woolf (1980) and a homogeneous dispersivity of 10 m is assumed for the latter (Frind, 1982; Gelhar et al., 1992) as follows:

The first term on the right-hand side of Eq. (9) represents the molecular diffusion including the tortuosity factor τ=ϕ^1.4 (Aachib et al., 2004), while the second term represents the hydrodynamic dispersion. The petrophysical parameters are assumed to be the same as those in the hydrothermal circulation model (Sect. 2.1; Table 1). The thermodynamic and transport properties of water are obtained through the hydrothermal circulation model (Sect. 2.1).

The steady-state values of F_r and F_p are obtained by solving $\partial F_{\mathrm{r}}/\partial t=\partial F_{\mathrm{p}}/\partial t=\mathrm{0}$ in Eqs. (5) and (6). The intrusion on the left-hand boundary is assumed to have an ¹⁸O∕total O mole ratio of fresh crust (F_m) that corresponds to 5.7 ‰ relative to SMOW (e.g., Holmden and Muehlenbachs, 1993), and the ocean above the top boundary is assumed to have a constant ¹⁸O∕total O mole ratio of seawater (F_sw). The other boundaries are impermeable for ¹⁸O fluxes via water (see Sect. 2.1). A finite difference method is used for equation discretization (first-order upwind and second-order central differencing schemes for the first-order and second-order spatial differential terms, respectively) and Newton's method is adopted to solve the difference equations. The calculation is conducted on the grid described earlier for the hydrothermal circulation model (Sect. 2.1). The ¹⁸O mole ratios are all reported in δ notation relative to SMOW, using $\mathrm{2.0052}\times {\mathrm{10}}^{-\mathrm{3}}$ as the ¹⁸O∕¹⁶O mole ratio of SMOW (Fry, 2006).

3.1 Application to the present day and model validation

The calculated flow geometry and temperature distribution (Fig. 1; note that Fig. 1c and d are plots of Fig. 1a and b, respectively, on logarithmic scales) are similar to those in previous studies (e.g., Cherkaoui et al., 2003). The 2D model results can be converted to associated mass and heat fluxes, assuming 10⁸ m ridge length (cf. Wolery and Sleep, 1976). The total heat flux from the system is 0.74×10¹² W, which is comparable to the observed cumulative heat flux within 1 Myr (corresponding to 30 km with $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{2}}$ m yr⁻¹ spreading rate) from the ridge axis, $\mathrm{0.4}(\pm \mathrm{0.3})\times {\mathrm{10}}^{\mathrm{12}}$ W (Stein and Stein, 1994). The total mass of water exchange is 1.2×10¹³ kg yr⁻¹, falling within the constrained range (3×10¹² to 1.6×10¹⁴ kg yr⁻¹) and close to the recommended value ($\mathrm{3}(\pm \mathrm{1.5})\times {\mathrm{10}}^{\mathrm{13}}$ kg yr⁻¹) by Elderfield and Schultz (1996), and is also comparable to the prediction from other hydrothermal circulation models (e.g., Norton and Knight, 1977; Fehn et al., 1983; Cherkaoui et al., 2003). The calculated 2D distributions of solid rock and porewater δ¹⁸O (Figs. 2a, d, 3a and d; note that Fig. 3 is a plot of Fig. 2 on logarithmic scales) are consistent with observations, especially at 30 km (1 Myr) from the ridge axis (black solid curve in Fig. 4; see Sect. 3.2 for paler curves). Where flow rate and temperature are highest near the ridge axis and close to the ocean–crust interface (at $<\sim \mathrm{200}$ m depths), porewater δ¹⁸O is in the range between 0 ‰ and 2 ‰ (Fig. 3d), which is comparable to the observed δ¹⁸O ranges for high-temperature hydrothermal fluids, e.g., 0.49 ‰–2.3 ‰, by Jean-Baptiste et al. (1997), and 0.3 ‰–1.4 ‰, by James et al. (2014). Away from the ridge axis, porewater δ¹⁸O is equivalent to the seawater value at the ocean–crust interface, but drops down to −4 ‰ to −14 ‰ within ∼200 m depths, which is again not inconsistent with observed porewater δ¹⁸O values (down to −8 ‰) at <500 m depths where water exchange is limited (e.g., Lawrence and Gieskes, 1981). Solid rock δ¹⁸O at 1 Myr (i.e., 30 km from the ridge axis) is characterized by high (<16 ‰) and low (>3 ‰) values at shallow ($\le \sim \mathrm{2}$ km) and deeper ($\ge \sim \mathrm{2}$ km) depths, respectively, which is consistent with observations of modern oceanic crust (e.g., Alt et al., 1986) and ophiolites (e.g., Gregory and Taylor, 1981, crosses in Fig. 4; see Fig. S3 in the Supplement). The ¹⁸O fluxes to the ocean from high- and low-temperature alteration are 3.0×10⁹ and $-\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{9}}$ mol yr⁻¹, respectively, well balanced, resulting in a net ¹⁸O flux of 0.2×10⁹ mol yr⁻¹, which is consistent with the suggestion of zero net flux by Gregory and Taylor (1981). Individual ¹⁸O flux values are also comparable to those suggested in the previous studies. As examples, ¹⁸O fluxes through high-temperature alteration have been suggested to be 4.5×10⁹ (Holland, 1984), 3.2×10⁹ (Muehlenbachs, 1998) and 2.8×10⁹ mol yr⁻¹ (Wallmann, 2001), while those through low-temperature alteration are $<-\mathrm{2.3}\times {\mathrm{10}}^{\mathrm{9}}$ (Lawrence and Gieskes, 1981), $-\mathrm{1.1}\times {\mathrm{10}}^{\mathrm{9}}$ (Holland, 1984), $-\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{9}}$ (Muehlenbachs, 1998) and $-\mathrm{0.9}\times {\mathrm{10}}^{\mathrm{9}}$ mol yr⁻¹ (Wallmann, 2001). The consistency between the model calculation and observations described above supports the validity of the model.

Figure 1Two-dimensional distributions of hydrothermal fluid flow (a and c) and temperature (b and d). Shown in (a) and (c) are logarithms of fluid velocity (m yr⁻¹) together with mass-based stream lines that are depicted with white curves. The same data in (a) and (b) are plotted on logarithmic scales in (c) and (d), respectively. Gray zones in (a) and (c) represent where rocks are impermeable below 6 km depth and/or with temperatures above the rock-cracking threshold (600 ^∘C).

Download

3.2 Evaluation of oceanic δ¹⁸O buffering capacity

As both the source and sink of oceanic ¹⁸O, hydrothermal systems can buffer oceanic δ¹⁸O but the strength of the buffering depends on the sensitivity of isotope exchange between oceanic crust and porewater to seawater δ¹⁸O (cf. Wallmann, 2001). As an extreme example, if the oxygen isotope fractionation between solid rock and porewater is independent of seawater δ¹⁸O, there should not be any feedbacks from the hydrothermal systems on changes in seawater δ¹⁸O, i.e., no oceanic δ¹⁸O buffering. The previous studies that support the strong oceanic δ¹⁸O buffering at midocean ridges assume linear relationships between oxygen isotope fractionation and seawater δ¹⁸O (e.g., Gregory and Taylor, 1981; Holland, 1984; Muehlenbachs, 1998). This assumption regarding the sensitivity to seawater δ¹⁸O is difficult to confirm through observations of geological records because seawater δ¹⁸O is not known. We can examine the response of rocks to seawater δ¹⁸O with the present model by adopting different values for seawater δ¹⁸O. We can then measure the buffering capacity of the system by plotting net ¹⁸O flux against seawater δ¹⁸O and calculating the slope value, i.e., $\partial \left(\text{net }{}^{\mathrm{18}}\mathrm{O}\text{ flux}\right)/\partial \left(\text{seawater }\mathit{\delta }{}^{\mathrm{18}}\mathrm{O}\right)$, as in, e.g., Muehlenbachs and Clayton (1976). For example, a hydrothermal system with a large negative slope value should exhibit a strong buffering of oceanic δ¹⁸O because a slight change in seawater δ¹⁸O makes the system introduce a large net ¹⁸O flux to the ocean that restores the change in seawater δ¹⁸O (e.g., Muehlenbachs and Clayton, 1976; Muehlenbachs, 1998).

Figure 2Two-dimensional distributions of solid rock and porewater δ¹⁸O (a–c and d–f, respectively) at seawater δ¹⁸O values of 0 ‰, −6 ‰ and −12 ‰ (a and d, b and e, and c and f, respectively).

Download

Figure 3As for Fig. 2, except plotted on logarithmic scales.

Download

Figure 4Solid rock δ¹⁸O as a function of depth at 1 Myr (30 km) from the ridge axis at seawater δ¹⁸O values of 6 ‰, 4 ‰, …, −12 ‰. Also plotted are Oman ophiolite data from Gregory and Taylor (1981, crosses) converting their reported depth y^′ (km) to $\mathrm{7.4}-y^{\prime }$ (km) to facilitate comparison.

Download

Most features of 2D distributions of solid rock and porewater δ¹⁸O are not significantly affected by changing seawater δ¹⁸O from 6 ‰ to −12 ‰ (Figs. 2–4). Although the porewaters close to the ocean have δ¹⁸O compositions close to those of seawater (e.g., Fig. 3), δ¹⁸O signatures of solid rocks are not linearly proportional to seawater δ¹⁸O (e.g., Fig. 4). The relative insensitivity at shallow depths ($\le \sim \mathrm{2}$ km; Figs. 2–4) can be explained by the kinetics of oxygen isotope exchange. The distance from the isotope exchange equilibrium can be measured by $\mathrm{\Omega }=F_{\mathrm{r}}(\mathrm{1}-F_{\mathrm{p}})/\mathit{\left\{}\mathit{\alpha }\right(\mathrm{1}-F_{\mathrm{r}}\left)F_{\mathrm{p}}\mathit{\right\}}$ (Eqs. 5 and 6; Fig. 5b and d; note that Fig. 5c and d are plots of Fig. 5a and b, respectively, on logarithmic scales), and the shallow regions show nonequilibrium states (Ω≠1) because of their relatively low temperatures (Fig. 1b and d). Therefore, despite the approximate isotopic equivalence between porewater and seawater (Figs. 2 and 3), solid rocks do not directly reflect seawater δ¹⁸O in their oxygen isotopic compositions. In deeper sections of oceanic crust, on the other hand, solid rocks attain isotope exchange equilibrium with porewater because of the relatively high temperatures and the correspondingly fast rates of isotope exchange (e.g., Fig. 5b and d). These deep regions, however, receive less water exchange than the shallow regions, as can be seen from the distribution of local water∕rock oxygen mole ratio η ($\equiv m_{\mathrm{f}}\left|\mathit{q}\right|/\mathit{\left\{}\right(\mathrm{1}-\mathit{\varphi }\left){\mathit{\rho }}_{\mathrm{m}}{m}_{\mathrm{s}}w\mathit{\right\}}$; Fig. 5a and c). Accordingly, deep oceanic rocks and porewaters are oxygen isotopically buffered by spreading solid rocks rather than by seawater, and δ¹⁸O values of deep porewaters are almost independent of seawater δ¹⁸O (Figs. 2d–f and 3d–f). Consequently, despite the isotope exchange equilibrium, deep oceanic rocks are almost completely insensitive to seawater δ¹⁸O. The combination of the above two factors (isotope exchange kinetics and ¹⁸O supply from spreading solid rocks) explains the insensitivity of midocean ridge systems to seawater δ¹⁸O (Figs. 2–4). Note that the kinetics of oxygen isotope exchange and ¹⁸O supply from rocks have been invoked to explain the positive and negative δ¹⁸O anomalies, respectively, relative to the pristine crustal δ¹⁸O value by Cathles (1983), which is not inconsistent with the two factors causing the partial decoupling between oceanic crust and seawater δ¹⁸O in the present model.

Figure 5Two-dimensional distributions of local water∕rock oxygen mole ratio (a and c) and degree of oxygen isotope exchange at a seawater δ¹⁸O value of 0 ‰ (b and d). The same data in (a) and (b) are plotted on logarithmic scales in (c) and (d), respectively. See the caption of Fig. 1 for the explanation of gray zones in (a) and (c).

Download

Consistent with the insensitivity of oceanic rocks to seawater δ¹⁸O (Figs. 2–4), the net ¹⁸O flux from the midocean ridge systems is a very weak function of seawater δ¹⁸O (Fig. 6). Compared to the large buffering capacity suggested in the previous studies, which can be recognized from large negative slopes in Fig. 6, the intensity of oceanic δ¹⁸O buffering exhibited by the present model is weak; the magnitude of slope ($-\mathrm{0.4}\times {\mathrm{10}}^{\mathrm{9}}$ $\mathrm{mol}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{‰}}^{-\mathrm{1}}$) is smaller than previously assumed by a factor of up to >7. The reason why the previous studies suggested strong buffering is because these studies assume isotope exchange equilibrium between porewaters and oceanic rocks and equivalence between porewater and seawater at the same time. These assumptions are not satisfied simultaneously in midocean ridge systems in general, as can be anticipated from Figs. 2–5. The exception to the previous studies is Wallmann (2001), which shows a relatively small negative slope value (Fig. 6). Note, however, that Wallmann (2001) simulated weak buffering in a different way, i.e., by adopting low and temperature-independent rate constants for oxygen isotope exchange.

Figure 6Net ¹⁸O flux to the ocean from hydrothermal systems as a function of seawater δ¹⁸O. Slope values are denoted near the lines from previous studies and the standard simulation in this study.

Download

3.3 Application to the Precambrian

As the tectonics of the Earth likely evolved through time, the midocean ridge systems in the Precambrian could have been quite different from those in the present day. Indeed, differences in ocean volume and the spreading rate of oceanic crust have been discussed (e.g., Kasting et al., 2006; Korenaga et al., 2017). Here, we examine the oceanic δ¹⁸O buffering during the Precambrian by utilizing possible ranges of the spreading rate ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ to $\mathrm{30}\times {\mathrm{10}}^{-\mathrm{2}}$ m yr⁻¹) and ocean volume (1–5 km of water depth at the ocean–crust interface) as well as seawater δ¹⁸O (−12 ‰ to 6 ‰) suggested in the literature (e.g., Phipps Morgan, 1998; Kasting et al., 2006; Jaffrés et al., 2007; Korenaga et al., 2017; Galili et al., 2019; Johnson and Wing, 2020).

Figure 7Two-dimensional distributions of hydrothermal fluid flow (a, c, e, g, and i) and temperature (b, d, f, h, and j) from simulations with different spreading rates and ocean depths. Logarithms of fluid velocity and mass-based stream lines are depicted in (a, c, e, g, and i). The panels in the first three rows assume different spreading rates ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ m yr⁻¹ in a and b, $\mathrm{9}\times {\mathrm{10}}^{-\mathrm{2}}$ m yr⁻¹ in c and d, and $\mathrm{30}\times {\mathrm{10}}^{-\mathrm{2}}$ m yr⁻¹ in e and f), while those in the last two rows assume different ocean depths (1 km in g and h; 5 km in i and j). Values of other parameters are the same as those in the standard simulation (e.g., Fig. 1). See the caption of Fig. 1 for the explanation of gray zones in (a, c, e, g, and i).

Download

Increasing the spreading rate results in increases in the heat flux and associated water exchange between the crust and ocean, which is consistent with modern observations (e.g., Baker et al., 1996; Bach and Humphris, 1999): 0.50×10¹², 0.99×10¹² and 1.25×10¹² W of heat flux and 2.4×10¹², 4.9×10¹³ and 8.5×10¹³ kg yr⁻¹ of water exchange at spreading rates of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$, $\mathrm{9}\times {\mathrm{10}}^{-\mathrm{2}}$ and $\mathrm{30}\times {\mathrm{10}}^{-\mathrm{2}}$ m yr⁻¹, respectively (Fig. 7a–f). On the other hand, differences in temperature and water-flow distributions caused by those in ocean depth (1–5 km) are relatively minor (0.74×10¹² W and (1.2–1.3)×10¹³ kg yr⁻¹; Fig. 7g–j). The distributions of solid rock and porewater δ¹⁸O predicted under these different alteration conditions are correspondingly modified but show features common to those shown in the previous sections with the standard parameterization. Porewater δ¹⁸O is close to seawater δ¹⁸O only along the ocean–crust interface and becomes quite different from seawater δ¹⁸O at deep depths (see the Supplement). Solid rock δ¹⁸O is relatively insensitive to seawater δ¹⁸O in general, showing the relative addition and removal of ¹⁸O at shallow and deeper depths, respectively (Fig. 8). The buffering capacity exhibited by models that reflect changes in the spreading rate and ocean depth is different from that in the standard case (Fig. 9) but still weaker than any of those assumed in the previous studies (cf. Fig. 6). Note that in the previous studies, the buffering capacity is assumed to increase linearly with the spreading rate (Gregory and Taylor, 1981; Holland, 1984; Muehlenbachs, 1998; Wallmann, 2001), and thus the slope values in Fig. 9 need be compared with those in the literature (Fig. 6) multiplied by a factor that accounts for changes in the spreading rate. As an example, the slope value of $-\mathrm{1.9}\times {\mathrm{10}}^{\mathrm{9}}$ $\mathrm{mol}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{‰}}^{-\mathrm{1}}$ with a spreading rate $\mathrm{30}\times {\mathrm{10}}^{-\mathrm{2}}$ m yr⁻¹ in Fig. 9 should be compared with $-\mathrm{6}\times {\mathrm{10}}^{\mathrm{9}}$ (Wallmann, 2001), $-\left(\mathrm{14}\text{–}\mathrm{27}\right)\times {\mathrm{10}}^{\mathrm{9}}$ (Gregory and Taylor, 1981), $-\mathrm{24}\times {\mathrm{10}}^{\mathrm{9}}$ (Muehlenbachs, 1998), and $-\mathrm{31}\times {\mathrm{10}}^{\mathrm{9}}$ $\mathrm{mol}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{‰}}^{-\mathrm{1}}$ (Holland, 1984) (cf. Fig. 6). Thus, the buffering with the present model, assuming $\mathrm{30}\times {\mathrm{10}}^{-\mathrm{2}}$ m yr⁻¹ of spreading rate, is much weaker than assumed in the previous studies. More details of changes in the system behavior are described below.

Figure 8Solid rock δ¹⁸O as a function of depth at 30 km from the ridge axis at seawater δ¹⁸O values of 6, 4, …, −12 ‰ from simulations with different spreading rates and ocean depths. Broken lines denote 1 ‰ and 16 ‰ between which observed δ¹⁸O of ophiolites and/or oceanic crust ranges (Gregory and Taylor, 1981; Barrett and Friedrichsen, 1982; Cocker et al., 1982; Elthon et al., 1984; Alt et al., 1986, 1995; Agrinier et al., 1988; Schiffman and Smith, 1988; Vibetti et al., 1989; Lécuyer and Fourcade, 1991; Stakes, 1991; Bickle and Teagle, 1992; Holmden and Muehlenbachs, 1993; Muehlenbachs et al., 2003; Alt and Bach, 2006; Furnes et al., 2007; Gao et al., 2012). Standard parameter values are assumed, except for the parameterization denoted at each panel.

Download

When spreading is slower, seawater penetrates deeper into the crust because of a weaker boundary heat flux (Eq. 4), and thus hydrothermal flows are more homogenized but less intense when compared to those in the standard case (Fig. 7a). Although the changes of crustal δ¹⁸O from the pristine δ¹⁸O value (5.7 ‰) become larger over the whole crust (Fig. 8a), the sensitivity of solid rock to seawater δ¹⁸O remains small (Fig. 8a). The small sensitivity is caused because oxygen isotope exchange at shallow depths of oceanic crust is still kinetically limited despite the elongated time duration for oxygen isotope exchange (3 Myr; cf. Figs. S13 and S16 in the Supplement) and because the deeper section is buffered by spreading rocks rather than circulating seawater, given the simultaneous decreases in the amount of water exchange and the spreading rate. The lower O supply from the mantle combined with the small sensitivity makes the buffering intensity weaker than that with the standard spreading rate (Fig. 9). With the spreading rate high, on the other hand, the hydrothermal penetration depth becomes shallower as a result of enhanced boundary heat flux (Eq. 4) and hydrothermal flows are more intense and more localized close to the ocean and the ridge axis (Fig. 7c and e). Accordingly, ¹⁸O depletion from solid rocks near the ridge axis at relatively high temperatures and water∕rock ratios is not recovered even after continued reactions at low temperatures. Nonetheless, the sensitivity of solid rock to seawater δ¹⁸O does not significantly differ from that in the other spreading rate cases (Fig. 8b–c) because the kinetics of oxygen isotope exchange and ¹⁸O buffering by spreading oceanic rocks remain effective to decouple the shallow and deeper oceanic crust from seawater δ¹⁸O, respectively. Note that the time duration for oxygen isotope exchange is reduced with the higher spreading rate (≥0.1 Myr) but can still be long enough to complete major changes in the oxygen isotopic composition of oceanic crust (cf. Fig. S15 in the Supplement). Despite the weak rock sensitivity to seawater δ¹⁸O, the larger O supply from the mantle makes the buffering capacity larger than that with the standard spreading rate (Fig. 9). Overall, the system as a whole exhibits nonlinear and reduced sensitivity to changes in the spreading rate as recognized from the relationship between the spreading rate and buffering capacity in Fig. 9 (cf. Gregory and Taylor, 1981; Holland, 1984; Muehlenbachs, 1998; Wallmann, 2001).

Figure 9Net ¹⁸O flux to the ocean from hydrothermal systems as a function of seawater δ¹⁸O from simulations with different spreading rates and ocean depths. See the legend for the types of symbol and line for individual simulations. Slope values are denoted near the lines.

Download

When assuming a shallow ocean with the ocean–crust interface at a 1 km water depth, hydrothermal circulation becomes slightly stronger and occurs slightly closer to the ocean (Fig. 7g) (cf. Kasting et al., 2006). In contrast, in a deep ocean with the ocean–crust interface at a 5 km depth from the sea surface, the mixing of water becomes slightly weaker and occurs slightly deeper (Fig. 7i). In both cases, the above changes in the flow geometry and intensity are relatively minor. The resultant oxygen isotope behaviors are thus similar to the standard case (Figs. 4 and 8d–e) and the oceanic δ¹⁸O buffering becomes only slightly weaker in the shallow- and deep-ocean cases than in the standard case (Fig. 9). Notably, oxygen isotope exchange through low-temperature hydrothermal alteration is relatively enhanced in the shallow ocean case (Fig. 9), which is not inconsistent with the prediction by Kasting et al. (2006). However, the magnitude of the ocean depth effect is much smaller than anticipated by Kasting et al. (2006), probably because of the settings in the hydrothermal circulation model (Sect. 2.1) in which 2D distributions of permeability and boundary heat flux associated with the spreading rate (Eq. 4) may exert a more significant control over the hydrothermal penetration depth (see Fig. 7 and the Supplement).

4.1 Interpretation of ophiolites

Simulations conducted in the present study suggest that oceanic rocks are not significantly affected by changes in seawater δ¹⁸O under any plausible alteration conditions (Figs. 4 and 8). Reported δ¹⁸O values of ophiolites and/or oceanic crust range from ∼1 to 16 ‰ (broken lines in Fig. 8). By comparison, the simulated solid rock δ¹⁸O values fall within this range for seawater δ¹⁸O values $\le \sim -\mathrm{10}$, $\sim -\mathrm{8}$ to 0, $\ge \sim -\mathrm{2}$ and $\ge \sim -\mathrm{2}$ ‰ at a spreading rate of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$, $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{2}}$, $\mathrm{9}\times {\mathrm{10}}^{-\mathrm{2}}$ and $\mathrm{30}\times {\mathrm{10}}^{-\mathrm{2}}$ m yr⁻¹, respectively (Figs. 4 and 8). Accordingly, we conclude that the constant seawater δ¹⁸O at 0 ‰ is neither a necessary nor a sufficient condition for explaining the relatively invariant δ¹⁸O records of ophiolites. As ophiolite δ¹⁸O profiles can be affected more by alteration conditions (e.g., the spreading rate and permeability; Sect. 3.3 and the Supplement) with the control by seawater δ¹⁸O remaining relatively weak, feedbacks between the alteration parameters (e.g., spreading rate and permeability) could be more important in reproducing ophiolite records. Therefore, ophiolites may be interpreted to indicate the insensitivity of oceanic rocks to seawater δ¹⁸O, realized by the feedbacks between alteration parameters, and not necessarily a constant seawater δ¹⁸O. The weak buffering (Figs. 6 and 9) accompanying the partial decoupling between the oceanic crust and seawater δ¹⁸O (e.g., Fig. 8) shows that this interpretation is more plausible than the constant seawater δ¹⁸O.

4.2 Controls of δ¹⁸O in the Precambrian oceans

Although the buffering of oxygen isotopes through hydrothermal alteration of oceanic rocks is weaker than previously assumed under any alteration conditions (Sect. 3.3), it could have been relatively strong when the spreading rate was $\ge \sim \mathrm{30}\times {\mathrm{10}}^{-\mathrm{2}}$ m yr⁻¹ (Figs. 6 and 9). Such high spreading rate conditions could have been possible only during the earliest period of Earth's history (>3.5 Ga;  Phipps Morgan, 1998) or even impossible according to Korenaga (2006). Excluding this earliest period, the buffering must have been weak despite the uncertainties in alteration conditions (Figs. 6 and 9). The weak buffering must have allowed variations of seawater δ¹⁸O through other surficial processes that exchange ¹⁸O with seawater, most likely through continental weathering (e.g., Walker and Lohmann, 1989). As surface environments likely have significantly changed through the eons (see the Introduction), ¹⁸O fluxes from continental weathering could have correspondingly varied and might have contributed to the oceanic ¹⁸O budget more significantly than from hydrothermal alteration of oceanic crust during certain periods in Earth's history.

Previous studies examining oceanic δ¹⁸O evolution on geological timescales have utilized box-modeling approaches to account for oxygen isotope exchange from both continental weathering and hydrothermal alteration of oceanic crust (Godderis and Veizer, 2000; Goddéris et al., 2001; Wallmann, 2001; Kasting et al., 2006; Jaffrés et al., 2007). Among them, the studies that assume strong buffering at midocean ridges have had difficulty in modifying oceanic ¹⁸O budget through other surficial processes, including continental weathering (Godderis and Veizer, 2000; Goddéris et al., 2001). Other studies that assume weak buffering instead have shown the possibility of significant oceanic δ¹⁸O changes through modification of simplified continental weathering parameters with the Earth's age (Wallmann, 2001; Kasting et al., 2006; Jaffrés et al., 2007). Revisiting these previous box model studies with the weak buffering suggested here (e.g., Figs. 6 and 9), and constructing a process-based model for continental weathering to be coupled with the present model, will lead to a better understanding of controls on oceanic δ¹⁸O during the Precambrian.

4.3 Comparison with other models

The results in the present study differ significantly from those by Lécuyer and Allemand (1999), who utilized prescribed distributions of temperature and water∕rock ratio along with an isotope fractionation equation for a half-closed system (a system that is open to the water reservoir but not to the rock reservoir; Gregory et al., 1989; see the Introduction). The assumed range for water ∕ rock ratio (≤10) by these authors is comparable to the simulated range in the present study (η≤164; Fig. 5a and c). States of equilibrium and nonequilibrium in the bulk-rock-based formulation of isotopes in the present study correspond to transport- and reaction-limited alteration states, respectively, in the mineral-based isotope formulation in Lécuyer and Allemand (1999); the difference in the formulation of isotopic composition of bulk rock should not cause any significant differences in the oxygen isotope behavior of bulk rock (cf. Bindeman et al., 2019). The temperature dependence of the kinetic constant for oxygen isotope exchange is described with the Arrhenius equation, based on the data sets by Cole and Ohmoto (1986) in Lécuyer and Allemand (1999) and by Cole et al. (1983, 1987) in the present study, and thus there should be little difference in the temperature dependence of isotope exchange kinetics between the two studies. The size difference in the calculation domain cannot explain the difference between the present study and Lécuyer and Allemand (1999), as additional numerical experiments on a wider calculation domain (300 km) with artificially imposing off-axis water exchange (up to 5.3×10¹⁵ kg yr⁻¹) yield essentially the same results as those with the standard calculation domain width (30 km; see the Supplement). Although a thorough comparison cannot be made between the two studies because Lécuyer and Allemand (1999) have not provided detailed model results, it is likely that the isotopic effect of spreading crust has not been explicitly considered by Lécuyer and Allemand (1999). Should this be the case, the sensitivity of oceanic rocks to seawater δ¹⁸O and the buffering capacity at midocean ridges could be overestimated by the model of Lécuyer and Allemand (1999). More generally, the models that do not explicitly account for the transport of solid rocks (e.g., Taylor, 1977; Norton and Taylor, 1979; Criss et al., 1987; Gregory et al., 1989; DePaolo, 2006) could overestimate the coupling strength between oceanic crust and seawater δ¹⁸O and the oceanic ¹⁸O buffering capacity of hydrothermal systems. Application of these models to ophiolites may also underestimate the uncertainty in reconstructed ancient seawater δ¹⁸O that is caused by the partial decoupling between oceanic crust and seawater δ¹⁸O (see Sect. 4.1 and the Supplement).