The hydrothermal power of oceanic lithosphere

Introduction Conclusions References


Introduction
The cooling of oceanic lithosphere over time as well as distance from ridges is a key constraint on plate tectonics.A predicted consequence of this cooling is that surface heat flux over young seafloor is relatively high and gradually diminishes with age.Although elevated heat flow on average has long been recognized over ridges (Bullard, 1952;Von Herzen, 1959;Von Herzen and Uyeda, 1963;Sclater, 2004), it was also found that measurements are highly scattered.Much subsequent work attempted to explain geophysical observations with models of lithospheric cooling and subsidence, and models were gradually refined by improved constraints on the geophysical properties of the upper mantle (e.g., Langseth et al., 1966;McKenzie, 1967;McKenzie and Parker, 1967;Le Pichon, 1968;Sleep, 1969;Sclater and Francheteau, 1970;Sclater et al., 1971).Eventually, it was recognized that measured heat flow close to ridge axes was far too low to be explained by a model that also explained seafloor subsidence accurately (Sclater and Francheteau, 1970).Following the work of Lister (1972Lister ( , 1974) ) and Bodvarsson and Lowell (1972), it became understood that the deficit between measured and predicted heat flow (as well as its scatter) originated in advective hydrothermal ventilation of heat between crustal basement and the oceans.This difference between model and measurement is therefore a proxy to estimating the lithospheric thermal power removed by hydrothermal circulation in oceanic lithosphere.Such constraints on the heat deficit are critically important for the understanding of numerous chemical and physical processes in the Earth.In addition to providing a direct constraint on the thermal budget of cooling oceanic crust and lithosphere (Mottl, 2003;Hasterok, 2013a), this is also the power available to drive chemical exchange between the crust and oceans (e.g., Seyfried et al., 1984;Spivack and Edmont, 1987;Nicolas et al., 2003;Staudigel, 2014) and to pass nutrients to sub-seafloor microbial communities (e.g., Jannasch, 1983Jannasch, , 1995;;Hessler et al., 1988;Tunnicliffe, 1991).In addition, the chemical alteration of oceanic crust by hy-drothermal fluids is thought to regulate the chemical budgets of subducting slabs (Schmidt and Poli, 2013;Ryan and Chauvel, 2013).Moreover, volatile cycling between the mantle and hydrosphere may impart a major control on secular changes in the efficiency of mantle convection over Earth history (Crowley et al., 2011).
Although the mentioned heat flow deficit exists in the "unfiltered" global database, Hasterok et al. (2011) and Hasterok (2013a, b) have demonstrated that the deficit is markedly reduced when only marine environments known to have thick sediment cover are considered.Such environments are expected to restrict advective heat transport because hydrothermal ventilation is generally confined to sites of outcropping basement, and thick sediments reduce the prevalence of outcrops.Moreover, these authors found that a compilation of heat flow estimates from young (< 30 My old) seafloor with thick sediment cover and extensive geophysical characterization agreed well with conventional lithospheric cooling models (e.g., model GDH1 of Stein and Stein, 1992).
Recently, models of the oceanic lithosphere have been developed that include constraints on the properties of the mantle and crust from mineral physics experiments (Hasterok, 2010;Grose and Afonso, 2013).These models predict that lithospheric heat flow over young seafloor may be significantly lower than that estimated with conventional models.Hasterok (2010) investigated lithospheric cooling models which account for differing properties of crust and mantle.Since the crust is relatively insulating, the effective heat transport properties of the lithosphere are "crust dominated" near ridge axes (a thermal blanketing effect yielding lower heat flow) and gradually approach mantle values over time.Furthermore, based on work characterizing crustal accretion and cooling, as well as hydrothermal heat transport (e.g., Chen and Phipps Morgan, 1996;Cherkaoui et al., 2003), Cochran and Buck (2001) and Spinelli and Harris (2011) showed that hydrothermal transport near ridge axes also results in low heat flow on ridge flanks away from the axis because, after the cessation of hydrothermal advection, the crustal geotherm rebounds from the effects of advective cooling.Finally, Grose and Afonso (2013) showed that axial hydrothermal circulation and crustal insulation together result in a compounded reduction of predicted heat flow in proximity to ridge axes.Since the most robust estimates of net ventilated hydrothermal power come from the difference of predicted and measured seafloor heat flow (Mottl, 2003), lower predicted conductive seafloor heat flow suggests that ventilated hydrothermal power in oceanic lithosphere must also decrease.The purpose of this paper is to examine these modified predictions of total ventilated hydrothermal power, and its spatial distribution, using recent models of lithospheric cooling.We will show that total ventilated hydrothermal power may be significantly less than conventional estimates, and that the fraction of hydrothermal power extracted on axis may be up to 50 % of the total.Thus, the extent of off-axis hy-drothermal circulation may be less than previously thought.Finally, as the resulting predictions of ventilated hydrothermal power (and its implications) depend on the accuracy of cooling models, we discuss constraints from global seafloor topography and several regions of young and geophysically well-characterized seafloor which are thought to constrain the lithospheric heat budget.

The power deficit
Heat lost from the cooling of oceanic crust and lithosphere must ultimately be transported through upper crustal layers and through (or around) sediment cover.If hydrothermal circulation occurs in the crust, and especially if hydrothermal fluids are exchanged between crust and oceans, measured heat flow can be substantially affected.Figure 1 shows an illustration summarizing perturbations of seafloor heat flow by hydrothermal circulation processes.We depict three possible thermal states in which oceanic crust may be regionally characterized.Firstly, if circulation in the crust is sluggish or is not unable to occur due to crustal impermeability, then lithospheric heat loss is largely conductive and measured heat flow at the surface should match the predictions of a conductive reference model.Secondly, if circulation does occur in the crust but fluids are not vented to the oceans, all heat removed from the surface occurs conductively, but closed circuit circulation will act as a high-conductivity layer, elevating surface heat flow above that expected by a model which does not include such a high-conductivity layer.This effect can be increased by increasing the thickness of the convective layer, or by hydrothermal exchange with deeper crust by advection along faults.We may call such a state superconductive.Finally, and most importantly, if vigorous hydrothermal circulation occurs in the crust and fluids are exchanged between the crust and oceans, the conductive heat flux above such an aquifer will be less than the average heat flux below it.This occurs because as heat is discharged to the oceans it is replaced by a recharge of cool water which is advected through a permeable aquifer.This reduces the temperature of aquifer as well as crust and sediment above it, lowering surface heat flow.This is the sub-conductive state in Fig. 1.Because measured heat flow over young seafloor is known to be systematically lower than that predicted by lithospheric models, it is believed that seafloor is dominantly in such a "sub-conductive" state.Also, as this deficit between measured and modeled heat flow is a consequence of fluid discharge and recharge, we can use it as a proxy to estimating the power lost to the oceans by ventilated hydrothermal circulation (e.g., Wolery and Sleep, 1976;Sleep and Wolery, 1978;Stein and Stein, 1994;Davis and Elderfield, 2004).
Following much prior work on this issue (see above references), we estimate the total power deficit by finding this deficit for each age and integrating over all ages as where Q H is the net power deficit, q m is the modeled heat flow, q o is the observed heat flow, A is the seafloor area, t is age, and t m is a maximum integration age.The three key variables which must be known are thus the measured and modeled heat flow, and the seafloor area-age distribution.We take the area-age distribution from the empirical model of Müller et al. (2008).The observed heat flow is based on the raw global heat flow database updated by Hasterok (2010).
As shown in Fig. 1, in sub-conductive regions there are locations of discharge and recharge.Ideally, we would wish to integrate heat flow over all possible seafloor area, including such anomalous sites.However, because characterizing heat flow around sites of discharge has been a major interest of marine heat flow surveys, such biases must be removed in order for the data to represent an average conductive heat flux from oceanic lithosphere.We remove measurements which are > 3000 W m −2 since such measurements indicate very close proximity to focused hydrothermal vents.We also remove high-resolution heat flow surveys with dense sampling over young seafloor with extensive thick sediment cover (e.g., Hobart et al., 1985;Davis et al., 1997Davis et al., , 1999) ) and an anomalous sampling of a mud volcano in the Barents Sea (Kaul et al., 2006).Also, all points for which seafloor age cannot be determined are removed.Hasterok (2013b) employed additional filters, including a thermal correction for thickening sediment (i.e., recently deposited sediment is initially at ocean-bottom temperature and must be gradually heated by conduction of lithospheric heat) and the removal of seafloor area with large igneous provinces.We do not include these for reasons of simplicity, our uncertainty in the accuracy of the thermal correction (which increases measured heat flow) and the possibility of other unaccounted corrections which would decrease measured heat flow or systematically increase predicted heat flow (e.g., thermal properties of sediment cover, internal heating sources, or non-plate-like reheating phenomena).Moreover, the thermal correction for sedimentation more strongly impacts observations over old seafloor, where hydrothermal ventilation is probably unimportant.Nevertheless, we will compare our results with those of Hasterok (2013b) to highlight the differences which ought to originate in our respective methodological choices.
In addition to constraints on the average seafloor heat flow, the deficit is calculated in respect to a given reference model of lithospheric heat flow.However, as illustrated in Fig. 1, if hydrothermal circulation occurs without ventilation, heat flow will be higher than that predicted by a "conductive" reference model.This super-conductive state occurs because an aquifer behaves as a high-conductivity layer.Thus, sub-conductive regions must also experience this phenomenon.Therefore, the proper reference model against which the power deficit should be taken must actually be a super-conductive reference.In other words, the reference model should ideally include a high-conductivity layer which will increase surface heat flux.Nevertheless, we do not include the effects of passive circulation for simplicity.While we will investigate a model with hydrothermal circulation near ridge axes, passive circulation does not occur in any of our reference models.Because hydrothermal circulation is thought to be largely confined to an aquifer which is only a few hundred meters thick (Spinelli and Harris, 2011), this fraction is likely small.On the other hand, if deep circulation is important, then our methods will underestimate ventilated heat loss.

Statistical analysis
Similar to Hasterok (2013b), we perform a Monte Carlo statistical analysis to estimate uncertainties in heat flow and the power deficit.The heat flow database (with our filters as described above) is divided into 1 My bins, and we construct cumulative distribution functions (CDFs) of heat flow for each bin.The fraction of observations (heat flow measurements) within the bin t having heat flow less than Q is given by the (discontinuous) CDF: where t is the bin interval (1 My), o i (q ≤ Q) is a measurement with heat flux q less than Q, and N o is the total number of observations in the bin t.Figures illustrating the cumulative distribution functions of heat flow for each 1 My bin are shown in Fig. 2. By randomly sampling the CDF (which is characterized by values between 0 and 1), we obtain samples of heat flow consistent with the frequency (or probability) of its occurrence in the data.Therefore, to estimate the power deficit via Eqs.( 1) and ( 2), (a) the CDF is randomly sampled once for each and every bin, up to a maximum age t m of 100 My; (b) we take the difference of the resulting statistical sample of the heat flow observations and the mean model heat flow in the bin; and we (c) multiply the difference by the total seafloor area which exists for the age bin.The resulting equation for the cumulative power deficit is the sum of the heat flow deficits for each bin multiplied by their respective seafloor areas, or where qo is the statistically sampled heat flow in the bin, and A t is the total seafloor area in the bin.We note that in our analysis we use the mean over the time interval to estimate q m for the bin.Hasterok (2013b) argued that the median is more appropriate, leading to different results (lower heat flow deficits).Finally, the statistics of the power deficit can be represented by probability density functions (PDFs) by repeating the integration of Eq. (3) 10 6 times.The results are shown in Fig. 4. Also, to illustrate how the cumulative heat flow deficit changes as the age of integration increases from 1 to 60 My, we record PDFs for each age bin.Again, note that these PDFs represent probabilities for the cumulative (Eq. 3) heat flow deficit, not the deficit in individual bins.

Reference models for seafloor heat flow
For plate and half-space models, the relationship between thermal properties of the lithosphere and seafloor heat flux may be given as where [ρC p ] is the volume heat capacity, T is the difference of mantle and surface temperature, D is the thermal diffusivity, and q a is the adiabatic part of the heat flux.The subscript "ei" denotes that the properties are "transient effective" properties: they are "effective" properties because they are weighted in terms of the depth-integrated change in temperature of the mantle, and they are "transient" because they are in terms of the instantaneous change in the depth-integrated temperature of the mantle (see Appendices).Moving the time t to the left-hand side, we can define the coefficient which may be referred to as a rate of heat flow diminution over time.This coefficient has been useful to characterize simple lithospheric cooling models (e.g., g(t < 20 My) = 510 mW m −2 My 1/2 for model GDH1; Stein and Stein, 1992) as it is basically a term isolating the thermal coefficients.In addition, this coefficient is particularly useful for analysis of lithospheric models with complicated heat transport properties.In such cases, the effective thermal properties have a characteristic time-dependence, thus yielding a timedependent g(t) coefficient.Subtle variations in these properties are much easier to discern in g(t) than in q sf .Cooling half-space models with no properties that depend on depth will have a constant g value, whereas g(t) for plate models will increase when the plate boundary is sensed, and g(t) will be a continuous function of time for any thermal model with depth-dependent properties.
We test three models of predicted heat flow, from which we calculate the heat flow deficit.The first plate model we use is from Hasterok (2013a), referred to here as H13.H13 has constant thermal properties and is thus taken as an optimal model prediction when the effects of thermal insulation and hydrothermal circulation are not considered.Model H13 has been constrained based on heat flow data only and predicts heat flow nearly equal to GDH1 (Stein and Stein, 1992) for ages < 50 My and slightly higher heat flow at older ages.Since H13 has constant thermal properties, it has a diminution coefficient which is constant over young (< 50 My old) seafloor.The second and third models we consider are from our previous work (Grose and Afonso, 2013), here referred to as GH and GHC.An illustration of initial and boundary conditions for these models is shown in Fig. 3.Both models GH and GHC have complicated thermal properties, including temperature-and pressure-dependent heat transport properties and thermal expansivity, thermal radiation heat transport, an initial geotherm calculated with an adiabat and latent heat of melting, and 2-D conductive heat transport.The physical details of the models are discussed in Appendix A. The properties of particular importance to models GH and GHC are the axial hydrothermal circulation and insulating properties of oceanic crust.The model name GHC indicates that this model includes both axial hydrothermal circulation and crustal insulation, while model GH only has axial hydrothermal circulation (in addition to the other properties noted above).
Axial hydrothermal circulation is modeled using a Nusselt (Nu) number approximation (e.g., Cochran and Buck, 2001;Spinelli and Harris, 2011), so that the effective thermal conductivity is Nu times the lattice value (Appendix A).High Nu occurs in lithosphere where the temperature is < 800 • C, the crust age is < 0.2 My, and we assume Nu = 10 (Fig. 3).These values are consistent with the "cracking front" limit for fast-spreading ridges estimated by Manning et al. (2000), and they produce a thermal structure near ridge axes which is in good agreement with predictions from seismic models for the East Pacific Rise (EPR; Dunn et al., 2000).Values which produce greater cooling rates (e.g., Spinelli and Harris, 2011;Maclennan et al., 2005;Cochran and Buck, 2001;Cherkaoui et al., 2003) have been shown to predict an excessively cool ridge axis (Grose and Afonso, 2013).
Crustal insulation is a consequence of the low thermal conductivity of crustal rocks.Models GH and GHC uses recent measurements of thermal properties, but we note that the resulting thermal conductivity is nearly equal to that estimated for basalts and gabbros by Zoth et al. (1988).
Models H13, GH, and GHC are investigated here for their predictions of ventilated hydrothermal power by means of the heat flow deficit.In addition, when discussing seafloor subsidence we will also consider the predictions of a model G and GC, which have all the properties of GHC except that model GC does not have axial hydrothermal circulation and model G has neither hydrothermal circulation nor crustal insulation.Specifically, model G is model RN1, model GH is model RN10, model GC is model RN1C, and model GHC is model RN10C from Grose and Afonso (2013).The key features of models GH and GHC (as well as model GC) are their strongly variable g(t) coefficients, shown in Fig. 5b.All models have been constrained with old-age seafloor heat flow.However, due to the effects of axial hydrothermal circulation and crustal insulation, g(t) (and thus heat flow) deviates significantly over young ages.Thus, comparing predictions of models H13, GH, and GHC should illustrate the consequences of axial hydrothermal circulation and crustal insulation on total ventilated hydrothermal power, and its spatial distribution.

Spatial power distribution
We refer to the deficit integrated to 1 My as the near-axial deficit.This deficit may be further divided into active and ridge flank fluxes.The active flux is the hydrothermal power vented over the ridge axis which is driven primarily by the emplacement and cooling of melts in the crust (Lister, 1982).After 0.2 My, active hydrothermal circulation ceases and the geotherm is in a state of conductive rebound (Cochran and Buck, 2001;Spinelli and Harris, 2011).While model H13 is a simple conductive plate model without hydrothermal circulation, we use 0.2 My as a cutoff value to estimate the active flux for this model also.
To obtain the active deficit we calculate the fraction of heat loss in the 0-0.2My bin as where the integral bounds are ages in millions of years (My).
From here the active flux 1 My), and the ridge flank flux is My). N A values for models H13, GH, and GHC are 0.435, 0.852, and 0.829, respectively.N A values are higher for models GH and GHC due to high q m from hydrothermal circulation.It may be suspected that a correction is needed for the change in observed heat flow between 0 and 0.2 and between 0.2 and 1.0 My; however inspection of the data (not shown) indicates that the heat flow distribution does not change significantly, and it is a small fraction of that predicted in any case.Thus our calculations assume that observed heat flow is the same over the entire width of the 1 My bin (as it is for all bins).Finally, the "passive" power deficit is ventilated hydrothermal transport driven by heat conducted into an upper crustal aquifer (Lister, 1982), resulting in long-term convective fluid exchange between the crust and oceans.The passive power deficit is the total power deficit after subtraction of the nearaxial deficit:

Results and discussion
Figure 4 shows the main results of our analysis of the heat flow deficit for different cooling models.PDFs of the timeintegrated power deficit are plotted (Fig. 4a-c) for maximum integration times between 1 and 60 My.Since the heat flow sampling is divided into 1 My bins, the PDF for the power deficit at 1 My is essentially the derivative of a single CDF and thus appears rough.Nevertheless, the modes (probability maxima) are well defined for all ages.For comparison, the means and medians of the PDFs are also calculated for each age. the low-probability realizations are skewed to low power estimates.For a different visualization, the mean, median, mode, and half-maximum bounds (HMBs) of the power deficit are also plotted as a function of age up to 100 My in Fig. 4df.This shows that the net deficit plateaus around 50 My for all models, where predicted and observed heat flow converge (Figs. 2 and 5).

Model H13
For model H13, our analysis predicts that the net power deficit is 7.8 TW, with half-maximum realizations falling between 3.8 and 10.1 TW.Compared to most previous studies (e.g., Sclater et al., 1980;Stein and Stein, 1994;Stein et al., 1995;Mottl, 2003;Spinelli and Harris, 2011) this is a low estimate.The estimates of 10.3 and 11.5 TW by Sclater et al. (1980) and Stein and Stein (1994), respectively, are even outside of the HMB for this cooling model.As model H13 is nearly equivalent to GDH1 (Stein and Stein, 1992), the difference is related to the heat flow database, seafloor area-age distribution, and statistical treatment.In addition, our mean (6 TW) and median (6.6 TW) net power deficit estimates are even lower than the mode (Fig. 4a).Although they are within the HMB, they are substantially lower than previous estimates.
It is notable that our estimate of 7.8 TW for H13 is equal to that estimated by Hasterok's (2013b) analysis using his filtered heat flow database.However, this is only coincidental as his filtered database removes seafloor with sediment thickness > 400 m, whereas we do not include such a filter.Moreover, even the more comparable unfiltered database of Hasterok (2013b) incorporates additional corrections not employed here, including (1) a thermal correction for thickening sediment, (2) removal of seafloor area with large igneous provinces, and (3) use of median reference model heat flow q m (we use means).With these additional corrections, Hasterok (2013a) predicted a net power deficit of 6.2 TW.As such, our method predicts a roughly 25 % greater net power deficit.Therefore, if Hasterok's (2013b) database and analytical techniques are preferred, our net deficit estimates should be reduced by 20 %, or slightly less than our median estimates.

Models GH and GHC
The power deficits for thermal plate models GH and GHC are shown in Fig. 4b and c.Model GH predicts a power deficit of 10 TW, with half-maximum realizations falling between 6.1 and 12.3 TW.This value is in good agreement with many previous estimates of the power deficit as discussed above.Inspection of Fig. 5 shows that predicted heat flow for GH and H13 is similar except near the ridge where GH heat flow becomes lower than H13.Thus the high heat flow deficit for model GH originates in active hydrothermal transport on ridge axes.For instance, model GH predicts 4-14 W m −2 for < 0.2 My.
Model GHC has similar properties to GH, except that GHC has an insulating oceanic crust, and all heat transport properties are allowed to vary from their experimental values in order to fit basin-scale geophysical observations (Grose and Afonso, 2013).The predicted net power deficit of model GHC is substantially lower than GH.We estimate 6.6 TW, with half-maximum bounds of 2.9 and 8.8 TW (Fig. 4c).Thus, if normal lithospheric cooling is better modeled by GHC, then the net power deficit is about 35 % lower than expected by model GH.Note, however, that this does not mean that the effect of crustal insulation is a 35 % decrease.This is because, again, the differences between models GH and GHC are both crustal insulation and a re-adjustment of all mineral physics properties to best fit geophysical observations.Since the effective thermal conductivity of model GHC was adjusted to be ∼ 10 % higher than GH (Grose and Afonso, 2013), the real effect of crustal insulation on the heat flow deficit is closer to a 45 % reduction.

Model H13
The spatial distribution of ventilated hydrothermal power predicted by models H13, GH, and GHC are highlighted at the top of their respective panels in Fig. 4a-c and are tabulated in Table 1.For H13, the power deficit at 1 My has a probability maximum at 2.5 TW with HMB between 1.9 and 2.8 TW (approximated after smoothing due to PDF roughness).The cumulative near-axial deficit is 32 % of the total.This may be divided into power on the ridge axis (0-0.2My) at 14 % of the total, and 18 % of the total on ridge flanks (0.2-1.0 My).The remaining deficit of 5.3 TW (68 % of the total) is due to passive ventilated hydrothermal circulation away from ridges.This prediction of 32 % of total ventilated hydrothermal power occurring over < 1 My old crust is only slightly higher than most previous estimates.Stein and Stein (1994) predict a 28 % near-axial deficit, Pelayo et al. (1994) predict about 23 %, and Mottl (2003) predicts 29 % occurring near the axis (Table 1).In addition to the use of an updated heat flow database, our slightly higher value may be attributed to our preference to modes of the power deficit rather than the means (Table 1).Using mean values, we calculate a near-axial power deficit closer to 20 % of the total.This is slightly lower than Spinelli and Harris's (2011) prediction that 25 % of the deficit occurs over near-axial seafloor (using Eq. 1 on a conduction-only model).

Models GH and GHC
As shown in Fig. 4b-c   is extracted near ridge axes (< 1 My).To our knowledge, this is a higher fraction than all previous estimates, including Spinelli and Harris's (2011) near-axial estimate of 40 %.Moreover, about 85 % of the axial deficit (40-45 % of the total) is active circulation on < 0.2 My old seafloor.For model GH, the near-axial deficit (5.2 TW) and the active deficit (4.4 TW) are high, leaving about 4.8 TW of heat removed by passive circulation.On the other hand, because the net power estimate of model GHC is significantly lower than GH, the near-axial power estimate (< 1 My) and the passive regime estimate are both only 3.3 TW.For the near-axial environment, this estimate is in good agreement with some previous investigators (e.g., Stein and Stein, 1994;Mottl, 2003;Spinelli and Harris, 2011; Table 1), but this passive estimate is much lower than previous estimates (Table 1).Some passive power estimates are more than twice our value (e.g., Stein and Stein, 1994;Pelayo et al., 1994;Mottl, 2003).The low passive power estimate for GHC originates primarily in the compounding effects of thermal rebound from active hydrothermal circulation and thermal insulation of oceanic crust.The hydrothermal model from Spinelli and Harris (2011) predicted a passive power budget of 5.4 TW, about 60 % higher than our estimate for model GHC (Table 1).Hasterok's (2013b) passive estimate of 3.9 TW is clos-est to our result.However, this value is for his unfiltered database which predicted a total deficit of 6.2 TW (Table 1).
If we used similar methods, we estimate that predictions for model GHC would be about 20 % lower, with a passive deficit of ∼ 2.6 TW and total deficit of ∼ 5.3 TW (80 % of 6.6 TW).Thus, while the purpose of this work is to demonstrate the relative importance of crustal properties on seafloor heat flow, methodological assumptions in data analysis are important for absolute estimates of ventilated hydrothermal power.

Heat flow constrained by topography
Our analysis has shown that insulating oceanic crust and hydrothermal circulation jointly impact estimates of hydrothermally mined energy in oceanic lithosphere, as well as its spatial distribution.This occurs because both effects result in the prediction of significantly lower heat flux over young (< 30 My old) seafloor compared to conventional models (Fig. 5), and active hydrothermal circulation elevates net heat flow on ridge axes.
Clearly, an important question is whether or not geophysical observations actually support lower heat flow over young seafloor.Although measured seafloor heat flow is contaminated by ventilated hydrothermal circulation processes resulting in a sub-conductive heat flow (Fig. 1), other geophysical observations such as seafloor subsidence are thought to be robust alternative constraints on lithospheric heat loss (Parsons and McKenzie, 1978;Sandwell and Poehls, 1980;Wei and Sandwell, 2006).Hofmeister and Criss (2005), based on the assumption that classical lithospheric cooling models do not fit the Earth, and that hydrothermal circulation cannot lead to what we have called a sub-conductive heat flow (Fig. 1), suggested that true heat loss from oceanic lithosphere is that actually measured (∼ 20 vs. ∼ 30 TW; Von Herzen et al., 2005).Wei and Sandwell (2006) attempted to show that, although measured heat flow does not match lithospheric heat loss as predicted by plate models, it is reflected in the rate of seafloor subsidence.As seafloor subsidence is related to the change in the integrated temperature of the upper mantle, subsidence may be linked to heat loss.Wei and Sandwell (2006) calculated seafloor heat flow based on a spatial integration of heat flows estimated by local subsidence rates found with the global seafloor depth grid of Smith and Sandwell (1997) and age grid of Müller et al. (1997).They calculated local seafloor heat flow by employing the equation where ρ is the lithospheric density, C p is the lithospheric specific heat; α is the lithospheric thermal expansivity; t (x, y) is the age as a function of spatial coordinates; q a is an additional heat flux which is extracted from sub-lithospheric mantle rather than the lithosphere (e.g., the adiabat); w is seafloor depth; and is the isostatic correction for seawater load, with mantle density ρ b and seawater density ρ w .
Although they found that their heat flux calculations were in good agreement with a simple half-space cooling model (where g = 480 mW m −2 My 1/2 ) and conventional net seafloor heat flux estimates, this agreement only occurred upon addition of a hidden lithospheric heat flux, q a , of 38 mW m −2 , or about 11 TW.Consequently, by removing this additional heat flux over global seafloor, we see that the analysis of Wei and Sandwell (2006) suggested that an empirical subsidence-based estimate of net seafloor heat flux is actually on the order of 20 TW, in agreement with Hofmeister and Criss (2005).
We suggest that the solution to this dilemma lay in the age dependence of the effective thermal properties of the lithosphere, especially the consequence of thermally insulating oceanic crust (Appendix B).For example, Goutorbe (2010) and Grose (2012) found that simple thermal plate models with temperature-dependent thermal properties (i.e., no dependence on age), optimally fitted to geophysical observations, required an effective thermal expansivity about 30-40 % lower than the experimental value for olivine.Since Wei and Sandwell (2006) performed no ad hoc adjustments to thermal properties, their use of the high experimental thermal expansivity for forsterite resulted in a low seafloor heat flux, and an 11 TW addition of heat was necessary for a reasonable result.On the other hand, our model GHC has two characteristics which result in different predictions.Firstly, the mineral physics thermal expansivity only needed to be reduced by 15 % to fit geophysical observations for model GHC.Secondly, the effective thermal expansivity is significantly lower near ridge axes compared to old seafloor since the properties of the lithosphere are crust dominated, rather than mantle dominated, at old ages.This second item is relevant because it means that seafloor heat flow at old age represents the cooling of lithosphere having a much higher effective thermal expansivity than over young lithosphere.Consequently, a given rate of subsidence over old seafloor indicates higher heat flux than the same subsidence rate over young seafloor.Fitting a model which has different properties over young versus old lithosphere results in a different view of general lithospheric behavior.
Consider that, similar to Eq. ( 7), seafloor heat flow and subsidence may be related by the equation (Appendix C) where b ei = dw/d √ t is the transient subsidence rate, with w the seafloor depth, and α ei the transient effective thermal expansivity.The age dependence in the rate of heat flow diminution g(t) is shown for our models GHC and GH, and Hasterok's (2013a) model H13, in Fig. 5b.While the heat flow predictions of model GHC are not as low as that measured, as expected by Hofmeister and Criss (2005), the difference between model and measurement is significantly lower than that estimated with simple plate models.Since this decrease is only partially compensated by high heat flow on ridge axes, the power transported by hydrothermal circulation is also markedly reduced.Accordingly, the problem of the heat flow deficit can be attributed to both complex thermal properties and ventilated heat loss.

The subsidence rate
Age dependence of the subsidence rate -as predicted by models GDH1 (Stein and Stein, 1993), H13 (Hasterok, 2013a) and our models GH and GHC -are shown in Fig. 6a compared to empirical estimates based on the global database of Hillier (2010).The empirical estimates are obtained from fitting a line (least-squares fit) through the data (0.1 My bins) in a sliding window of width δt 1/2 .Figure 7 shows the predictions of several lithospheric models compared to the data of Hillier (2010) as well as ridge flank topography of the  (Grose and Afonso, 2013), the classical model GDH1 (Stein and Stein, 1992), and the new heat-flow-constrained plate model of Hasterok (2013b), compared with subsidence rates estimated from the global depth data set of Hillier (2010).Red-yellow colors correspond to small sliding windows, and blue-violet colors correspond to large sliding windows over which subsidence rates are determined using a least-squares fit.The small discontinuities around 20 My for GHD1 and H13 are due to the imperfect fit of their author's respective equations for seafloor depth.(b) Comparison of model predictions of model GHC, GC, and G. Models GC and G have not been used to calculate hydrothermal power loss but are included here to clarify the role of hydrothermal circulation and crustal insulation on seafloor subsidence.
East Pacific Rise from Cochran and Buck (2001).To explore the sensitivity of estimates to the sampling window size, we show estimates for δt 1/2 between 0.2 and 2.0 My 1/2 .Due to the roughness of the data in small bins, the variance is high when the window is small and decreases as the window becomes larger.
Comparison of model predictions and empirical estimates shows that models GH and GHC both fit the general trends in almost all of the data well, while models H13 and GDH1 substantially overpredict subsidence rates for the youngest lithosphere (< 5 My).Model H13 does not fit the data well because it was fit to the depth curve predicted by the model of McKenzie et al. (2005), but we include it for completeness.The empirical subsidence rate clearly has a rising trend between near-zero age and about 30 My, and then decreases gradually in accordance with seafloor flattening.It is notable that empirical subsidence rates using large δt 1/2 tend to rise near zero age and thus appear to be in better agreement with model GH rather than GHC.However, this is due to the loss of resolution due to binning.Model GHC is a superior fit to ridge flank subsidence.The estimates with δt 1/2 ∼ 0.2 My −1/2 indicate that the subsidence rate is ∼ 150 m My −1/2 near the ridge axis (∼ 0.5 My), increases to ∼ 350 m My −1/2 around 30 My, and finally decreases gradually from the effects of seafloor "flattening" to great age.To explain the contributions of axial hydrothermal circulation and crustal insulation, we compare subsidence rates from three models in Fig. 6b.Along with model GHC, we include previously unmentioned models G and GC.As indicated by their names, model G has neither crustal insulation nor hydrothermal circulation, and model GC has crustal insulation but no axial hydrothermal circulation.Subsidence rates for model G are similar to GDH1 (Stein and Stein, 1992) over young seafloor, with a roughly constant subsidence rate of ∼ 340 m My −1/2 .Model GC decreases the subsidence rate to about ∼ 270 m My −1/2 near ridge axes, rising gradually over about 20-30 Ma until flattening.Model GHC further decreases the subsidence rate near the ridge axis to about 110 m My −1/2 .An important observation to consider is that, while the effect of axial hydrothermal circulation on the ridge flank subsidence rate is markedly greater than that for crustal insulation, the hydrothermal circulation effect is confined to the youngest ages while crustal insulation is more persistent.
An important question is whether or not the decreasing subsidence rate in proximity to ridge axes, and the corresponding misfit of model GDH1, is a real reflection of isostatic balance, or if other contributions, such as flexural effects, are important.If other processes are important, this could indicate that observation and model prediction of a low subsidence rate on ridge flanks is not actually related to crustal insulation and hydrothermal circulation.Cochran (1979) showed that gravity anomalies are present over ridge axes.However, while the anomalies appear somewhat significant over Atlantic ridges, they are small and confined to the immediate vicinity of ridge axes over the EPR.Consequently, while elastic sources of deviation from isostasy may be present on the ridge axis, to our knowledge there is no compelling evidence to believe that non-isostatic effects are responsible for the low subsidence rate on ridge flanks (> 0.2 My).
Figure 7b shows a close-up of predicted topography for several lithospheric models compared to the data of Hillier (2010) and EPR data from Cochran and Buck (2001).
Models with no crustal insulation or axial hydrothermal circulation do not fit the low subsidence rates of ridge flanks, nor the axial rise for ages < 0.5 My 1/2 .Models with only hydrothermal circulation or only crustal insulation improve the fit near ridge axes but still fail to provide a satisfactory fit to the data.On the other hand, model GHC provides a satisfactory fit to seafloor topography.A remaining misfit occurs at zero age, where axis depth of model GHC is about 150 m greater than that in the data of Cochran and Buck (2001).This can be explained by the buoyancy of melt.The depth change of the seafloor due to solidification of melt, with isostatic correction and using an exponential solidification rate, can be given as where n = 0.04 My is a time coefficient for the solidification rate, ρ m = 2750 kg m −3 (Stolper and Walker, 1980;Stolper et al., 1981;Hooft and Detrick, 1993) is the density of melt, ρ s = 2950 kg m −3 is the density of solidified melt (ocean crust), φ m is the melt volume, and other terms are defined previously.The volume of melt beneath ridges (depth-integrated φ m in Eq. 10) can be estimated based on the seismic tomography study of a segment of the EPR by Dunn et al. (2000).Their seismic inversion models are consistent with a depth-integrated melt column of either ∼ 0.8 or 1.9 km, depending on whether the pore texture consists of thin films or spheres, respectively.From Eq. ( 10) this predicts w m = −80 and −190 m, respectively, at t = 0, in good agreement with the difference with model GHC (Fig. 2b).
w m predicted by Eq. ( 10) are shown in Fig. 7b.Although Cochran and Buck (2001) argued that the axial rise and low subsidence rate over ridge flanks can be explained by axial hydrothermal circulation, our models indicate that axial circulation alone may not be able to explain these features.On the other hand, all major topographic features characteristic of normal fast-spreading seafloor are explained by model GHC.
It may be suggested that model GH can be made to better fit ridge flank topography by increasing the rate of hydrothermal heat removal (e.g., Nu > 10).This may improve the topographic fit, although we have already noted (Sect.2.4) that a consequence of this (as Cochran andBuck, 2001, andSpinelli andHarris, 2010, used Nu = 20) is a significantly colder near-axial environment which Grose and Afonso (2013) showed is not consistent with evidence from seismic models of the EPR by Dunn et al. (2000).Moreover, in Sect. 5 we will show that model GH (let alone such a model with additional axial hydrothermal cooling power) is not consistent with estimates of heat available for release in the process of crustal accretion.

Note on old-age topography and seafloor flattening
Although the old-age behavior of the lithosphere does not affect our calculations of ventilated hydrothermal power extraction since the heat flow deficit only exists over young seafloor, the fit of models to old-age seafloor is a major constraint on general lithospheric properties.Consequently, the better fit of model GDH1 to the old-age (>50 My) subsidence rate in Figs. 6 and 7a might seem to indicate the superior explanatory powers of GDH1, or at least a problem with the properties of model GHC.However, if the model of lithospheric cooling is meant to represent the normal cooling behavior of the lithosphere, this superior fit may be treated as a flaw rather than a success.Following Crosby et al. (2006) and Hillier and Watts (2005), we consider it likely that seafloor around 100-130 Ma (10-11.5 Ma 1/2 ) is anomalous.Of particular note, the subsidence rate becomes negative around these ages, which cannot be accomplished by passive cooling processes and is not known to occur by means of small-scale convection beneath old seafloor (Zlotnik et al., 2008;Afonso et al., 2008).Thus, we consider models GH and GHC to better reflect the normal behavior of oceanic lithosphere.

Constraints from the thermal budget of crustal cooling
Based only on the fit to the global data, we cannot discount the possibility that crustal insulation is an unimportant contribution.On the other hand, we may cast doubt on the feasibility of model GH for a thermodynamic reason, thus requiring a contribution from crustal insulation.It is notable that the estimate of near-axial circulation for model GH is so high (5.2TW, Table 1).Based on the limited heat budget for ocean crust formation and cooling, Mottl (2003) argued that axial cooling cannot be more than about 3.1 TW.This indicates that model GH, due to high effective thermal conductivity and axial boundary conditions around the shallow axial magma chamber, may be extracting more heat than can realistically be released by advection and crystallization of magmas.Model GHC, on the other hand, only transports about 3.3 TW of heat, consistent with Mottl's arguments.In addition, we previously showed (Grose and Afonso, 2013) that model GHC predicts ridge thermal structure that is in good agreement with a seismic model over the East Pacific Rise by Dunn et al. (2000).As this fit to a seismic model reflects the Nu number, the seismic model supports the choice of Nu ≈ 10 rather than significantly higher values suggested elsewhere (e.g., Cochran and Buck, 2001;Spinelli and Harris, 2011).On the other hand, Han et al. (2014), based on the observation of off-axis magma lenses in regions Dunn et al. (2000) expected to be cool, suggested that Dunn's model may be inaccurate.However, it is not clear what ambient thermal structure is consistent with the presence of off-axis magma lenses, as these may be anomalous, even if frequent, features.We stress that models GH and GHC use a simple Nu-number approximation of hydrothermal transport and therefore can only represent the average behavior of ridges both along and across axes.Moreover, if the thermal structure of Dunn et al. (2000) is too cold, the necessary adjustments to model GHC may be small, except in direct proximity to the axial magma lens.Such corrections may also be applied to model GH, but they will be larger and may be at the cost of good fit to seafloor subsidence.
In summary, if there is no effect of crustal insulation, the fit to seafloor subsidence is slightly compromised and unrealistic amounts of heat are extracted on axis.Thus, we suggest that both insulating oceanic crust and a moderate amount of axial hydrothermal circulation are important for estimates of ventilated hydrothermal power.

Constraints from high-resolution heat flow surveys
A primary goal of this work is to show that crustal insulation strongly affects lithospheric cooling as well as the amount and distribution of hydrothermal power loss.Precise constraints on any of these three items are impossible without knowledge of the conductive lithospheric heat loss over young seafloor.Hasterok et al. (2011) and Hasterok (2013b) attempted to use the global heat flow data set to constrain deep lithospheric heat flow over young seafloor by applying special filters to the data.These authors showed that, by removing seafloor regions with thin sediment, the residual heat flow over young seafloor markedly increased.Because extensive thick sediment cover acts as an impermeable boundary over the crust, fluid exchange between the crust and oceans should become less effective with sediment thickness and heat transfer by conduction becomes dominant (Lister, 1972).Specifically, by removing oceanic regions with sediment cover > 400 m, measured heat flow is in good agreement with model H13 (and GDH1) for ages > 25 My.However, for younger ages a significant deficit remained (Fig. 5).Interestingly, while the deficit between filtered measurements and model H13 is large, the deficit with model GHC is substantially less pronounced (Fig. 5).As the main result of our paper shows, this implies a lower net power deficit.Nevertheless, the remaining deficit indicates that a filtered global database still fails to remove the effects of ventilated hydrothermal circulation on the seafloor heat flow.Consequently, the conductive heat flow still cannot be constrained over young seafloor with global heat flow data.
An alternative approach is to examine specific regions which have been studied extensively enough to demonstrate that the effects of ventilated hydrothermal circulation are small in certain areas.Among the global data, Hasterok et al. (2011) recognized four sites on young seafloor which have thick sediment cover and have been extensively surveyed.These include the Juan de Fuca Ridge flank (Davis et 1997,1999), the Costa Rica Rift (CRR) flank (Davis et al., 2004;Hobart et al., 1985;Langseth et al., 1988), the Gulf of Aden (Cochran, 1981;Lucazeau et al., 2008Lucazeau et al., , 2010)), and the Cocos Plate (Hutnak et al., 2008).Our analysis of the heat flow data for these sites is discussed in the following section, and the resulting estimates are shown in Fig. 5. Subsequently, we discuss interpretations of the data in relation to the different estimates of young-age lithospheric heat flow by models H13, GH, and GHC.
6.1 Heat flow at high-resolution sites

Gulf of Aden
The Gulf of Aden is a rifted margin between Africa and the Arabian Plate which separated at 34 Ma with the onset of seafloor spreading around 18 Ma (Leroy et al., 2012).Lucazeau et al. (2010) reported high-quality heat flow measurements along eight seismic profiles near the margin of the Arabian Plate (near Dhofar), seven of which are aligned with the direction of spreading.They correct their heat flow measurements for sedimentation rate (we use the average of their two cases), topography, and heat refraction.These profiles extend from the continental domain, through the ocean-continent transition, and onto oceanic lithosphere.The small variance in heat flow along profiles led Lucazeau et al. (2008) to conclude that effects of hydrothermal circulation are not important.Of their 162 measurements along 8 profiles, we use 40 points from 6 profiles located on seafloor with ages known from magnetic anomalies (d 'Acremont et al., 2010).Comparison of the model age grid of Müller et al. (2008) and magnetic anomaly isochrons suggests that ages are overestimated by 5-10 My near the continental margin.Therefore, we neglect two profiles from seafloor on the east side of the Socotra Hadbeen fracture zone, since we cannot confidently determine precise ages.The 40 measurements are on seafloor 16-17.6My old and are plotted in Fig. 5 along with their mean and standard deviation (114 ± 10 mW m −2 ).It is likely that this survey examines anomalous seafloor.The site characterizes the early stages of rifting margins, the onset of seafloor spreading, and any thermal consequences of abutting a continental margin.Moreover, based on an examination of heat flow and thermomechanical modeling, Lucazeau et al. (2008) suggested that an intense (300 • C) thermal anomaly below the ocean-continent transition may be likely.If this is the case, then the reported heat flow values may be elevated above that of normal seafloor.

Cocos Plate
Hutnak et al. ( 2008) performed a regional survey of heat flow over Cocos Plate seafloor with ages 18-24 My.The region is blanketed with thick (400-500 m) sediments except for unevenly spaced sites of outcropping basement.Distributed throughout the region are colocated heat flow and seismic-reflection profiles, some of which extend from outcrop sites and others are interspersed about.This heat flow survey revealed a bimodal areal variation in surface heat flow for low (∼ 30 mW m −2 ) and high (∼ 110 mW m −2 ) heat flow areas, a pattern which the authors explain by low-temperature hydrothermal discharge and recharge among outcrops in the low-heat-flow areas, and "warm" hydrothermally inactive crust for the high-heat-flow areas.Their estimate of 97-120 mW m −2 is plotted in Fig. 5 with the age range 18-24 My.

Juan de Fuca Ridge
The Endeavour segment of the Juan de Fuca Plate is heavily sedimented and has been extensively studied with colocated heat flow (Davis et al., 1997(Davis et al., , 1999)), seismic reflection profiles (Rosenberger et al., 2000), and geochemical study from nine Ocean Drilling Project (ODP) boreholes (Elderfield et al., 1999) over a 80 km transect in the direction of spreading.Heat flow measurements are shown in Fig. 5.Note that measurements cluster around 1 and 3.6 My.The scattered black line is an empirical calculation based on the relationship between sediment thickness, basement temperature, and surface heat flow (Davis et al., 1999).Davis et al. (1999) estimated that a sediment correction of +15 % was necessary to estimate basement heat flow.However, we use the +6 % correction from Pribnow et al. (2000) which accounts for thermal anisotropy of sediment.As the data are scattered, we calculate distance-weighted averages of heat flow from the measurements for ages (A) 1.0-1.56My and (B) 3.34-3.6My.

Costa Rica Rift
The geophysical environment of ∼ 6.5 My old seafloor around ODP Hole 504B on the CRR has been characterized with a high-resolution (∼ 1 km spacing) gridded survey of heat flow (Davis et al., 2004;Hobart et al., 1985;Langseth et al., 1988) and seismic reflection profiles (Swift et al., 1998).Figure 8a shows a sediment thickness map produced from the seismic reflection profiles (12 545 points) and a bicubic spline.Figure 8b shows a bicubic spline of heat flow.Inspection of Fig. 8a and b shows that there is some correlation between sediment thickness and heat flow, although the relationship is rough (Swift et al., 1998).Using only the measurements, heat flow for the region has a mean of 229 ± 46 (1σ uncertainty) and median of 218 (194,250) mW m −2 (interquartile range uncertainty).However, since measurement coverage has a greater density around sites of elevated heat flow and thin sediment cover, the statistics are biased to elevated values.The statistics for the bicubic spline are shown as a probability density function in Fig. 8c; they are also compared to a coarser PDF for the measurements alone.The mean of the spline is 211 ± 35 mW m −2 , the median is 203 (186, 229) mW m −2 , and the mode is 190 mW m −2 (173, 217) (half-maximum uncertainty).If we consider that sampling coverage is extensive enough to cover spatial heterogeneity, the mean of the spline may be preferred.However, we note that there are many areas where closely spaced points reveal exceptional lateral gradients, even away from areas of thin sediment or evidence of anomalies in basement topography (Fig. 8).Moreover, recall that measurements on the Juan de Fuca flank have a much smaller spacing (∼ 250 m) and show substantial scatter.A similar scatter may be normal for the CRR site, so that the mean of the spline may not accurately characterize the true mean.We thus take the mode of the spline as the lowest reasonable statistical tendency.

Gulf of Aden and Cocos Plate
The goal of the examination of high-resolution heat flow surveys is to evaluate whether or not the deeper (sub-crustal) lithospheric heat flux can be constrained by close regional inspection.Heat flow surveys of the Cocos Plate and Gulf of Aden regions are in best agreement with model GHC, although the statistical bounds (standard deviations or halfmaximum bounds) are spread over all models.Thus, while these surveys indicate that the conductive models are a better indication of lithospheric heat flow compared to the global heat flow data (with or without a filter for sediment thickness), predicted heat flows for models over seafloor aged 16-24 My are too similar to clearly differentiate models.On the  8).The black line is model GHC with oceanic crust varied between 0 and 10 km.The black dashed line is the same model, except hydrothermal circulation is not allowed to occur below the insulating layer.The red line is model GHC with the thickness of oceanic crust varied and no hydrothermal circulation on the ridge axis.All predictions are for 5.9 My old seafloor.
other hand, because the Gulf of Aden might be anomalously warm (Lucazeau et al., 2008), normal lithospheric heat flow for this age may be somewhat lower, which can improve the agreement with model GHC but worsen the agreement with H13 and GH.

Juan de Fuca
Comparison of the Juan de Fuca (Endeavour segment) heat flow to the models generally suggests best agreement with model GH.Ignoring the tail to low heat flow for ages < 1 My (Fig. 5), the young-age end of the data is in best agreement with model GH, while the bin of data around 3.5 My is in between model GHC and GH.The conventional hydrogeological interpretation of heat flow along the Endeavour flank includes recharge over the young (< 0.66 My old) unsedimented seafloor and discharge at a basement high around 1.3 My, about 20 km from the recharge site (Davis et al., 1997(Davis et al., , 1999;;Newman et al., 2011).Borehole measurements and correlations of surface heat flux and sediment column thickness indicate that the temperature of crustal basement is roughly constant over this region.
Elevated heat flow in this region may be explained by (1) spatial sampling bias, (2) a source of additional heat flow, or (3) different properties of the lithosphere than expected by model GHC.Possible sources of elevated heat flow are (1) continued deep hydrothermal circulation, perhaps associated with faulting as suggested by Nedimovic et al. (2009); (2) heat release from hydration of the crust (Lowell and Rona, 2002); (3) heat transported from younger seafloor and discharged, consistent with models (Davis et al., 1997(Davis et al., , 1999;;Newman et al., 2011); or (4) microbial thermogenesis, which is difficult to quantify but could be significant if nutrient supplies are adequate (D. LaRowe, personal communication, 2012).Davis et al. (1989Davis et al. ( , 1997Davis et al. ( , 1999) ) suggested that advective heat loss around 3.5 My old seafloor is large enough to depress regional heat flow below that of basement heat flux.This would mean that heat flow predicted by model GHC is too low for this region.However, if additional heat has been introduced into the upper crustal aquifer by exchange with the deeper crust, perhaps driven by faults in this rough basement region (Davis et al., 1997), then measured heat flow may be greater than that conducted through the deeper crust, not less.If model GHC is correct for this region, measured heat flow suggests that measurements are elevated due to advective heat exchange between the upper and lower (or deeper) crust.In other words, this region may be super-conductive (Fig. 1), albeit with moderate ventilated discharge.

Costa Rica Rift
The mode of 190 mW m −2 (173, 217; half-maximum bounds) calculated for the bicubic spline over the extensive 2-D Costa Rica Rift survey is in good agreement with models GH and H13 but is significantly higher than the ∼ 170 mW m −2 estimate from model GHC at 6 My.At least part of the reason for this is that the CRR does not consist of "normal" seafloor.In this region, the thickness of oceanic crust has been constrained to be ∼ 5 km thick, significantly thinner than other examples of seismically normal oceanic crust (Becker et al., 1989).Thus, the crustal insulation effect in this region should be lower.We have tested this by performing a sensitivity analysis with model GHC wherein the thickness of the insulating layer is varied between 0 and 10 km (Fig. 9).Model GHC with a 5 km thick crust, or GHC-CRR, predicts heat flow ∼ 182 mW m −2 .This is better agreement, although it remains barely within the lower bound of uncertainty using a somewhat generous statistical technique.An additional contribution may be related to the reduced vigor of axial hydrothermal circulation due to a thinner crust.As an end-member case, we test predictions of model GHC with no hydrothermal circulation and variable crust thickness (red line).The resulting model is in good agreement with the mode, although this is a maximum.Lastly, the dashed line is a test of GHC-like models in which the crust thickness is varied and the Moho is treated as a maximum depth limit for axial hydrothermal penetration.Therefore, realistic models for the CRR likely occur between the dashed line and the red line for 5 km thick crust, or about 182-189 mW m −2 .While this is a remarkable agreement with our statistical analysis, the completeness of sediment cover and apparent absence of basement outcrops and ventilated circulation in the region (Davis et al., 2004) may suggest that the median or mean, which is about 200-210 mW m −2 , may precisely capture lithospheric heat flux.Under this interpretation, CRR heat flow remains significantly elevated above the predictions of model GHC.We suggest that this is in fact the case, such that lithospheric heat loss in this region is noticeably higher than predicted by model GHC.Specifically, we suggest that this region is characteristically superconductive due to deep hydrothermal circulation as depicted in Fig. 1. Davis et al. (2004) found that deep borehole thermal gradients decrease from apparently near-conductive values in the upper crust (∼ 200 mW m −2 ) to much lower values (< 100 mW m −2 ) at depths greater than ∼ 700 m from the basement surface.Davis et al. (2004) suggested that this is due to hydrothermal convection in deeper layers of the crust.If hydrothermal redistribution of heat occurs on such a scale in the crust, this will essentially result in the appearance of a higher effective thermal conductivity of the crust by raising its Nu number.Consequently, GHC may predict low heat flow because it does not use a high Nu number in crust older than 0.2 My.

Summary of high-resolution sites
Overall, the four specific sites discussed above show heat flow elevated above the predictions of the sediment-filtered global database of Hasterok (2013b) and are in rough agreement with all models considered.As previously recognized (Davis et al., 1997(Davis et al., , 1999;;Lucazeau et al., 2008;Hutnak et al., 2008;Hasterok et al., 2011), this demonstrates that the low scattered heat flow over young seafloor is due to ventilated hydrothermal circulation, and that careful geophysical characterization can allow the lithospheric heat budget to be at least partially revealed on young seafloor.Heat flow estimates from the Cocos Plate (Hutnak et al., 2008) and the Gulf of Aden (Lucazeau et al., 2008(Lucazeau et al., , 2010) ) are in marginally better agreement with model GHC.However, these sites are located on seafloor where predicted heat flow from all models is not significantly different (Fig. 5), and the Gulf of Aden might not be considered normal seafloor.The Juan de Fuca Ridge flank is geophysically well characterized, but there is strong evidence of ventilated discharge over much of the sampled area (Davis et al., 1997(Davis et al., , 1999)).The elevated heat flow in this area may reflect this discharge, persistent deep hydrothermal circulation along faults (Nedimovic et al., 2009), heat release from hydration (Lowell and Rona, 2002), or advection from younger seafloor.Heat flow on the Costa Rica Rift is probably the most important data point as it (1) samples seafloor young enough to potentially differentiate models, (2) is heavily sedimented with no evidence of thermally significant ventilated transport, and (3) is well characterized for heat flow and basement depth.CRR heat flow is higher than our preferred model GHC, which may be partly explained by thin oceanic crust.Probably, however, the mean or median of heat flow in this region, which is significantly higher than predicted by model GHC, is a good indicator of lithospheric conduction.The high heat flow in this region may be attributed to deeper hydrothermal transport as suggested by Davis et al. (2004).The Costa Rica Rift is in superconductive state as depicted in Fig. 1.An important unanswered question then regards the "normality" of such deeper transport in the crust.Such transport may be regarded as a process of passive hydrothermal circulation (albeit sealed from exchange with oceans) which is not included in model GHC.If this is normal for global oceanic crust, then model GHC may ultimately overestimate the effect of crustal insulation on lithospheric cooling.
Nevertheless, because only two sites considered here are located on < 10 My old lithosphere, and only the Costa Rica Rift appears to provide a precise constraint on lithospheric heat flow, it is difficult to judge the validity of these models using any existing heat flow measurements.Ideally, detailed surveys such as that performed near the Costa Rica Rift should be performed over several other well-sedimented near-ridge regions globally.Even still, it is not clear if heat flow measurements will succeed in providing precise constraints on normal lithospheric heat loss near ridge axes.We thus expect that continued development of our understanding of the relationship between heat loss and seafloor topography may be instrumental (e.g., Parsons and McKenzie, 1978;Sandwell and Poehls, 1980;Wei and Sandwell, 2006).Since heat flux may be calculated if the effective thermal expansivity and volume heat capacity are known (Eq.9), future efforts should focus on linking detailed models of near-ridge environments (e.g., Cherkaoui et al., 2003;Maclennan, 2008;Craft and Lowell, 2009;Theissen-Krah et al., 2011) with comprehensive mineral physics and hydrogeological models of the crust and lithosphere (e.g., Davis et al., 2004;Afonso et al., 2007Afonso et al., , 2008;;Hasterok, 2010;Goutorbe and Hillier, 2013;Grose and Afonso, 2013).

Conclusions
We have estimated the power of ventilated hydrothermal heat transport, and its spatial distribution, using a set of recent plate models which highlight the effects of hydrothermal circulation and crustal insulation.The most important conclusion of our study is that a model with both of these effects predicts that the difference between measured and modeled heat flow is significantly lower than most previous estimates.Consequently, the total heat vented to the oceans by hydrothermal circulation is lower, and we also find that the fraction of that vented is higher on ridge axes.
Our estimate of ventilated hydrothermal power for a model with constant thermal properties is similar to the recent analysis of Hasterok (2013a), predicting a net power deficit of 7.8 TW, 34 % of which is extracted near the ridge axis (< 1 My).The effect of axial hydrothermal circulation alone is a higher net power deficit (10 TW), and about 50 % of the hydrothermal heat flux occurs near ridge axes (< 1 My).Finally, the effect of crustal insulation with hydrothermal circulation is a markedly lower net power deficit (6.6 TW), with no relative change to the heat flow distribution (50 % near the axis).If median or mean estimates are preferred against probability maxima, total ventilated hydrothermal power estimates are about 10 or 20 % lower than these estimates, respectively.
Many physical and chemical processes in the Earth may be affected by the above predictions.As less heat is transported by hydrothermal circulation, this may also imply that less fluid is circulated in the crust, or that such fluids have a lower average temperature.The lower off-axis advective heat flux also suggests that off-axis "diffusive" hydrothermal circulation is not as vigorous as previously thought.These reduced energy constraints must affect chemical exchanges between the crust and oceans, including the passing of nutrients to sub-seafloor microbial communities and the alteration of oceanic crust.In turn, these effects should impact chemical budgets in the subduction factory and the secular chemical evolution of the mantle.Also, the crustal insulation effect may have broader implications for the thermal evolution of the Earth.Crustal insulation reduces the present-day global heat flux by about 2-3 TW.This reduction will multiply into the past as a warm mantle generates systematically thicker crust and greater insulation.
Finally, the important question remains: does our cooling model with axial hydrothermal circulation and crustal insulation represent the average behavior of oceanic lithosphere?With this question in mind, we have studied model fits to empirical estimates of seafloor subsidence and geophysically well characterized sites of heat flow measurements.While a model with both hydrothermal circulation and crustal in-sulation (GHC) best fits global average seafloor subsidence, site-specific heat flow can be explained with a model which does not include crustal insulation (GH).Additional detailed heat flow surveys of < 10 My old seafloor will be helpful in the choice of best models on the basis of heat flow.The heat budget of cooling oceanic crust, however, is in best agreement with model GHC, as model GH predicts a probably unrealistic extraction of energy on the ridge axis.Accordingly, we find that the cooling regime of the near-axial environment and basin-scale oceanic lithosphere is best explained by models with insulating oceanic crust and axial hydrothermal circulation.for mantle with seafloor ages > 2 My, where k is the thermal conductivity; ρ is the density; C p is the isobaric specific heat; T is temperature; t is time; u is the horizontal velocity of the plate (half-spreading velocity); and x and z are the horizontal and vertical Cartesian coordinates, respectively.The surface boundary is always 0 • C, and the lower boundary uses a constant temperature taken from the geotherm beneath the ridge (zero age).The region of "unchanging temperature" near the ridge axis in Fig. 3 is kept constant, while the mantle outside of this region continuously migrates horizontally at the velocity u.The initial geotherm in this region of unchanging temperature is where T p is the mantle potential temperature (1325 • C), A is a constant adiabat (0.5 • C km −1 ), φ m is a total melt fraction, and S is the latent heat of melting (335 kJ kg −1 ).The melt fraction is estimated based on isentropic batch melting from Asimow et al. (2001).To predict the density of the mantle and crust as a function of pressure, we use the third-order Birch-Murnaghan equation of state (EOS): where P is pressure, K 0 is the bulk modulus at P = 0, K T is the pressure derivative of the isothermal bulk modulus, and is the isothermal volume change.With a thermal correction (Anderson et al., 1992;Anderson and Isaak, 1993;Kumar and Bedi, 1998), the thermal EOS in terms of bulk density is ρ(P , T ) (A5) = ρ(P ) , where δ T = 6 is the Grünneisen parameter, K 0 = 130 GPa, K T = 4.8, and α(T ) = α 0 + α 1 T is the temperaturedependent thermal expansivity (α 0 = 1.64 × 10 −5 and α 1 = 1.32 × 10 −8 for the crust, and α 0 = 2.83 × 10 −5 and α 1 = 0.785 × 10 −8 for the mantle).
The heat transport properties involve constraints on the thermal diffusivity and specific heat, which are strongly temperature-dependent.The specific heat is 1.611-12.479/T 0.5 -1728477/T 3 for the mantle and 1.578-13.11/T 0.5 -763984/T 3 kJ kg −1 K −1 for the crust (Robie and Hemingway, 1995;Berman and Aranovich, 1996).Thermal diffusivity is calculated as  (Pertermann and Hofmeister, 2006;Hofmeister andPertermann 2008, Branlund andHofmeister, 2012).Thermal radiative heat transport is included by arithmetic addition of a radiative thermal conductivity from Hofmeister (2005) given by where h is Planck's constant, d is grain size (1 cm), c is the velocity of light, n (1.7) is the index of refraction, A(ν) is the spectral absorption coefficient, k b is Boltzman's constant, and φ is a coefficient equal to 1 or 0, depending on the values of A(ν).We found the solution with the above parameters and fitted the solution with two Gaussian terms to more easily employ k rad (T ) in our models (Grose and Afonso, 2013).The solution to Eq. (A7) uses spectra of olivine at temperature as explained by Hofmeister (2005).Finally, axial hydrothermal circulation is included in models GH and GHC using a Nusselt (Nu) number approximation (e.g., Cochran and Buck, 2001;Spinelli and Harris, 2011), so that the effective thermal conductivity in crust where T < 800 • C is 10 times the lattice value (plus the radiative contribution).Specifically, this may be stated as where k e (z, x) is the effective thermal conductivity in space and k L (T , P ) = D(T , P )ρ(T , P )C p (T ) is the lattice or "phonon" part of the thermal conductivity computed using functions for thermal diffusivity, density, and specific heat discussed earlier.The maximum temperature of 800 • C for high Nu is used as it corresponds with the limit of the cracking front of hydrothermal fracturing in the lower crust of fast-spreading ridges estimated by Manning et al. (2000).Note that radiative heat transport both in the crust and in the zone of hydrothermal transport because the mean temperature can still be high enough for thermal emission, although we note that the abundance of plagioclase and clinopyroxene in the crust, as well as chemical and microstructural alteration, may affect the role of radiation.In our models hydrothermal circulation is extinguished at 0.2 My.The values of Nu = 10 times the lattice value, and a 0.2 My cutoff of axial circulation was used since these conditions resulted in a model which predicted near-ridge thermal structure consistent with seismic models for the East Pacific Rise (Dunn et al., 2000).For example, Spinelli and Harris (2001), using Nu = 20 and a cutoff of 0.1 My, predicted a cooling rate which is too high to fit the seismic evidence (Grose and Afonso, 2013).
which has the same form as Eq.(C2) except that physical coefficients are also transient effective properties.In comparison, the seafloor heat flux is q sf (t) = k dT (t) dz z=0 , (C5) which, for half-space cooling, may be related to properties of the lithosphere as q sf (t) = [ρC p ] ei (t) T ei (t) D ei (t) π t + q a , (C6) where q a is the adiabatic part of the heat flux.Moving time to the left-hand side, we can define the variable which may be referred to as a diminution rate for surface heat flow, analogous to the transient subsidence rate b ei (t).
From the above it can be seen that the seafloor heat flux is not easily related to net seafloor subsidence or even the net subsidence rate.However, the heat flow diminution rate may be related directly to the transient subsidence rate as b ei (t) 2 (t)α ei (t) = g(t) [ρC p ] ei (t) , (C8) which has a similar form to previous derivations of the relationship between the heat content of the lithosphere and seafloor topography (e.g., Parsons and McKenzie, 1978;Sandwell and Poehls, 1980;Wei and Sandwell, 2006).

Figure 1 .
Figure 1. Figure showing the relationship between hydrothermal circulation patterns in the crust and resulting heat flow measured at the surface.Ventilation results in a net reduction in seafloor heat flux (sub-conductive heat flow), and vigorous hydrothermal circulation without ventilation results in a higher measured heat flux (super-conductive heat flow).The black squiggly lines indicate conductive heat transport, and the blue convection lines indicate hydrothermal patterns.The sites of recharge and discharge are located at basement highs in these examples.

Figure 2 .
Figure 2. (a) Cumulative distribution functions of heat flow for 1 My bins to 100 My compared to the three thermal plate models used in this study.(b) Cumulative distribution functions of heat flow in 1 My bins between 1 and 100 My.

Figure 3 .
Figure 3. Boundary and initial conditions for models GH and GHC.The temperature in the region outlined in white dashes varies with depth but remains constant over horizontal distance.The distance between the surface and the top of the axial boundary condition is 1.4 km.All material around this space migrates at a constant rate away from the ridge axis.The top figure shows the boundary conditions of the 2-D part of the model (0-2 My), and the bottom figure shows the conditions of the 1-D part of the model (2-170 My).After Fig. 2 in Grose and Afonso (2013).

Figure 4 .
Figure 4. Probability density functions of the heat flow deficit determined from Monte Carlo analysis for models (a) H13 (Hasterok 2013a), (b) GH, and (c) GHC(Grose and Afonso, 2013).The red PDF (left panels) represents numerical integration only within the first 1 My bin (near-axial power deficit), and the blue PDF is the heat flow deficit integrated to 60 My (net power deficit).The filled circles, triangles, and squares indicate the mode, median, and mean of the PDFs, respectively.The white filled circles with error bars indicate the mode and half-maximum bounds for the 1 My (near-axial power deficit) and 60 My (net power deficit) PDFs.The red, green, and blue bars above the graphs (left panels) indicate the fraction of active, flank, and passive advective power, respectively.Hasterok's (2013b) "unfiltered" estimate is indicated in panel (a).The right-side panels show the mean, median, mode, and half-maximum uncertainty of the power deficit as a function of age for each model.

Figure 5 .
Figure 5. Predicted heat flow (a) and heat flow diminution rate (b) for models H13, GH, and GHC as a function of age compared with global and site-specific data.The site-specific data points are discussed in the text.The global sediment filtered and unfiltered data shown here are fromHasterok (2013b).The highlighted region is the uncertainty (1σ ) of the sediment-filtered global data set.Note that the misfit between models GHC and the data occurs becauseGrose and Afonso (2013) fitted model GHC to heat flow data which did not include a correction for continuous sedimentation.

Figure 6 .
Figure 6.(a)Subsidence rates for models GH and GHC(Grose and Afonso, 2013), the classical model GDH1(Stein and Stein, 1992), and the new heat-flow-constrained plate model ofHasterok (2013b), compared with subsidence rates estimated from the global depth data set ofHillier (2010).Red-yellow colors correspond to small sliding windows, and blue-violet colors correspond to large sliding windows over which subsidence rates are determined using a least-squares fit.The small discontinuities around 20 My for GHD1 and H13 are due to the imperfect fit of their author's respective equations for seafloor depth.(b) Comparison of model predictions of model GHC, GC, and G. Models GC and G have not been used to calculate hydrothermal power loss but are included here to clarify the role of hydrothermal circulation and crustal insulation on seafloor subsidence.

Figure 7 .
Figure 7. Predictions and observations of seafloor topography and subsidence on ridge flanks.(a) Seafloor topography of our models compared with model GDH1(Stein and Stein, 1992), 0.1 My binned global "normal" depth data with 1σ standard deviation fromHillier (2010), and EPR ridge depth data fromCochran and Buck (2001).(b) Close-up on the model predictions and observed topography on ridge flanks.In addition, the bottom left corner shows an additional cause of ridge elevation due to the buoyancy of melt beneath the ridge based on two models of seismic tomography fromDunn et al. (2000).For ages < 0.2 My 1/2 there is no subsidence in our models since energy transfer is supported by latent heat release.Figure modified afterGrose and Afonso (2013).

Figure 8 .
Figure 8. Sediment thickness and heat flow for the Costa Rica Rift.(a) Sediment thickness map constructed from a bicubic spline of seismic reflection profiles (dotted line, where dots are data points; Swift et al., 1998).Open circles are sites of heat flow measurement.(b) Heat flow map constructed from a bicubic spline of measurements (black dots).(c) Non-normalized probability density functions (PDFs) of heat flow.The red PDF is from the bicubic spline in panel (c), whereas the black line is the PDF using only measurements.

Figure 9 .
Figure9.Predicted seafloor heat flow for model GHC with varied crust thickness and parameters affecting hydrothermal circulation are varied, compared to measured heat flow (mean, median, mode, and half-maximum bounds) for the Costa Rica Rift (Fig.8).The black line is model GHC with oceanic crust varied between 0 and 10 km.The black dashed line is the same model, except hydrothermal circulation is not allowed to occur below the insulating layer.The red line is model GHC with the thickness of oceanic crust varied and no hydrothermal circulation on the ridge axis.All predictions are for 5.9 My old seafloor.

Table 1 .
and Table1, both models GH and GHC predict that about 50 % of total hydrothermal power Comparison of heat flow deficit estimates for different seafloor domains for this study compared to previous studies.