Fractal density and singularity analysis of heat flow over ocean ridges

Peak heat flow occurs at mid-ocean ridges and decreases with the age of the oceanic lithosphere. Several plate models, including the Parsons and Sclater (PSM) model, Global Depth and Heat (GDH1) model and Constant Heat flow Applied on the Bottom Lithospheric Isotherm (CHABLIS) model, have been used to predict heat flow in the ocean lithosphere. The discrepancy between the predicted and measured heat flow in the younger lithosphere (i.e. younger than 55 Myr) influenced by local hydrothermal circulation has been used to estimate hydrothermal heat flux and investigate hydrothermal processes. We can modify the cooling models by substituting the ordinary mass density of lithosphere by fractal density with singularity. This new model provides a modified solution to fit the observed heat flow data used in other models in the literature throughout the age range. This model significantly improves the results for prediction of heat flow that were obtained using the PSM, GDH1 and CHABLIS models. Furthermore, the heat flow model does not exhibit special characteristics around any particular age of lithosphere. This raises a fundamental question about the existence of a “sealing” age and accordingly the hydrothermal flux estimation based on the cooling models.


Results
Several global datasets have been used in the literature for plate cooling model validation and construction. For example, the PSM model was constrained using about 3300 heat flow determinations 20 , GDH1 by about 5500 data 8 and more recently these models were re-evaluated by Hasterok et al. 28 using an updated global heat flow database including more than 14,000 heat flow determination. In the above works, the cooling models are usually fitted with two functions with age separation at various cut-off values, such as 120 Myr in PSM 20 , 55 Myr in GDH1 model 25 , 80 Myr in CHABLIS model 21 29 . Although these models ensure two functions fitted to the observed data in two age ranges exhibit continuity at the cut-off age, the two separate models are irrational in the sense that the derivatives of the functions with respect to age are discontinuous. Several recent studies used the updated database 27 and subsets by filtering or sediment-correction in model calibration 3,9,[27][28][29][30][31][32] . A few sets of site specific data in young seafloor, such as heat flow measured in bare rocks have also been applied 29,41 . Hasterok 29 used filtered and corrected values from a large database to evaluate several models and suggested a new H model which significantly reduces the issue of hydrothermal deficit on < 25 Myr seafloor. Unfiltered data for old age seafloor, filtered data and several high resolution site-specific heat flow data in young age seafloor were further combined to form a high quality dataset with minimal hydrothermal disturbance for model calibrations 30 .
Here, I utilize the above datasets used in the literature for PSM, GDH1 and H models to evaluate the new power-law model for heat flow derived on the basis of fractal density. The details of the derivation of the model will be described in the Methods section. It will be demonstrated that the power-law model, ct −(1/2+Δα) , where c is a constant and Δ α is the singularity of fractal density of lithosphere, fits the observed heat flow datasets very well throughout the observed age range. The implication of fractal density inclusion into heat flow modelling will be presented in the Discussion section.
The data used by Stein and Stein 8 . Stein and Stein 8 compiled 5,539 data points at sites younger than 166 Myr, with good coverage of the main oceanic plates. Figure 1 shows the averaged heat flow data in 2-Myr bins (as dots) and one standard deviation concerning the mean value. The data was fitted using the PSM and GDH1 plate models. To reanalyse the heat flow data using the power-law model introduced in the current paper, the 83 dots were digitised to be converted to numerical data and were then fitted using the least-squares method. The results include the function q(t) = 151.81t −0.22 , a coefficient of determination R 2 = 0.671, and a Student's t-value = 8.1 (n = 83 samples). The results show that the power-law model with an exponent − 0.22, or Δ α = − 0.28 fits the observed heat flow data very well throughout the age range. Unlike the results obtained by other models such as PSM, GDH1, CHABLIS and H models, the single power-law model fitted to the data over the whole age range does not exhibit special characteristics around any particular age, such as t = 55 Myr. Any order-derivatives of the power-law heat flow function are continuous everywhere except at t = 0, where heat flow exhibits a singularity. 25 . The second dataset is the averaged heat flow data used in Stein and Stein 25 . The data is the same dataset as used by Stein and Stein 8 except the average heat flow calculated in 15 variable unequal age bins including several small bins younger than 6 Myr. Several of the individual values were obtained by averaging the results of detailed surveys on crust younger than 6 Myr to avoid biasing the data. The data is shown in Fig. 2a

The data used by Stein and Stein
The data were fitted using the least-squares method as well as the power-law model. The equation is shown in Fig. 2a. The results (Fig. 2a,b) show that the power-law function with exponent − 0.186, or Δ α = − 0.314, fits the heat-flow data very well throughout the age range, with R 2 = 0.966 and Student's t-value = 13.5.
The data used by Hasterok et al. 28 . The third dataset is from Hasterok 28 who have compiled the most updated and complete heat flow database nearly 15,000 oceanic heat flow measurements. Based on the database, the authors created a new heat flow model which was further used to estimate the conductive heat flow loss through the seafloor. To minimize the effect of ventilated hydrothermal circulation on heat flow and maximize the number of data points used in model fitting, a subset of heat flow measurements was extracted by filtering with > 400 m of sediment cover and located > 60 km from the nearest seamount. In total, 5023 measurements were retained which is about 35% of the initial pre-filtered dataset. The effect of hydrothermal circulation on average heat flow obtained with the new dataset is substantially reduced (Fig. 3a) 29 . However, a small heat flow deficit compared to lithosphere cooling models persists as ages < 50 Myr, especially < 25 Myr. It was interpreted that some of this deficit could result from incomplete thermal rebound and uncertainties in sediment thickness 29 . A combined dataset including unfiltered, filtered, and site-specific data was used to improve goodness of fit to GDH1 model (Fig. 3b) not only for heat flow but also for seafloor depth 29 . In this section, I will use the heat flow datasets 28,29 to evaluate the power-law model.
Both the average heat flow per 2.5 Mys bin calculated from the filtered and unfiltered data 28 are shown in Fig. 3a. The data were fitted using the least-squares method in accordance with the power-law model. The results show that the power-law function fits the average heat flow data very well throughout the age range. The curves fitted to unfiltered data (not shown) and filtered data result in, q(t) = 186.57 t −0.271 with R 2 = 0.81 for unfiltered data, and q(t) = 224.03 t −0.307 with R 2 = 0.91 for filtered data, respectively. Accordingly, the fractal density singularity are calculated as Δ α = − 0.229 and Δ α = − 0.193, respectively. Figure 3b shows the results based on the median heat flow per 2 Myr bin calculated from a combined dataset with filtered, unfiltered and site-specific measurements. Considering the skewed heat flow data distribution and test of importance of choose of non-Gaussian versus Gaussian statistics on the dataset, the median rather than mean statistics were used on assessment of heat flow distribution and estimation of the total power deficit 30 . With the median heat flow data, the least square fit with power-law function is well up to 100 Myr. There would be a slight under-estimation if the power law function is extrapolated to fit the data for whole age range 180 Myr. The median heat flow data from combined unfiltered, filtered and site-specific measurements < 100 Myr were fitted and shown in Fig. 3b, which gives q(t) = 433.63 t −0.445 and Δ α = − 0.055 with R 2 = − 0.97.

Discussions
The three forgoing examples indicate that peak heat flow measured over the oceanic ridges exhibit singular behaviour, following a power law relationship with age or distance from the ridges. The decay trend of heat flow with age can be approximately described by power-law models with negative exponent (− 1/2 − Δ α ). The exponent was estimated ranging from − 0.186 to − 0.307 for average heat flow within 2 or 2.5 Myr bins and − 0.445 calculated for median heat flow per 2 Myr bin. The single power-law model fits the heat flow data in each example throughout the available age range (except for median heat flow data with correction in age range younger than 100 Myr). The results suggest that there are no specific characteristics of measured heat flow at any particular age. The power-law function is continuous, and any order-derivatives of power-law function are also continuous except at age zero at which it exhibits singularity. According to separation of active hydrothermal circulation near the ridge axis and passive circulation away from ridges 25,42,43 the current result may question the presence of a "sealing" age or simply gradual change from a predominantly "active" to "passive" phenomena as the lithosphere age increases. Stein and Stein 25 concluded hydrothermal flow decreases with age due to reduced crustal porosity and permeability. The heat flow fraction depends primarily on crustal age and secondarily on sediment distribution; it increases with age because of reduced flow in the crust due to decreased crustal porosity and hence permeability 25,44 . The overlaying sediment may reduce, but not eliminate the effects of water flow. Therefore, this result indicates that, if it exists, the "sealing" age, representing the age at which the hydrothermal heat flux becomes minor, may not cause abrupt change in heat flow pattern with, for example, discontinuity of derivatives of the heat flow function.
The three examples used for analysis all give negative singularity of fractal density (Δ α ∼ − 0.314 to − 0.193 for average heat flow and − 0.055 for median heat flow) implying significant negative anomalies of fractal density around ridges. The variation of the singularity Δ α around value − 0.25 estimated from these datasets is due to selection and filtering/correction of the heat flow data. It is also affected by the decision of averaging of data with  25 ), the dashed line denotes the fitted curve obtained using the GDH1 model and the solid lines symbolize the fitted curves obtained using the presented model. different scheme of partition the age intervals. For example, in the first example, equal interval of 2.5 Myr was used to calculate average heat flow which gives relatively more data points in old age seafloor than young age seafloor. Therefore the least squares fit would be influenced by the points in old age seafloor. This gives a relatively small singularity Δ α = − 0.28. In contrast, the second dataset uses small number of variable unequal age bins (15 bins) including several small bins younger than 6 Myr, this dataset gives relatively more influence of young data on the least square fitting which can lead to larger value of singularity Δ α = − 0.314. The number of data and data quality in the dataset can also affect the estimation of the model, for example, the third dataset with far more data points gives singularity value Δ α = − 0.229 for unfiltered data and Δ α = − 0.193 for filtered data. The reduced value of singularity estimated from the filtered dataset is understandable since it has few data points from relatively young age which in turn increase the influence of data from relative old age. In general, more restrictive filters applied to heat flow measurements lessen the effect of heat flow data with young age on the heat flow modelling, consequently the weaker the singularity value. This is common in fractal and multifractal modelling for estimating multiple scale singularities that the actual values of singularity could be affected by the range of scale used for calculation. Future study may need to focus on individual ridges to explore the spatial variation of fractal density calculated from heat flow data from different oceans. Ideally, the singularity of fractal density of the crust should be independently estimated by associating average density and scale from measurements about density data of lithosphere. The singularity of fractal density can then be used in the cooling model as parameter. However, to estimate singularity of fractal density of lithosphere over ridges one needs a large number of samples taken from the surface and subsurface to provide database to estimate fractal density considering 3D heterogeneity of lithosphere around ridges. An indirect way to characterize the lithosphere's density might be by inversion of gravity or seismic data since these measurements are related to the density of the lithosphere. These types of data are also affected by other factors. The spatial distributions of oceanic topography and lithosphere density have been found to exhibit multifractality with multiple scale singularities 45 . However, these types of multifractality of lithosphere density have not been considered in heat flow dynamic modelling.
Several authors have discussed the influence of gravity and density of the lithosphere on hydrothermal circulation occurrences in lithosphere over ridges. For example, previous near-bottom gravity data indicated a porosity of 30% over the first few hundred meters of the upper layer of the oceanic crust 41 . High porosity is supported by low seismic velocities 32 . High porosity upper crustal around ridges is also reflected by the low mass density of lithosphere which can be indirectly reflected by negative free-air gravity anomalies. Studies [46][47][48] have demonstrated that the gravity around the oceanic ridges and transform faults are incredibly complex with positive and negative  residual mantle Bouguer gravity anomalies. The mass density of lithospheric over these areas has significant variability due to change of rock porosity, serpentinization of mantle peridotite, faults and fractures, and/or crustal thickening. Watts 46 stated that the most striking of the features found in gravity anomalies over oceanic rifts are large amplitude free-air gravity anomaly low over oceanic rifts in pacific, Indian, and Atlantic oceans. Lambeck et al. 49 reported that gravity anomalies over ridges can be explained by the thermal expansion of the lithosphere at the ridges, suggesting that gravity measurements can be used as constraints on the lithospheric models in the same way as heat flow and topography measurements. To summarize, the above suggestions indicate that the lithosphere around ocean ridges, especially the uppermost crust, is very complex not only in porosity, mass density, composite, topography, seismic and gravity properties, and texture and structures but also the irregular geometry and multi-scale of hydrothermal circulations within the lithosphere. Some but not all of these properties have been accounted in heat flow modelling. It was commented that the deviations of measured heat flow from the plate model prediction may be due to temperature and pressure effects on physical properties/processes not accounted for by the plate model such as lithospheric strength, heat production, and thermal rejuvenation 28 .
Understanding the hydrothermal flow process should ideally reconcile models based on large-scale average crustal properties, such as heat flow or seismic velocity, with site measurement of water flow and crustal physical property 25 . The concept of fractal density proposed in this paper can describe the non-linear scaling property of mass density and characterize the complexity of density with change of support. The new solution derived from the standard conductive cooling model with assumed fractal density may consist of two components: heat flow solution with fractal density ρ α and other constant parameters (such as ∝t −1/2 ) and the component with singularity of fractal density (∝t −Δα ). As explained in the Methods section, the singularity index Δ α has three possibilities: if Δ α = 0, corresponding to ordinary non-singular density, then the solution of cooling model remains unaffected; if Δ α > 0, corresponding to "convex" anomaly of density over ridge, then the solution of cooling model will show an enlarged exponent (0.5 + Δ α > 0.5), in this case the measured heat flow tend to deviate from the ordinary solution; and if Δ α < 0, corresponding to "concave" anomaly of density over ridge, then the new solution of cooling model will show a reduced exponent (0.5 + Δ α < 0.5), in this case the measured heat flow tend to converge to the ordinary solution.
It must be mentioned that the value of singularity of fractal density should be estimated separately from the heat flow measurement so that it can be used to assess the effect of fractal density on the prediction of heat flow by cooling models, and accordingly to judge the influence of hydrothermal circulation on heat flow prediction. Without a clear understanding and good estimate of the singularity of fractal density, it is impossible to differentiate the influence of hydrothermal flow flux and the influence of fractal density. Another complication of this issue is that the singularity of fractal density may vary from ridge to ridge; therefore, fractal density should be investigated in various specific regions around ridges. Considering this, the hydrothermal flux's estimation conducted in the literature based on ordinary cooling model and measurements of heat flow can be questioned. Nevertheless, introduction of fractal density into cooling model for heat flow prediction in ocean ridge might have opened a path for future research. I need to point out that the work about insertion of fractal density in cooling model of the current paper is only introduced for modification of heat flow modelling. Future study should also take into account influence of fractal density on substance of floor which is an important constraint of cooling model.

Methods
Non-linear processes are widespread in nature. Extreme physical, chemical and biological processes occurring in the Earth's crust often result in a singularity of energy release or mass accumulation 50 . In this paper, we show that the end products of these types of singular processes such as hydrothermal systems at mid-ocean ridges can be characterised as fractal densities based on a power-law relationship between density and scale.
Fractal density and power-law model. Since density is a fundamental physical parameter involved in the thermal models (PSM, GDH1 and CHABLIS) used for prediction of heat flow at the mid ocean ridges, we will first introduce a new definition of fractal density with nonlinear property which was ignored by the traditional dynamics systems. The principle of density was discovered by the Greek scientist Archimedes approximately 2000 years ago. Density has become a foundational property of mass and energy and a well-known physical concept with a variety of applications in nearly all fields of study. The density of material or energy is defined as its mass or energy per unit volume. Therefore, density often has units of mass over volume (e.g., g/cm 3 , kg/m 3 ) or energy over volume (J/cm 3 , w/L 3 ).
To calculate the mass density of an object (ρ ), one takes the mass contained in a volume (m(v)) divided by the volume (v): where ρ is the average density of an object, which becomes independent of the volume only if the density of the object is uniform. However, if the object has heterogeneous properties, the density must be calculated using the derivative of the mass over volume: The above density exists only if the limit converges when the volume becomes infinitesimal. In this paper, we will demonstrate that the limit in Eq. (2) does not always converge for complex objects with fractal properties. A new concept of fractal density is defined here as the limit of the following relation (3) if there is a parameter α (a positive value) so that the following limit converges 51,52 : Scientific RepoRts | 6:19167 | DOI: 10.1038/srep19167 The value of ρ α can be considered to be the generalised ordinary density ρ because the relation in Eq. (2) becomes a special case of Eq. (3) when α = 3. The fractal density defined in Eq. (3) has units of the mass ratio to a fractal set of α dimensions, for example, g/cm α or kg/m α . Similarly, the units of fractal energy density can be J/ cm α or w/L α . Combining Eqs. (2) and (3) yields the following relationship between the ordinary density and the fractal density: Thus, the ordinary density obeys a power-law relationship with the volume with the following properties 51,52 : In accordance with these properties, the ordinary density shows a singularity when α ≠ 3. The concept of fractal density and Eqs. (3) and (4) can be extended to other situations, for example, areal and linear densities. A general model associating the fractal density and the ratio of mass and scale (linear size of an E-dimensional set) can be expressed as follows: This equation indicates that both ordinary density and mass follow power-law relations of scale ε. A power-law function is determined by two parameters: the fractal density ρ α and the exponent-singularity index α (fractal dimension), or E-α ; the latter corresponds to the co-dimension of fractal density. The singularity index measures the deviation of the fractal dimension from the dimension of normal density. These two parameters (ρ α and α ) can be estimated from observed data by measuring the intercept and slope of the log-log plot of m against ε, which plots as a straight line. Singularity is observed in non-linear processes including cloud formation 53 , rainfall 54 , hurricanes 55 , flooding 56 , landslides 57 , earthquakes 58 , mineralisation 59,60 , and solar wind turbulence 61 .
Multiplicative cascade processes and formation of fractal density. Multifractal cascade models play a fundamental role in quantifying turbulent intermittency and other non-linear processes 62 . Here, we will use a relatively simple 1-dimensional de Wijs multiplicative cascade model 63  The quantity α (ξ ) represents the singularity index in the context of multifractals 64,65 . Note that when ε n → 0, the ordinary density of bar segments with length ε n and measures μ (ε n ) approaches either infinity or zero unless α = 1, which corresponds to d = 1-d = 0.5. In other words, if d ≠ 0.5, then α ≠ 1, and the ordinary density approaches zero if α < 1 or approaches infinity if α > 1. The dimension of fractal density varies from one location to the next along the bar L. The dimensions of the fractal densities of the two segments at the ends of the bar are α (0) = − log(d)/log(2) ≠ 1 and α (1) = − log(1-d)/log(2) ≠ 1.
Fractal heat density and power-law model. Here we will show that "singular" locations with fractal densities where α < E may indicate an anomalous release of energy caused by hydrothermal events at mid-ocean ridges. According to heat-flow dynamics, the volumetric heat content (H in units of Jm −3 ) can be represented by the equation 2 where c p is the calorific capacity, ρ is the mass density, and T is the temperature. In lithosphere cooling models introduced in the literature for modelling hydrothermal event including heat flow over mid oceanic ridges, the mass density ρ is often assumed to be ordinary density either with constant value of mantle density (3.0 g/cm 3 or 3.3 g/cm 3 ) or treated as linear function of depth or nonlinear function of temperature 66 . Similarly, the values of c p can be taken as constant such as 0.25 cal/g°C or nonlinear function of temperature. The standard analytical lithosphere cooling model ignores the variation of ρ and c p with temperature and treats these parameters as constants. It was pointed out by McKenzie et al. that the effects of such variations of ρ and c p with temperature are small and they are easily included in a numerical scheme 66 . The plate model with temperature dependent parameters fits the variation of heat flow and depth with age at least as good as those of the analytic model with constant and adjustable coefficients. One should be reminded that these models involve ordinary mass density ρ whether as constant value, linear function of depth or nonlinear function of temperature. For example, it can have the normal density unit of g/cm 3 . However, as pointed out in forgoing discussions of this paper, many studies indicate the density variation of lithosphere over ridges is complicated; thus, fractality of density should be considered. Several studies indicated 43-45 that free-air gravity anomaly low over oceanic rifts could reflect that the mass density of lithospheric over these areas are of significant variability due to change of rock porosity, serpentinization of mantle peridotite, and/or crustal thickening. The increase of serpentinization of the uppermost mantle lithosphere can also change the density of ocean crust. The spatial distributions of oceanic topography and lithosphere density have been found of multifractality with multiple scale singularities 44 . The property of multifractals and multiple scale singularities of lithosphere density with respect to measuring scale have not yet been taken into account in any lithosphere cooling models. In this paper, we attempt to consider fractal density of lithosphere (6), the volumetric heat content in a volume of linear size ε should obey the relationship where H α is a fractal thermal density with units of J/m α , and, accordingly, the ordinary thermal density H should contain a component expressed as power-law relationship of scale ε. Given the new definition of fractal density and fractal heat density, the heat flow q expressed in units of volume can be rewritten where k is the thermal diffusivity which has been treated as constant in the standard analytic model and as a function of temperature in numerical modeling 66 , and q α is the heat flow calculated based on the fractal heat density. Therefore, the ordinary heat flow can be treated as the average heat flow per unit volume, obeying the power-law relationship with the volume. This relation (11) indicates that the ordinary heat flow q should be modified by the term tackling the singularity of the fractal density. According to the prevailing paradigm of the PSM, the GDH1 and CHABLIS that model the heat flow measured at the surface of a mid-ocean ridge can be described by a one-dimensional function q(t) where t is the age of the lithosphere. The distance from the ridge x can be calculated from the age of the lithosphere according the constant speeding rate. Accordingly, q(t) can be expressed either as the distance from the ridge or the lithosphere's age. Considering the one-dimensional problem, the scale ε can be treated as two times that of x, ε = 2x. Combining the ordinary heat flow function, q α (t), to be derived from the standard plate model with constant k and c p except that fractal density ρ α is in the form of singularity (11), we can derive the new function as Consequently, the power-law model (12) can be formulated as a function of age with exponent −(0.5 + Δ α ) (Δ α = E − α , E = 1 for one dimensional problem). For CHABLIS model the function (12) becomes approximation when x → 0 and Δ α < 0. This function was fitted to the datasets in Figs 1 -3 and they give estimated value Δ α = − 0.28, − 0.314, − 0.229, − 0.193 and − 0.055, respectively. These values are consistently less than zero, indicating strong negative singularity of mass density of lithosphere would be expected over the ocean ridge.