Quantifying the buffering of oceanic oxygen isotopes at ancient midocean ridges

To quantify the intensity of oceanic oxygen isotope buffering through hydrothermal alteration of the oceanic crust, a 2D hydrothermal circulation model was coupled with a 2D reactive transport model of oxygen isotopes. The coupled model calculates steady-state distributions of temperature, water flow and oxygen isotopes of solid rock and porewater given the physicochemical conditions of oceanic crust alteration and seawater δ¹⁸O. Using the present-day seawater δ¹⁸O under plausible modern alteration conditions, the model yields δ¹⁸O profiles for solid rock and porewater and fluxes of heat, water and ¹⁸O that are consistent with modern observations, confirming the model's validity. The model was then run with different assumed seawater δ¹⁸O values to evaluate oxygen isotopic buffering at the midocean ridges. The buffering intensity shown by the model is significantly weaker than previously assumed, and calculated δ¹⁸O profiles of oceanic crust are consistently relatively insensitive to seawater δ¹⁸O. These results are attributed to the fact that isotope exchange at shallow depths does not reach equilibrium due to the relatively low temperatures, and ¹⁸O supply via spreading solid rocks overwhelms that through water flow at deeper depths. Further model simulations under plausible alteration conditions during the Precambrian showed essentially the same results. Therefore, δ¹⁸O records of ophiolites that are invariant at different Earth ages can be explained by the relative insensitivity of oceanic rocks to seawater δ¹⁸O and do not require constant seawater δ¹⁸O through time.

1 Introduction

Hydrothermal alteration of oceanic crust at midocean ridges works as the dominant source and/or sink of several elements and/or isotopes in the ocean (e.g., Wolery and Sleep, 1976; Elderfield and Schultz, 1996). Notably, oxygen isotopes have been considered to be primarily controlled by isotope exchange at midocean ridges. The observation of oceanic crustal δ¹⁸O has revealed that oceanic rocks gain and lose ¹⁸O through low- and high-temperature hydrothermal alteration at shallow and deeper depths, respectively. The close balance between the addition and removal of the heavy isotope, together with the huge oxygen supply from the mantle, has led to a hypothesis that the water–rock interactions at midocean ridges have buffered oceanic δ¹⁸O at the present-day value (0 ‰ relative to standard mean ocean water – SMOW) throughout the Earth's history (e.g., Muehlenbachs and Clayton, 1976; Gregory and Taylor, 1981; Holland, 1984; Muehlenbachs, 1998). Age-invariant δ¹⁸O records of ophiolites (ancient oceanic crust) have been argued to support the hypothesis (e.g., Holmden and Muehlenbachs, 1993).

In contrast to ophiolites, sedimentary rocks (e.g., carbonates, cherts, phosphorites and shales) have shown secular δ¹⁸O increases with the Earth's age; sedimentary rocks in the Precambrian are depleted in ¹⁸O by as much as 10 ‰ compared to modern sediments (e.g., Shemesh et al., 1983; Karhu and Epstein, 1986; Bindeman et al., 2016), and the exception to this has rarely been reported (e.g., Blake et al., 2010). In theory, if constituting minerals of sedimentary rocks were formed in equilibrium with seawater, and if the later diagenetic or metamorphic modification of δ¹⁸O was negligible, δ¹⁸O records of sedimentary rocks can be utilized to infer surface temperatures of the past, using temperature-dependent isotope fractionation factors and assuming the δ¹⁸O value of seawater from which the constituting minerals formed. Indeed, several authors have suggested hot climates in the Precambrian from ¹⁸O-depleted sedimentary records assuming the present-day seawater δ¹⁸O (e.g., 70±15 ^∘C at 3.2–3.5 Ga and 50–60 ^∘C through the later Precambrian; Knauth and Lowe, 2003; Knauth, 2005). Such high temperatures are apparently at odds with glacial records observed through the Precambrian (e.g., Catling and Kasting, 2017, Chap. 11). If one excludes hot climates based on the glacial records, one must conclude that sedimentary δ¹⁸O records do not reflect seawater δ¹⁸O (e.g., instead reflecting diagenetic or metamorphic overprints) and/or that seawater δ¹⁸O has changed through time (e.g., Walker and Lohmann, 1989). The latter conclusion has recently been supported by Proterozoic–Phanerozoic marine iron (oxyhydr)oxides, whose δ¹⁸O shows a similar secular evolution to those in sedimentary rocks despite its weak temperature dependence (Galili et al., 2019). On the other hand, the secular change in seawater δ¹⁸O apparently conflicts with the result derived from ophiolites: an invariant seawater δ¹⁸O resulting from strong oceanic δ¹⁸O buffering at midocean ridges. Overall, one cannot rely on sedimentary δ¹⁸O to reconstruct ancient surface temperatures and discuss associated climate forcings unless δ¹⁸O records of sedimentary rocks and ophiolites can be explained at the same time (cf. Godderis and Veizer, 2000; Vérard and Veizer, 2019).

Several hypotheses have been put forward to reconcile apparently conflicting invariant and variant δ¹⁸O records of respective ophiolites and sedimentary rocks. Perry et al. (1978) pointed out that oceanic rocks in the Precambrian were more mafic than in the present day and this could have resulted in a more intense low-temperature alteration of oceanic crust, removing more ¹⁸O from seawater than today. Walker and Lohmann (1989) argued that the midocean ridges could have been above sea level in shallower Precambrian oceans. Resultant subaerial water–rock interactions at low temperatures could then have removed a large amount of ¹⁸O from the Precambrian oceans. Even if midocean ridges were below the sea level, lower pressures on the seafloor in the shallow oceans could have brought less water into the oceanic crust, and water–rock interactions at midocean ridges could have been dominated by those at low temperatures at shallower depths, again resulting in lower oceanic δ¹⁸O in the Precambrian (Kasting et al., 2006). The absence of biogenic sediment cover during the Precambrian could also have contributed to more significant low-temperature seafloor alteration and a lower seawater δ¹⁸O (e.g., Jaffrés et al., 2007). The above arguments, however, give possible explanations only for low oceanic δ¹⁸O during the Precambrian; they do not explain age-invariant ophiolite δ¹⁸O records. Wallmann and colleagues have indicated that the invariant δ¹⁸O records of ophiolites can be explained by decoupling ancient oceanic crust from the contemporaneous seawater δ¹⁸O (e.g., Wallmann, 2001; Jaffrés et al., 2007). However, the simple box models considered by these authors cannot explain whether and/or how this decoupling could have been made possible (Jaffrés et al., 2007).

To resolve the issue, a process-based approach that simulates oxygen isotope behavior during oceanic crust alteration is indispensable. Lécuyer and Allemand (1999) developed an oxygen isotope exchange model, which utilizes prescribed distributions of temperature and water∕rock ratio within oceanic crust and an equation for half-closed systems (e.g., Gregory et al., 1989), and is thus not entirely process based. Using this model, Lécuyer and Allemand (1999) concluded that oceanic δ¹⁸O cannot change from the present-day value by more than 2 ‰ due to the strong buffering exhibited by their modeled isotope exchange at midocean ridges. This conclusion contrasts with that of Wallmann (2001), who adopted a box model to argue that the buffering must have been weak and that the oxygen isotopic composition of the ocean has evolved through the Phanerozoic. These conflicting conclusions from different models emphasize that a mechanistic understanding of oxygen isotope behavior during oceanic crust alteration is lacking.

The present study develops a process-based model to simulate oxygen isotope behavior during hydrothermal alteration of oceanic crust and presents its application to the Precambrian. A 2D hydrothermal circulation model presented in the literature is combined with a 2D reactive transport model of oxygen isotopes. The coupled model can thus yield 2D distributions of solid rock and porewater δ¹⁸O reflecting alteration conditions, including seawater δ¹⁸O and hydrothermal fluid circulation. After confirming the model's validity by comparing the model results that assume the present-day seawater δ¹⁸O with modern observations, we examine the intensity of oceanic δ¹⁸O buffering by hydrothermal alteration of oceanic crust at midocean ridges by changing seawater δ¹⁸O. The buffering quantification is then conducted under different physicochemical conditions, which could have been the case during the Precambrian, to give insights into how oceanic δ¹⁸O could have been affected by water–rock interactions at midocean ridges during the Earth's early eons.

2 Methods

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

$$\begin{array}{}\text{(1)}& \mathrm{\nabla }\cdot \mathit{q}=\mathrm{0},\end{array}$$

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):

$$\begin{array}{}\text{(2)}& \mathit{q}=-{\displaystyle \frac{k}{\mathit{\nu }}}\left(\mathrm{\nabla }P-{\mathit{\rho }}_{\mathrm{f}}\mathit{g}\right),\end{array}$$

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:

$$\begin{array}{}\text{(3)}& \left\{\mathit{\varphi }{\mathit{\rho }}_{\mathrm{f}}{c}_{\mathrm{p}}^{\mathrm{f}}+\left(\mathrm{1}-\mathit{\varphi }\right){\mathit{\rho }}_{\mathrm{m}}{c}_{\mathrm{p}}^{\mathrm{m}}\right\}{\displaystyle \frac{\partial T}{\partial t}}=\mathrm{\nabla }\cdot \left(-\mathit{q}{c}_{\mathrm{p}}^{\mathrm{f}}T+\mathit{\kappa }\mathrm{\nabla }T\right),\end{array}$$

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⁻¹):

$$\begin{array}{}\text{(4)}& J_{\mathrm{b}}={\mathit{\rho }}_{\mathrm{m}}w\left(T_{\mathrm{m}}-T\right)c_{\mathrm{p}}^{\mathrm{m}},\end{array}$$

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 1 Symbols and their definitions and values.

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left(\mathrm{1}-\mathit{\varphi }\right){\mathit{\rho }}_{\mathrm{m}}{m}_{\mathrm{s}}{\displaystyle \frac{\partial F_{\mathrm{r}}}{\partial t}}=-\left(\mathrm{1}-\mathit{\varphi }\right){\mathit{\rho }}_{\mathrm{m}}{m}_{\mathrm{s}}w{\displaystyle \frac{\partial F_{\mathrm{r}}}{\partial x}}\\ \text{(5)}& {\displaystyle }& {\displaystyle }-{\mathit{\rho }}_{\mathrm{b}}{m}_{\mathrm{s}}{m}_{\mathrm{f}}{k}_{\text{ex}}\left\{F_{\mathrm{r}}\left(\mathrm{1}-F_{\mathrm{p}}\right)-\mathit{\alpha }\left(\mathrm{1}-F_{\mathrm{r}}\right)F_{\mathrm{p}}\right\},\text{ and }{\displaystyle }& {\displaystyle }\mathit{\varphi }{\mathit{\rho }}_{\mathrm{f}}{m}_{\mathrm{f}}{\displaystyle \frac{\partial F_{\mathrm{p}}}{\partial t}}=\mathrm{\nabla }\cdot \left(-m_{\mathrm{f}}\mathit{q}{F}_{\mathrm{p}}+\mathit{\varphi }{\mathit{\rho }}_{\mathrm{f}}{m}_{\mathrm{f}}D\mathrm{\nabla }{F}_{\mathrm{p}}\right)\\ \text{(6)}& {\displaystyle }& {\displaystyle }+{\mathit{\rho }}_{\mathrm{b}}{m}_{\mathrm{s}}{m}_{\mathrm{f}}{k}_{\text{ex}}\left\{F_{\mathrm{r}}\left(\mathrm{1}-F_{\mathrm{p}}\right)-\mathit{\alpha }\left(\mathrm{1}-F_{\mathrm{r}}\right)F_{\mathrm{p}}\right\},\end{array}$$

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:

$$\begin{array}{}\text{(7)}& k_{\text{ex}}={\mathrm{10}}^{-\mathrm{8.5}}\exp \left\{-{\displaystyle \frac{E}{R_{\mathrm{g}}}}\left({\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{278}}}\right)\right\},\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }{\mathrm{10}}^{\mathrm{3}}\ln \mathit{\alpha }=& {\displaystyle }\mathit{\beta }\left({\displaystyle \frac{\mathrm{6.673}\times {\mathrm{10}}^{\mathrm{6}}}{T^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{10.398}\times {\mathrm{10}}^{\mathrm{3}}}{T}}-\mathrm{4.78}\right)\\ {\displaystyle }& {\displaystyle }\cdot \exp \left({\displaystyle \frac{\mathrm{1}-\mathit{\beta }}{R_{\mathrm{g}}T}}\right)-{\displaystyle \frac{\mathrm{2.194}\times {\mathrm{10}}^{\mathrm{6}}}{T^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{15.163}\times {\mathrm{10}}^{\mathrm{3}}}{T}}\\ \text{(8)}& {\displaystyle }& {\displaystyle }+\mathrm{4.72}+\mathrm{1.767}\left(\mathrm{2}\mathit{\beta }-\mathrm{1}\right),\end{array}$$

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:

$$\begin{array}{}\text{(9)}& D=\mathit{\tau }\mathrm{6.9}{\displaystyle \frac{T}{{\mathit{\rho }}_{\mathrm{f}}\mathit{\nu }}}\sqrt{{\displaystyle \frac{\mathrm{18}}{\mathrm{20}}}}+\mathrm{10}{\displaystyle \frac{\left|\mathit{q}\right|}{\mathit{\varphi }{\mathit{\rho }}_{\mathrm{f}}}}.\end{array}$$

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 Results

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 1 Two-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).

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 2 Two-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).

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

Figure 4 Solid 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.

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 5 Two-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).

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 6 Net ¹⁸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.

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 7 Two-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).

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 8 Solid 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.

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 9 Net ¹⁸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.

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 Discussion

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).

5 Conclusions

The 2D reactive transport model of oxygen isotopes combined with 2D hydrothermal circulation simulations enables us to predict distributions of temperature, water flow and oxygen isotopes of solid rocks and porewaters within oceanic crust based on mass, momentum and energy conservations. The model assuming the present-day seawater δ¹⁸O reproduced those distributions consistent with modern observations, supporting the model's validity. The intensity of seawater δ¹⁸O buffering at midocean ridges was evaluated by calculating the net ¹⁸O flux as a function of seawater δ¹⁸O. The buffering intensity predicted by the model is significantly weaker than previously assumed under any plausible alteration conditions during the Precambrian. The weak buffering is realized because isotope exchange equilibrium is not reached in shallow sections of oceanic crust with low temperatures, and because ¹⁸O supply via spreading solid rocks exceeds that through hydrothermal circulation in deeper high-temperature sections. Consistent with the weak buffering, oceanic rocks are insensitive to seawater with respect to oxygen isotopes. Thus, ophiolites can alternatively be explained by the insensitivity of oceanic rocks to seawater δ¹⁸O that could have evolved through the Precambrian because of the weak buffering at midocean ridges.