A Magma Accretion Model for the Formation of Oceanic Lithosphere: Implications for Global Heat Loss

A simple magma accretion model of the oceanic lithosphere is proposed and its implications for understanding the thermal field of oceanic lithosphere examined. The new model (designated VBA) assumes existence of lateral variations in magma accretion rates and temperatures at the boundary zone between the lithosphere and the asthenosphere. Heat flow and bathymetry variations calculated on the basis of the VBA model provide vastly improved fits to respective observational datasets. The improved fits have been achieved for the entire age range and without the need to invoke the ad-hoc hypothesis of large-scale hydrothermal circulation in stable ocean crust. The results suggest that estimates of global heat loss need to be downsized by at least 25%.


Introduction
Detailed understanding of large-scale variations in the thermal field of the oceanic lithosphere provides important constraints on deep tectonic processes.
Nevertheless, thermal models of the lithosphere proposed to date have failed to provide a satisfactory account of some of the important features of large-scale variations in oceanic heat flow. For example, both the Half-Space Cooling [55] and Plate [30] models predict heat flow much higher than the observed values, for young (ages less than 55 Ma) ocean crust. Also, the magnitudes of heat flow anomalies associated with the mid-ocean ridge systems are systematically lower by a factor of 6 at younger ages than those predicted by thermal models proposed in the current literature [41]. In addition, the widths of thermally anomalous zones associated with the spreading centers are narrower (less than 23 Ma) than those calculated (~66 Ma) for a wide range of plausible model parameters. Such discrepancies between model predictions and observational data have given rise to the so-called "oceanic heat flow paradox", for which no satisfactory solution has been found for over the last forty years. The common practice in the current literature is to consider the paradox as originating from eventual perturbing effects of possible regional scale hydrothermal circulation in the ocean crust not accounted for in conventional heat flow measurements (e.g. [41], [54], [56], [27], [60]).
There are however dissenting views on the subject matter of hydrothermal circulation on regional scales [20]. Direct experimental evidences presented to date have confirmed the existence of only isolated pockets of hydrothermal circulation in the central valley and in the rift flanks of spreading centers (e.g. [12], [13], [17], [13], [29], [34]). No direct experimental evidence has so far been presented that point to the existence large scale circulation systems operating in stable ocean crust. Most of the arguments presented to date in favor of the supposed existence of regional-scale convection systems, in parts of ocean floor at large distances from spreading centers, are based on indirect inferences (e.g. [56], [2], [3], [14]).
In the current work, we present a new model of oceanic lithosphere that can overcome the above-mentioned problems and present a satisfactory solution for the heat flow paradox, without the need to invoke the ad-hoc hypothesis of largescale hydrothermal circulation in stable ocean crust. To place the new model in context, we summarize the main characteristics and inherent limitations of currently accepted thermal models of the lithosphere. Next, the characteristics of thermal fields associated with upwelling of asthenospheric materials are outlined and its compatibility with the new model features examined. Following this, details of the new model fits to observational data on heat flow and bathymetry are presented, along with results of numerical simulations exploring the influence of model parameters. We point out, in addition, that empirical relations such as those proposed for GDH reference models are unnecessary. Implications of the new model results for understanding regional scale variations in global heat flow are discussed and the need to downsize the current estimates of global heat loss emphasized.

Current Models of the Lithosphere
Thermal models of the lithosphere, with wide acceptance in the current literature, may be classified as falling into essentially two generic groups: In the HSC model the basic assumption is that the temperature of the medium at origin time (t = 0) has a constant value T m for all depths. This constant T m then holds for all time at infinite depth. The lithosphere is considered as boundary layer of the mantle convection cells, arising from near surface conductive cooling. The lithosphere (in other words, the boundary layer) grows in thickness continuously as it moves away from the up-welling limb of the mantle convection system. Analytical expressions for temperature variation of the boundary Layer may be obtained as solution to the one-dimensional heat conduction equation [7].
The boundary layer approach has been successful in accounting for first order features in variation of oceanic heat flow with age (e.g. [38], [39], [44]). Nevertheless, this model cannot be considered as satisfactory for several reasons.
To begin with the HSC model is strictly valid for heat flux arising from cooling of a stagnant body and not one in which lateral movements occur in response to thermal convection. This is in direct contradiction with one of the essential ingredients of thermal convection, that of lateral movements. In addition, the model predicts infinite heat flow at the ridge axis (the well known problem of singularity in heat flow at time of origin) and heat flow values about five-fold higher than those observed in regions close to the ridge axis. The problem of high model heat flow for young ocean crust is a direct consequence of specific boundary and initial conditions imposed in the Half-Space Cooling model. For example, the assumption of constant temperature for the region beneath the boundary layer imply that asthenosphere resemble an isothermal "magma-ocean". The assumption itself is incompatible with the well-known characteristics of natural convection systems (e.g. [8], [11], [47]). The relatively low heat flow values predicted for older segments of the lithosphere is yet another characteristic feature of the Half-Space Cooling model. It is a consequence of the assumption (e.g. [45], [47]) that boundary layer growth (in other words, the steady advance of the solidification front) takes place exclusively due to heat loss from the upper surface (see discussion in the next section).
It was pointed out by McKenzie [30] that the difficulty with low model heat flow at large distances from the spreading centers can be overcome by assuming that the thickness of the lithosphere at large distances from the ridge axis approaches a constant value. This came to be known later as the Plate model and was adopted in several later studies ( [32], [39]). In regions close to the ridge axis the interval considered as "Plate" also includes the underlying asthenospheric wedge. The Plate model also assumes that basal temperature is constant beneath vast stretches of stable ocean lithosphere (in other words, asthenosphere is isothermal) and that magma injection occurs only at the central plane of the ridge axis. The Plate model approach has been successful in accounting for the first order features in variation of heat flow with age at large distances from the spreading centers. However, as in the case of the HSC model, it also predicts heat flow values much higher than the observed ones for ocean crust with ages less than 55 Ma.
The assumptions in the Plate model that basal temperature and thickness of the lithosphere are constant rely on the argument that lateral movements of surface layer take place over large nearly isothermal cores present in mantle convection systems. While these may be true of oceanic lithosphere away from spreading centers they can hardly be considered as representative of the thermal structure in regions close to the ridge axis, where non-isothermal conditions are likely to prevail at the base of the lithosphere. In particular, the assumption of constant basal temperature in zones overlying upwelling limbs of asthenosphere contradicts the vast body of observational evidences on temperature variations in intrusive magmatic and thermal metamorphic processes (e.g. [5], [28], [52], [61]).
In an apparent attempt to minimize problems of this type it has been proposed ( [39], [53]) that the relations derived from the HSC and Plate models may be combined in such a way that their characteristic constants are compatible with theoretical estimates of heat flow for oceanic crust with ages less than 55 Ma and with experimental heat flow data for ages greater than 55 Ma. In other words, the heat flow -age relations become hybrid in character, a result of the sequential use of solutions of the HSC and Plate models for separate age intervals. The selection of age ranges has been somewhat arbitrary, based on the best match to the data.
The hybrid models of Stein and Stein [53] has since then been accepted in the relevant literature as Global Depth Heat Flow (GDH) reference models.
Yet, significant discrepancies continue to exist between the hybrid model values and observational heat flow data, for oceanic lithosphere with ages less than 55 Ma. The current consensus is that such differences arise from perturbing effects of supposed regional scale hydrothermal circulation in the ocean crust, believed to be unaccounted for in conventional heat flow measurements in the oceanic crust ( [39], [41], [56], [60]). The origins of some of the problems with the hybrid versions can be traced back to the boundary conditions imposed in the Half-Space Cooling and Plate models. As pointed out recently by Hamza et al [16] GDH reference models imply discontinuities in the deep temperature fields of the lithosphere. In fact the GDH model requires an artificial change in heat flow at 20 My in order to fit the bathymetry data while another artificial change at 55My is necessary to fit the heat flow data [16]. In view of such limitations, the hybrid model approach can hardly be considered as a satisfactory alternative to understanding the thermal structure of the lithosphere.

Magma Accretion Model of the Oceanic Lithosphere
We consider now a new thermal model of oceanic lithosphere that can overcome some of the shortcomings of the HSC and Plate models discussed in the previous section. Following the premises of the previous models we also assume that lithosphere represents the boundary layer of mantle convection and that its temperatures are always at or below the melting temperature. In developing the new model it is assumed that the growth of this boundary layer, in regions away from the ridge axis, is determined not only by the cooling effects of surface heat loss but also by mass and energy exchange processes taking place at the bottom boundary of the lithosphere. In particular, we consider that the effects of basal  [37], [38], [49]). This is the classical case of boundary layer growth in stagnant fluids ( [8], [9], [11], [51]).
However, in cases where accretion is determined also by temperature variations where L 0 is the stable thickness of the plate at large distances from the ridge axis and η an appropriate scaling factor which may be considered as a measure of the change in the degree of basal accretion. Note that at L(x) = 0 at x = 0 (i.e. at the ridge axis) and L(x) = L 0 at x = ∞ (i.e. in stable ocean basins). Also, the second member on the RHS of equation (1), given by: may be considered as the thickness of the column of asthenospheric material between the base of the solid lithosphere and the level of the asthenosphere in stable ocean basins. According to equation (1b) the height of this asthenospheric column decreases from L m = L 0 at x = 0 to L m = 0 at x = ∞.
The growth of boundary layer is also affected by the presence of lateral temperature variations in the asthenosphere. The assumption of lateral temperature variations is compatible with the vast body of observational evidences on temperature fields in magmatic and thermal metamorphic processes (e.g. [5], [28], [52], [61]) and experimental data on thermal structures of convective plumes ( [35], [36], [59]). In discussing geophysical constraints on mantle temperatures Solomon [50] refers to mineral thermometric data for primary magmas of deep origin penetrating the oceanic lithosphere. The conclusion of Solomon [50] is that the asthenospheric temperatures near the spreading centers are in the range 1200 to 1250 0 C, nearly 200 degrees higher than the corresponding values beneath stable ocean basins. Lateral temperature variations are, in fact, a rule rather than an exception in many of the tectonothermal processes in the upper crust. In addition, there are no physically plausible reasons to believe that asthenospheric upwelling take place under isothermal conditions.
In the present case, we assume that the temperature variation in the asthenosphere, along a horizontal plane at the depth corresponding to the base of the stable lithosphere, is best represented by a relation of the type: where T a is the temperature of the asthenosphere in the upwelling regions, T b (x) its temperature at distance x, T m its temperature at large distances from the ridge axis and c 1 a scaling constant. Note that equation (2) describes the variation of temperature in the upwelling zone of the asthenosphere. The temperature T b (x) may also be considered as the basal temperature of the interface zone between the asthenosphere and the lithosphere. However, it should not be confused with the temperature at the base of the lithosphere, this latter one remains constant at the solidification temperature T S and is independent of the distance from the ridge axis.
As mentioned earlier, the main consequence of basal accretion and lateral temperature variations is an increase in the rate of "migration" of the solidification isotherm to larger depths relative to those encountered for isothermal fluids ( [9], [51]). As a result the width of the zone of partial melting is narrower in VBA model compared to those in HSC and Plate models. A schematic illustration of this fundamental difference between the VBA and Plate models is illustrated in Figure   (  and variable temperature asthenospheric upwelling. Note that the magma injection zone, whose geometry is determined by the depth to the solidification isotherm, is wider for constant temperature case relative to that for variable temperature.
In the following sections we consider the mathematical basis of the new VBA model and compare the model predictions against observational heat flow and bathymetry data for oceanic regions. In addition, we also compare VBA model values of heat flow and bathymetry with those derived from the half-space cooling and Plate models.

Theoretical Formulation
Consider first the problem of two dimensional heat transfers in a rectangular plate of thickness L moving with velocity v in the horizontal (x) direction. In the model discussed by McKenzie [30] both the thickness of the plate and its basal temperature are assumed to be constants. In the VBA model of the present work these parameters are considered as space dependant variables. As discussed in the previous section, the form of variation of lithosphere thickness (L) is assumed to be determined by equation (1) while the systematic decrease in the temperature beneath its base (ie: at the top of the asthenosphere) is assumed to be determined by equation (2). Impositions of these conditions however make the heat transfer problem under consideration non-linear, for which there are no easy analytical solutions. One of the convenient ways of overcoming problems of this type is to make use of the standard method of piece-wise approximation. In this approach, the medium is assumed to be composed of a system of discrete elements, the spatial domains of which are chosen to be sufficiently small that the effects of changes in parameter values (in the present case, the plate thickness and the fusion temperature) may be considered as negligible within each individual element. The relevant differential equation for any particular element of this system is: where ρ is density, C P the specific heat, T the temperature, t the time, λ the thermal conductivity, g the rate of heat generation and x and z the horizontal and vertical coordinates respectively. The origin of the coordinate system is fixed at lower left corner of the rectangular element under consideration. The subscript j refers to the discretization index, which assumes values 0, 1, 2, 3, …n, n being the number of elements. The boundary conditions are: where ψ i are the eigen functions, µ i the eigen values and N the norm of the solution. Note that equations (5), (6) and (7) are derived from the corresponding equations (A28), (A29) and (A30) of the dimensionless variables in the Appendix.
The terms f i and a 1 are given by the relations (see Appendix): In equation (11) Pe is the Peclet number, given by the relation: The solutions (5), (6) and (7) 6) and (7). It must however be pointed out that the solutions obtained for the system of discrete elements are coupled in the sense that the solution for element (j+1) is derived from the solution for the previous element (j). The equation relating the temperature profile on the right hand side of the i th set of elements constituting the lithospheric block of thickness L (at lateral position X) to the left hand side of the (i+1) th set of elements constituting the adjacent block of thickness L+∆L (at lateral position X+dX) is given by: Transfer of boundary temperature profile from a block of thickness L to the next one with larger thickness L+ΔL leads to a reduction in the value of the temperature gradient. It is a consequence of three parallel processes: a) The lithospheric block of larger thickness is positioned over a region of asthenosphere with relatively lower temperatures; b) There is a reduction in the rate of basal accretion of magma; and c) Loss of heat by the lithospheric block in the vertical direction towards the surface.
Note that the thermal effects of the first and second processes were not taken into account in previous models of the lithosphere (HSC and Plate models), a consequence of the assumption of isothermal asthenosphere in these models.

VBA Model Fit to Observational Heat Flow Data
We now make a comparative analysis of the VBA model predictions with results of heat flow measurements in oceanic regions. Following earlier studies (e.g. [53]) we also consider data for oceanic lithosphere with ages less than 160 Ma. The parameter values used in model calculations are given in Table (1), which are essentially identical to those used in earlier studies, allowing thereby direct comparison. Our value of thermal conductivity (λ) is slightly larger, a compromise consistent with results of modern measurements ( [40], [19]). The average value of thermal diffusivity (κ ) used (0.8mm 2 /s in Table -1) is below the average from 298 to 1300K of modern data for a 60 olivine -30 clinopyroxene -10% garnet composition ( [4], [18]). Our use of lower values than modern averages is consistent with the lower lithosphere exerting a stronger role in controlling heat flow, as it is both more insulating than the upper layers and because heat from the magma must first traverse these deepest lithosphere. The variation of VBA model heat flow with age, determined on the basis of equation (7) and the parameter values in Table (1), is illustrated in Figure (

Comparison between Plate and VBA Models
where T a is the temperature of the lateral boundary, L the plate thickness and κ the thermal diffusivity of the medium.  (3) in a form similar to that proposed recently by Cardoso and Hamza [6] for the variable basal heat input problem: There are several important differences between the solutions (11) and (12).
To begin with we note that the source strength of the transient component (determined by the pre exponential terms in the summation series of equation (11) Note that, with the exception of the argument of the sinusoid, equation (13) is similar to the relation derived by Royden and Keen [43] for a lithosphere undergoing stretching due to intrusions. The relation for heat flux, derived from (13) is: where λ is the thermal conductivity. For n > 1, the higher order terms in the summation series on the right hand side of equation (14a) are practically negligible.
It may therefore be simplified as: An interesting aspect of equations (13) and (14)  Note that the derivation of equations (13) and (14) assumes constant λ, and hence, κ represents the average thermal diffusivity over the lithosphere. Although Equations (13) and (14) pertain to the limit z = L, suggesting that room temperature values for λ may be the most appropriate, it can also be argued that average values are consistent with a constant κ-λ derivation. Average values were used in all previous Plate and Boundary Layer models.

VBA Model Fit to Bathymetry Data
Fits to data for ocean floor bathymetry variations rather than that for surface heat flow is often considered as a relatively more rigorous test of thermal models of the lithosphere. The relation for bathymetry in VBA model has been developed following the isostatic compensation scheme discussed in earlier studies (e.g. [30], [48], [39]). Consider, for example, the mass balance relations for two columns: one situated at the ridge and the other one away from the ridge. For a constant transverse section of the lithosphere the relation for isostatic balance between these columns is: Where ρ w is the density of sea water, d r the elevation of the ocean ridge above the final level to which the lithosphere subsides and ρ ast the density of the asthenosphere. The terms e(x) and ε(x) represent, respectively, the elevation of the sea floor and the elevation of the base of the lithosphere, dh the thickness of infinitesimal volume element where the temperature change is taking place, ρ 0 the reference density of the base of the lithosphere, β the volumetric expansion coefficient and ∆T the temperature difference between the base of the lithosphere and the volume element. Equation (15) has been derived assuming that The integration in equation (15) is to be carried out over the entire thickness of the lithosphere whose basal temperature is T i . Developing the integral on the RHS of equation (15), after substitution for the temperature from equation (13), we have: On the other hand, the integral I 2 in equation (18) where is the constant of proportionality which has units of m -1/2 . Substituting (20) in (18), changing the limits of integration and solving for the integral:  (23) In comparing the bathymetry results of VBA model with the observational datasets we make use of the same data sets ( [57], [26], [48]) and procedure as that employed in the previous study of Stein and Stein [53], allowing thereby direct comparison. The bathymetry data employed refer to values averaged over 2-Ma bins and these have been used also in our comparative analysis. A comparison of VBA model fit to this bathymetry dataset, where the calculations make use of the parameter values given in Table ( 1), is presented in Figure (

Discussion and Conclusions
In contrast to HSC and Plate modes, which are closely related regarding the assumptions made and in their implementation, the newly proposed VBA model of oceanic lithosphere assumes that both the thickness and the temperature of the magma rich basal segment vary with distance from the ridge axis. Estimates of  The discrepancy of previous models of global heat flow from the experimental values has been hypothesized to originate from regional scale hydrothermal circulation in oceanic crust. As mentioned earlier, the validity of the hypothesis of regional scale hydrothermal circulation in oceanic crust is questionable, in view of available information on the thermal and hydrological characteristics of the ocean crust ( [20], [16]).
A problem of related interest is the Global Heat Loss. The VBA Model leads to estimates of regional heat flow lower than those derived from previous models, in agreement with recent results of higher degree harmonic representation of global heat flow (Hamza et al, [15] 9 -Estimates of root-mean-square misfit between the VBA Model values and experimental heat flow data are relatively much better than those found for the previous models. Given the uncertainty in marine heat flow measurements and the quality of the fit, there appears to be no need to invoke the hypothesis of regional scale hydrothermal circulation in oceanic crust; 10 -The VBA Model of the present work leads to estimates of regional heat flow that are significantly lower than those derived from previous thermal models of the lithosphere. The new estimates are in reasonable agreement with the results of higher degree harmonic representations of global heat flow (Hamza et al., [15]); and 11 -The current estimates of global heat loss need to be downsized by at least 25%, in support of recent assessments ( [20], [15], [16]) and the view that the Earth is in quasi-steady thermal state.

Acknowledgments
We thank Prof. In addressing the problem described by equations (3) and (4)  This allows us to rewrite the differential equation (3) for the system of n discretized elements as: where j is the discretization index, X 0 > 0 and Pe is the Peclet number given by: The interval [X j , X j+1 ] can be chosen to be sufficiently small that the Peclet number may be considered constant within any specific interval. We assume steady state conditions and consider that heat production term G is negligible. In this case the differential equation and the boundary conditions become: The purpose of condition (A4b) is to avoid the well known problem of singularity at the position X = 0.
We assume that θ (X, Z) may be expressed as: The problem in K (Z) is: Hence the solution of problem in K (Z) is: The problem in u (X, Z) is: The condition (A7b) is necessary for the solution u(X, Z) to be compatible with the condition (A4e).
We now admit that the solution of the problem in u (X, Z) may be expressed as: The integral on the RHS of (11) is zero for i ≠ j (eigen functions are orthogonal, see Özisik, 1980). For i = j, the integral leads to the norm N i , associated with the problem: Consequently, the unknown C i (X) is given by: It is obvious that the transform and its inverse in equation (A14) are: This procedure transforms the partial differential equation into a system of ordinary differential equations. We adopt the following strategy: