Deep-water cycling and the Magmatic History of the Earth

8 Earth is a magmatically active planet. Magmatism connects Earth’s interior to its atmosphere, hy9 drosphere, and biosphere through cycling of volatiles, greenhouse gasses, and nutrients [18]. Earth’s 10 magmatic history is intertwined with its thermal and tectonic evolution. How magmatism has evolved 11 and been maintained in the face of planetary cooling remains an open question. We address this question 12 using data-constrained deep-water cycling and thermal history models. We track magmatic potential 13 using a homologous temperature: the ratio of upper mantle to melting temperatures. After an initial 14 decline, homologous temperature is buffered at a nearly constant value from roughly 2.5-2.0 Ga to the 15 present day. Melt buffering reflects two factors: 1) The dependence of melting temperature on water 16 content [21], and 2) The dependence of mantle viscosity on temperature and water content [15, 31, 27]. 17 The latter allows solid Earth evolution to self-regulate via feedbacks that keep mantle viscosity at a near 18 constant value. Self-regulation occurs even though the mantle remains far from thermal equilibrium, 19 consistent with heat flow data. The added feedback from water-dependent melting allows magmatism 20 to be co-buffered over geological time. This indicates that coupled thermal and water cycling feedbacks 21 have maintained melting on Earth and associated volcanic/magmatic activity. Magmatic self-regulation 22 affects not only the lifetime of geological activity on Earth but also, to the degree that planetary life 23 connects to volcanic activity, the maintenance of conditions favorable for life. 24


28
It has long been noted that the temperature of Earth's shallow mantle is remarkably close to the melting 29 temperature of rock [39]. That proximity (Figure 1a) is critical to Earth's current volcanic activity. It 30 could be a coincidence, in which case our planet's magmatic/volcanic activity will decline as it continues 31 to cool. More likely, it could reflect some form of feedback(s) that allow the Earth's cooling and magmatic 32 potential to be co-buffered. Magmatic/volcanic regulation over geologic time has not generally been from thermal equilibrium. This, in turn, has been used to argue that mantle convection is not self-52 regulated, which has implications not only for understanding our own planet's evolution, but also for 53 comparative planetology [25]. Although this argument is robust for thermal self-regulation, it does not 54 rule out self-regulation altogether. 55 The ability of a planet to self-regulate depends on a relationship between the vigor of mantle con-56 vection, as characterized by the mantle Rayleigh number (Ra), and convective heat flux (N u). That 57 relationship is given by and ρ is density, α is thermal expansivity, g is the acceleration due to gravity, ∆T is the superadiabatic 60 driving temperature, Z is the thickness of the convecting layer, κ is the thermal diffusivity and η is the (3)

Conservation of energy leads to
where C p is specific heat, V is mantle volume, H is mantle heat production, and Q s is surface heat flow.

95
Conservation of mantle water content leads to Following the assumption that viscosity remains statistically steady, relative to the time scale over which 97 significant changes occur in internal heat generation, leads to an estimate for the Urey ratio given by where η χ = ∂η ∂χ and η T = ∂η ∂T . If R exceeds D, then the Earth can be out of thermal equilibrium and low 99 values of U r are viable without requiring a weak, or negative, relationship between surface heat loss and 100 Ra. The analysis of Crowley et al.
[7] is significant in motivating our work, as it re-opens the possibility 101 of planetary self-regulation. It did not, however, show that the Earth followed a self-regulated path, 102 nor did it address magmatic evolution. In what follows we will do so using data-constrained thermal 103 evolution models that allow for water cycling.

105
Thermal history models have multiple free parameters and initial conditions. This can require millions 106 of model paths to map parameter space, with the vast majority of paths falling outside of observational 107 constraints [38,36]. The problem can be bridled via data assimilation. Here we develop and employ 108 a novel data assimilation method applied to coupled thermal history and deep-water cycling models.

109
The method directly builds in data constrained thermal history trajectories over geologic time (a full the mantle water capacitance. This, however, defines an upper limit within some uncertainty. We know 132 of no physical reasoning demanding that the mantle remain at this limit, and the majority of our results 133 fall within or below their uncertainties. Figure 2d shows the evolution of mantle viscosity from successful 134 models. In the absence of water, an expectation would be a monotonically increasing viscosity due to 135 mantle cooling. However, successful models show an increase in mantle viscosity for roughly the first half 136 of Earth's history followed by a milder decline. The rollover coincides with the change from net mantle 137 dewatering to rewatering. Physically, this corresponds to a switch from hot and dry subduction to cold 138 and wet subduction that cycles larger volumes of water into the mantle. passive margins that plate speeds have not been decreasing over geologic time and could, within data 145 uncertainty, even be increasing [1]. It is also consistent with the conjecture, based on observational constraints on water transport beneath Japan arcs, that deep-water cycling could stabilize and prolong 147 mantle convection and the associated geological activity of the Earth [16].

148
Another measure of self-regulation, beyond Ra, is the homologous temperature (T H ). We define T H 149 as the ratio of two depth-dependent profiles: the mantle geotherm and the wet solidus. Figure 1a shows 150 these profiles in black and blue, respectively. As the two profiles change with depth, we define T H at the 151 maximum distance between the two curves (see Methods). In Figure 1a where k is the thermal conductivity of the mantle, Z is the depth of the convecting layer, and ∆T is 250 the temperature difference driving convection. The latter is the difference between the surface temper- where Ra is and α is the mantle thermal expansivity, κ is the mantle thermal diffusivity and η is the mantle viscosity.

255
Viscosity depends on temperature and water content according to where η o is a scaling constant, c 1 , c 2 and c 3 are empirically determined constants [27], r is the water 257 fugacity exponent, A cre is a material constant, E is the activation energy for creep and R is the universal 258 gas constant. In Equation 14, χ m has units H/10 6 SI.

259
Combining Equations 7 to 13, we find that the change in mantle temperature evolves according to Rearranging to isolate mantle viscosity, Equation 15 becomes Equations 14 and 16 have η isolated. As such, we can use these equations to estimate χ m . We use h i , which will determine the rate of radiogenic heating for that model. Given the parameters in Table 1 267 and the constraints laid out here, χ m remains the only unknown in Equations 14 and 16. We initially  Table 1 lists the uncertainties we considered for each control point. We drew random We determined the constants α i by using the control points P o , P m1 and P f to form a system of three Solving for α i required a system of equations based on four control points. To account for this, we 288 introduced P m2 such that t m2 = t m1 − τ mt and T m2 = T m2 − τ mT . Here τ mt and τ mT represent offset 289 times and temperatures. These allow for flattening of the thermal trajectory after an initial temperature 290 decline, which can occur when water and thermal cycles effect mantle viscosity [37]. 291 We required that each thermal trajectory fit the data of Herzberg et al. [13] and Condie et al. [5] 292 within some measure of goodness. We used a reduced chi-squared statistic, the chi-square (χ 2 ) per degree 293 of freedom (ν). We adopt χ 2 as traditionally defined: This cumulatively measures the error (σ(t)) normalized difference between the data (D(t)) and modeled 295 thermal trajectory (M (t)). For measuring the goodness of fit we included all data points from both 296 data sets. This gave us a total of 38 data points (N d ). The definition for degrees of freedom is: 297 ν = N d − N α + N P given the number of parameters (N α ) and control points (N P ). We only kept 298 thermal trajectories that had χ 2 /ν <= 1. Using this method, we found a median accepted value of 0.98.

299
As this is nearly unity, so the thermal paths approximate the data error variance without over-fitting.

300
of an exponentially damped sine wave, which is of the form where A f is the amplitude of the sine wave, λ f is the decay constant, P f is the period and t is time, in 303 billion of years. We set A f to one percent of the initial mantle potential temperature, λ f to 0.22 Gyr −1 304 and P f to 1 Gyr. We then superimposed T f on top of a path defined as above. We still enforced the 305 condition that χ 2 /ν <= 1. We do no pretend to know what path fluctuations follow. They likely follow 306 something much more complicated than presented here. However, we believe that the qualitative form of 307 our findings will hold. An exact description of the fluctuations is beyond the scope of this paper. How to 308 account for them in forward modeling was covered by Seales et al. [38]. Regardless, choosing any other 309 path would find the same qualitative conclusions presented herein.
The convective profile is the mantle adiabat. We used a linearized version of Mckenzie and Bickle [30] 316 above to convert from T m to T p . We can also use this adiabat to construct the convective element of the 317 geotherm according to The conductive and convective profiles intersect at the base of the lithosphere (H L ): In our analysis we compared the geotherm to the solidus. We used the dry solidus of Hirschmann where the temperature is given in Kelvin and A 1 , A 2 , A 3 , K and γ are calibration constants.

323
Generally a homologous temperature is defined as the ratio of actual temperature to melting tem-324 perature. We follow this convention and define the homologous temperature (T H ) as the ratio of the 325 geotherm temperature to the solidus temperature. Both temperatures vary with depth, so we chose the 326 depth that maximized the thermal distance between them -the base of the lithosphere: This quadratic has two roots, one above the surface (unphysical) and one at depth. The one at depth 331 defines (H T ) and is given by We found the base of the melt zone (H B ) by equating the convective profile (Equation 21) with the solidus (Equation 24). Grouping like terms and gathering gives the quadratic This quadratic has two roots. The physical root is the shallower of the two. The deeper is an artifact of 333 the chosen solidus structure rather than anything physical. We define the base of the melt zone (H B ), 334 then, as Data constrain the solidus to a depth of 300 km. We set this as a hard maximum limit for H B .

337
The distance between the the solidus and geotherm determines the amount of melt produced. The