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

Abstract

A magma accretion model of oceanic lithosphere is proposed and its implications for understanding its thermal field examined. The new model (designated Variable Basal Accretion—VBA) assumes existence of lateral variations in magma accretion rates and temperatures at the boundary zone between the lithosphere and the asthenosphere. However, unlike the previous thermal models of the lithosphere, the ratio of advection to conduction heat transfer is considered a space dependent variable. The results of VBA model simulations reveal that the thickness of the young lithosphere increases with distance from the ridge axis, at rates faster than those predicted by Half-Space Cooling models. Another noteworthy feature of the new model is its ability to account for the main features in the thermal behavior of oceanic lithosphere. The improved fits to bathymetry have been achieved for the entire age range and without the need to invoke the ad-hoc hypothesis of large-scale hydrothermal circulation. Also, use of VBA model does not lead to artificial discontinuities in the temperature field of the lithosphere, as is the case with GDH (Global Depth Heat Flow) reference models. The results suggest that estimates of global heat loss need to be downsized by at least 25%.

1. 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 [1] and Plate [2] 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 [3]. 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., [3–7]).

There are however dissenting views on the subject matter of hydrothermal circulation on regional scales [8].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., [9–13]). No direct experimental evidence has so far been presented that point to the existence of 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., [5, 14–16]).

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 large-scale 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.

2. 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:(i)Half-Space Cooling (HSC) Models,(ii)Constant Thickness Plate Models.

In the HSC model the basic assumption is that the temperature of the medium at origin time (t = 0) has a constant value 𝑇_𝑚 for all depths. This constant 𝑇_𝑚 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 [17]. The boundary layer approach has been successful in accounting for first order features in variation of oceanic heat flow with age (e.g., [18–20]). 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 fivefold 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 implies that asthenosphere resembles an isothermal “magma-ocean”. The assumption itself is incompatible with the well-known characteristics of natural convection systems (e.g., [21–23]). 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., [23, 24]) 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 [2] 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 [19, 25]. 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 nonisothermal 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., [26–29]). The available experimental data on thermal structures of convective plumes also point to the existence of lateral variations in convective systems (e.g., [30–32]). Other problems associated with the HSC and Plate models have been discussed in recent works (e.g., [8, 33–35]).

In an apparent attempt to minimize problems of this type it has been proposed [19, 36] 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 C. A. Stein and S. Stein [36] have 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 [3, 5, 7, 19]. 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. [37] 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 55 My is necessary to fit the heat flow data [37]. 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.

3. 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 magma accretion and lateral temperature variations of the asthenosphere play important roles in the formation of the lithosphere. The new model is designated hereafter as the Variable Basal Accretion model, abbreviated VBA.

The basal accretion may take place as a result of pressure and temperature variations in the ascending magma column, compositional changes occurring during up-flow and differential rates of migration of volatile components. It is well known in fluid dynamics studies [38, 39] that the material and thermal exchange processes occurring in the transition zone immediately below the depth of solidification isotherm have variable degrees of viscous coupling with those occurring in regions of fully developed flows. This line of reasoning leads to the deduction that the interaction between the asthenospheric flow and the lower boundary of the lithosphere becomes less significant as the plate moves away from the ridge axis. Consequently, the amount of heat advected into the basal parts of the lithosphere decrease with distance from the spreading center.

In the present context of developing a new thermal model of the lithosphere the main interest is in examining the effect of basal accretion on the thickness of the lithosphere and on the surface heat flux. If accretion can be considered as a consequence exclusively of conductive heat loss from the upper surface of the lithosphere, an approximate description of the boundary layer growth and the ensuing heat flow variation at the surface can be provided on the basis of the well-known inverse square root relation of the time elapsed [17, 18, 40, 41]. This is the classical case of boundary layer growth in stagnant fluids [21, 22, 42, 43]. However, in cases where accretion is determined also by temperature variations (and ensuing chemical or gravitational differentiation processes) in the stagnant layer (between the solid and liquid parts) the rate of boundary layer growth differs from that in the classical case. Quantitative assessments of such changes in accretion are difficult in view of the limited knowledge of the thermal state of matter and of the chemical and gravitational processes operating at the lithosphere-asthenosphere interface. For the purposes of the present work we assume that the variation in thickness of the lithosphere (𝐿) with distance x may be represented by a relation of the type:𝐿(𝑥)=𝐿₀1[1−(1+𝜂√𝑥)], (1a) where 𝐿₀ 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 𝐿(𝑥)=0 at 𝑥=0 (i.e., at the ridge axis) and 𝐿(𝑥)=𝐿₀ at 𝑥=∞ (i.e., in stable ocean basins). Also, the second member on the RHS of (1a) and (1b), given by𝐿_𝑚𝐿(𝑥)=⁰(1+𝜂√𝑥), (1b) 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 (1b) the height of this asthenospheric column decreases from 𝐿_𝑚=𝐿₀ at 𝑥=0 to 𝐿_𝑚=0 at 𝑥=∞.

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., [26–29]) and experimental data on thermal structures of convective plumes [30–32]. In discussing geophysical constraints on mantle temperatures Solomon [44] refers to mineral thermometric data for primary magmas of deep origin penetrating the oceanic lithosphere. The conclusion of Solomon [44] is that the asthenospheric temperatures near the spreading centers are in the range 1200 to 1250^∘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 takes 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:(𝑇_𝑎−𝑇_𝑏(𝑥))(𝑇_𝑎−𝑇_𝑚)=1−exp(−𝐶₁𝑥), (2) where 𝑇_𝑎 is the temperature of the asthenosphere in the upwelling regions, 𝑇_𝑏(𝑥) its temperature at distance 𝑥, 𝑇_𝑚 its temperature at large distances from the ridge axis, and 𝑐₁ a scaling constant. Note that (2) describes the variation of temperature in the upwelling zone of the asthenosphere. The temperature 𝑇_𝑏(𝑥) 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 𝑇_𝑆 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 [42, 43]. As a result the width of the zone of partial melting is narrower in VBA model compared to that in HSC and Plate models. A schematic illustration of this fundamental difference between the VBA and Plate models is illustrated in Figure 1. The VBA model prediction for a narrower width of the magma injection zone seems to be supported experimental heat flow data for oceanic regions. For example, the width of heat flow anomalies in mid-ocean ridge zones usually has dimensions much smaller than that predicted by the Half-Space Cooling (HSC) and Plate models.

Figure 1 

Schematic illustration of solidification isotherms (𝑇𝑆) for cases of constant 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.

3.1. 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 (𝑥) direction. In the model discussed by McKenzie [2] 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 (𝐿) is assumed to be determined by (1a) and (1b) while the systematic decrease in the temperature beneath its base (i.e., at the top of the asthenosphere) is assumed to be determined by (2). Impositions of these conditions however make the heat transfer problem under consideration nonlinear, 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 piecewise 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𝜌𝐶_𝑝[𝜕𝑇𝜕𝑡+𝑣𝜕𝑇𝜕𝜕𝑥]=𝜆[²𝑇𝜕𝑥²+𝜕²𝑇𝜕𝑧²𝑥]+𝑔_𝑗<𝑥<𝑥_𝑗+1,𝑡>0,0<𝑧<𝐿, (3) where 𝜌 is density, 𝐶_𝑃 the specific heat, 𝑇 the temperature, 𝑡 the time, 𝜆 the thermal conductivity, 𝑔 the rate of heat generation, and 𝑥 and 𝑧 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 𝑗 refers to the discretization index, which assumes values 0,1,2,3,…𝑛,  𝑛 being the number of elements. The boundary conditions are𝑇(𝑥,𝑧<0)=𝑇_𝑏(𝑥),𝑇(𝑥,𝑧=0)=𝑇_𝑆,𝑇(𝑥,𝑧=𝐿)=0,𝑇(𝑧,𝑥=0)=𝑇_𝑎, (4a)(4b)(4c)(4d) In (11) Pe is the Peclet number, given by the relation:Pe(𝑥)=𝜌𝐶_𝑝𝑣𝜆𝐿(𝑥)𝑥_𝑗<𝑥<𝑥_𝑗+1. (9b)

The solutions (5), (6), and (7) are similar to the respective relations derived by McKenzie [2] for the constant temperature Plate model. An important difference is the presence of the coefficient a₁ in the exponential terms. The value of this coefficient (see (9a) and (9b)) depends on the Peclet number, which in turn depends on the plate thickness L (see (1a), (1b), and (10)). In implementing the discretization scheme the size of the elements is made sufficiently small, that the Peclet number may be considered as constant within each individual element, allowing thereby use of the solutions (5), (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 (𝑗+1) is derived from the solution for the previous element (j). The equation relating the temperature profile on the right-hand side of the 𝑖th set of elements constituting the lithospheric block of thickness L (at lateral position X) to the left-hand side of the (𝑖+1)th set of elements constituting the adjacent block of thickness 𝐿+Δ𝐿 (at lateral position 𝑋+𝑑𝑋) is given by𝑇(𝑋_𝑗,𝑍_𝑗)_𝑖=𝑇(𝑋_𝑗+1,𝑍_𝑗)_𝑖+1. (10)

Transfer of boundary temperature profile from a block of thickness L to the next one with larger thickness 𝐿+Δ𝐿 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.(c)There is 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.

At this point a brief remark on the scaling constant 𝜂 is in order. Note that the numerical value of this parameter is inversely proportional to the rate of change in basal accretion. Thus, a null value of 𝜂 means that the thickness of the lithosphere is independent of the distance from the ridge axis which implies that the Peclet number is constant. Nonzero positive values of 𝜂 refer to the more general case considered in the present work, where the rates of magma accretion fall off systematically with distance from the ridge axis. Thus, in the VBA model, small values of 𝜂 mean that the total amount of accreted material is large and consequently the approach to stable thermal conditions of the lithosphere is relatively slow. On the other hand, large values of 𝜂 mean higher rates of initial accretion, but relatively rapid approach to the final stable conditions. In such cases, decrease of heat flow with age takes place at rates higher than those predicted by the HSC and Plate models. The decrease in thickness of the column of asthenospheric material with distance from the ridge axis (which is inversely proportional to the increase in thickness of the lithosphere), for different values of the basal magma accretion factor 𝜂, is illustrated in Figure 2.

Figure 2 

Decrease in thickness of the column of asthenospheric material as a function of the distance from the ridge axis,. The numbers on the curves are values of the parameter (𝜂) in (1a) and (1b). See text for details.

Results of numerical simulations indicate that the VBA model values of heat flow for the case 𝜂=0 are practically identical to the values derived from the Plate model of McKenzie [2]. This is not altogether surprising since the Plate model assumes that magma accretion occurs only at the ridge axis. It is clear that the Plate model envisaged by McKenzie [2] is a particular end-member case of the more general class of VBA models.

The intervals chosen in discretization of VBA model simulations are in the range of 10 m to 1 km. For gradual changes in the thickness of the lithosphere the computational accuracy of the results obtained in this piecewise approximation is not overly sensitive to the size of the interval chosen for discretization. On the other hand, the approach has the advantage that the effects of lateral temperature variations arising from compositional changes, which determine the variability in the lower boundary condition (see (2) and (4a)), can be taken into account. A number of numerical simulations were carried out as part of sensitivity tests of model response to parameter values given in Table 1. An example of the results is given in Figure 3, for the case of uniform physical properties of the lithosphere. It reveals that the temperature field is characterized by smooth and continuous changes and devoid of the presence of discontinuities and inversions, provided that the changes in temperatures at the top of the asthenosphere are within reasonable limits.

Table 1 

Values of parameters used in numerical simulations of the VBA model of the oceanic lithosphere.

Figure 3 

Distribution of isotherms in the oceanic lithosphere, derived from the VBA model of the present work. The numbers on the curves are temperatures in degrees centigrade.

3.2. 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., [36]) 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 [47, 48]. The average value of thermal diffusivity (𝜅) used (0.8 mm²/s in Table 1) is below the average from 298 to 1300 K of modern data for a 60 olivine −30 clinopyroxene −10% garnet composition [49, 50]. 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 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 (7) and the parameter values in Table 1, is illustrated in Figure 4. Also included in this figure are the mean heat flow values reported by C. A. Stein and S. Stein [36] for ocean crust with ages of up to 100 Ma.

Figure 4 

Comparison of VBA model values with experimental heat flow data for the oceanic lithosphere. The model values for basal accretion rates 𝜂 between 0.5 and 2 are indicated by the set of blue, black, and green curves. The red colored curve refers to the case 𝜂=0. It represents the lower end-member case of the family of VBA models and is identical to that derived from the Plate model [2].

As mentioned earlier, the VBA model leads to a family of solutions depending on the value of the factor 𝜂 that determines the basal accretion rate. The model curves for values of 𝜂 in the range of 0.5 to 2 bracket most of the experimental heat flow data for ocean crust with ages less than 55 Ma. For ages greater than 55 Ma the model curves, independent of the value of 𝜂, tend towards an asymptotic limit. It is clear that VBA Model is capable of providing vastly improved fit to marine heat flow data, relative to that which can be achieved within the framework of HSC and Plate models.

We now examine the dependence of the VBA model predictions on the plausible variations in the values of the main parameters. For this purpose, a number of numerical simulations were carried out for a range of values of the thickness of the plate (𝐿), temperature of the stable segment of the asthenosphere (𝑇_𝑖), and thermal diffusivity (𝜅). Results of some of the numerical simulations are illustrated in the set of panels in Figure 5.

(a)

(b)

(c)

Figure 5 

Results of numerical simulations illustrating the dependence of VBA model response to changes in the values assumed for plate thickness (a), basal temperatures (b) and thermal diffusivity (c). The dark squares are the mean oceanic heat flow values.

Note that changes in thickness of the lithosphere (Figure 5(a)) and its basal temperature (Figure 5(b)) have only minor influence on the VBA model results. However, the value of thermal diffusivity is found to have a marked influence on the model predictions, as can be seen in the results illustrated in Figure 5(c). The diffusivity values compatible with experimental heat flow data lie in the range of 11 to 43×10⁻⁶ km²/yr, which is in reasonable agreement with the values of thermal diffusivity of representative mantle minerals [34, 47]. The highest and lowest values of 𝜅 used in our test are extremes observed at 25^∘C and 1027^∘C. This exercise shows that the temperature dependence of 𝜅, had it been taken into account in a much more complex model, would provide model heat flow values that are consistent with experiments. Reproduction of observational data by the VBA model is not predicated on the specific parameters used in our calculations: values of 𝐿, 𝜅, 𝜆, and 𝑇_𝑖 that reasonably represent the lithosphere all lead to excellent agreement of the model with the observables.

3.3. Comparison between Plate and VBA Models

The classical solutions for transient temperature distributions in a plate with constant boundary temperatures have been discussed extensively in the literature [17, 40]. The Plate model of McKenzie [2] assumes that the transient temperature at the left lateral boundary is equal to the difference between the bottom boundary temperature and the steady state geotherm. The transient temperature in this case is given by the relation (McKenzie, [2]):𝑇_𝑡(𝑧,𝑡)=2𝑇_𝑎𝜋^∞𝑛=0(−1)^𝑛𝑛sin(𝑛𝜋𝑧𝐿)exp(−𝑘𝑛²𝜋²𝑡𝐿²), (11) where 𝑇_𝑎 is the temperature of the lateral boundary, 𝐿 the plate thickness, and κ the thermal diffusivity of the medium.

A direct comparison between the transient component of VBA model solution (second term on the RHS of (5)) and the solution for Plate Model (11) is not straightforward. It is more convenient in this context to adopt the Fourier series solution approach of McKenzie [2] and recast the solution of (3) in a form similar to that proposed recently by Cardoso and Hamza [51] for the variable basal heat input problem:𝑇_𝑡(𝑧,𝑡)=2𝑇_𝑎𝜋^∞𝑛=11𝑛sin(𝑛𝜋(𝐿−𝑧)𝐿𝜂√𝑥)𝑒^{−(𝜅𝑛2𝜋2𝑡)/𝐿2}, (12) where 𝜂 is a scaling constant (with dimensions of m^−1/2), and 𝑥 the distance from ridge axis. It is fairly straightforward to show that solution (12) satisfies the initial and boundary conditions of the relevant heat conduction equation of the Plate model. Results of numerical simulations show that the transient temperatures calculated using the solution (12) are nearly identical to the corresponding values obtained from the second term on the RHS of (5), provided that an appropriate value is chosen for the scaling constant 𝜂. In the present case, the value chosen for 𝜂 is 0.4. The advantage of using (12) is that a direct comparison is now possible against the solution obtained for the constant temperature Plate model of McKenzie [2].

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 preexponential terms in the summation series of (11)) in the Plate model is independent of the distance from ridge axis. Consequently, the magnitude of the transient perturbation in the Plate model remains relatively high for large distances from the ridge axis. On the other hand, the source strength of transient perturbation in VBA model (determined by the pre exponential terms in the summation series of (12)) is dependent on the distance from ridge axis, and in addition the argument of the sine function is scaled down by the factor (𝐿−𝑧)/(𝜂√𝑥). Thus, the magnitude of transient perturbation is relatively small for all times in the VBA model compared to that of the Plate model.

Another important difference between the solutions (11) and (12) may be illustrated by considering the vertical distribution of transient temperature components. As can easily be verified from (11) the transient component in the Plate model has a maximum in the upper parts of the lithosphere but is zero at its top and bottom boundaries, a consequence of the imposed constant temperature boundary conditions. Initially, the vertical distribution is asymmetric with respect to the central plane of the lithosphere, as illustrated in Figure 6. The asymmetry is a consequence of the fact that the transient component at the left lateral boundary is set as the difference between the bottom boundary temperature and the steady-state geotherm. Nevertheless, the overall shape of the transient response is quite different from that expected for perturbations arising from high-temperature intrusions penetrating the base of a relatively cold plate [17, 40, 52, 53].

Figure 6 

Vertical distribution of transient temperatures in the Plate model (McKenzie, [2]), for the parameter values listed in Table 1. The numbers on the curves are temperatures in degrees centigrade.

The fact that the transient component is absent for all times at the base of the lithosphere is an indication that the Plate model is inadequate in providing a satisfactory description of heat transfer processes at the lithosphere—asthenosphere boundary, mainly in regions close to the ridge axis. On the other hand, the vertical distribution of the transient component in the VBA model, illustrated in Figure 7, reveals that it is relatively large at the base of the plate and zero at the top boundary. Its overall shape is thus similar to that expected for perturbations generated by a relatively high-temperature intrusive penetrating the lower boundary of the lithosphere.

Figure 7 

Vertical distribution of transient temperatures in the VBA model of the present work, for the parameter values listed in Table 1. The numbers on the curves are temperatures in degrees centigrade.

Including the steady-state component in (12) leads to the complete solution for temperature:𝑇(𝑧,𝑡)=𝑇_𝑎𝑧(1−𝐿)+2𝑇_𝑎𝜋×^∞_𝑛=11𝑛sin(𝑛𝜋(𝐿−𝑧)𝐿𝜂√𝑥)𝑒^{−(𝜅𝑛2𝜋2𝑡)/𝐿2}. (13) Note that, with the exception of the argument of the sinusoid, (13) is similar to the relation derived by Royden and Keen [54] for a lithosphere undergoing stretching due to intrusions. The relation for heat flux, derived from (13) is𝑇𝑞(𝑡)=𝜆^𝑎𝐿2(1+𝜂√𝑥^∞𝑛=1𝑒^{−(𝜅𝑛2𝜋2𝑡)/𝐿2}), (14a) where 𝜆 is the thermal conductivity. For 𝑛>1, the higher-order terms in the summation series on the right-hand side of (14a) are practically negligible. It may therefore be simplified as𝑇𝑞(𝑡)=𝜆^𝑎𝐿2(1+𝜂√𝑥𝑒^{−(𝜅𝜋2𝑡)/𝐿2}). (14b)

An interesting aspect of (13), (14a), and (14b) is that these exhibit features similar (though not identical) in character to those of the HSC and Plate models, depending on the time scale chosen. For example, for short times (i.e., for ages less than 55 Ma), the transient parts of the solutions in (13), (14a), and (14b) are determined mainly by the multiplication factors of the summation, the exponential terms themselves being much less significant. In this case, the temperature and heat flow variations are inversely proportional to the square root of the distance (equivalently, the age of the oceanic crust). Hence the behavior of the thermal field is similar to that predicted by the HSC model. On the other hand, for large times (i.e.,: for ages greater than 55 Ma), the transient part of the solution is determined mainly by the exponential terms. In this case the temperature and heat flow variations are similar to those predicted by the Plate model.

Note that the derivation of (13), (14a), and (14b) assumes constant 𝜆 and, hence, 𝜅 represents the average thermal diffusivity over the lithosphere. Although (13), (14a), and (14b) pertain to the limit 𝑧=𝐿, 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.

4. 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., [2, 19, 55]). 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(𝑒(𝑥)−𝑑_𝑟)(𝜌₀−𝜌_𝑤)=𝜌₀𝛽^𝑒(𝑥)_𝜀(𝑥)Δ𝑇𝑑ℎ, (15) where 𝜌_𝑤 is the density of sea water, 𝑑_𝑟 the elevation of the ocean ridge above the final level to which the lithosphere subsides, and 𝜌_ast the density of the asthenosphere. The terms 𝑒(𝑥) and 𝜀(𝑥) represent, respectively, the elevation of the sea floor and the elevation of the base of the lithosphere, 𝑑ℎ the thickness of infinitesimal volume element where the temperature change is taking place, 𝜌₀ the reference density of the base of the lithosphere, 𝛽 the volumetric expansion coefficient, and Δ𝑇 the temperature difference between the base of the lithosphere and the volume element. Equation (15) has been derived assuming that 𝜌_ast≈𝜌₀.

The integration in (15) is to be carried out over the entire thickness of the lithosphere whose basal temperature is 𝑇_𝑖. Developing the integral on the RHS of (15), after substitution for the temperature from (13), we have^𝑒(𝑥)_𝜀(𝑥)Δ𝑇𝑑ℎ=𝐼₁−𝐼₂, (16) where𝐼₁=^𝑒(𝑥)_𝜀(𝑥)(𝑇_𝑎−(𝑇_𝑎𝑧(1−𝐿𝐼)))𝑑ℎ,₂=^𝑒(𝑥)_𝜀(𝑥)(2𝑇_𝑎𝜋^∞𝑛=01𝑛sin(𝑛𝜋(𝐿−𝑧)𝐿𝜂√𝑥)𝑒^{−(𝜅𝑛2𝜋2𝑡)/𝑎2})𝑑ℎ. (17)(18) The integral 𝐼₁ in (17) represents contraction of the lithosphere in zones without significant magma injection. The limits of the integration are 𝜀(𝑥)=0 and 𝑒(𝑥)=𝐿. The result is straightforward𝐼₁=𝑇_𝑎𝐿2. (19) On the other hand, the integral 𝐼₂ in (18) represents the transient thermal effect of magma injection and its solution can be obtained by introducing a change of variable that takes into consideration lateral increase in the thickness (h) of the solid lithosphere. Since for any specific depth z there are different values of h, the relation between z and h may be expressed asℎ=𝐿−(𝐿−𝑧)𝜎√𝑥, (20) 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:𝐼₂=−2𝐿𝑇_𝑎𝜎𝜂𝜋^2∞𝑛=11𝑛²cos(𝑛𝜋𝑧𝐿)₀^{𝐿+𝑒(𝑥)}𝑒^{−(𝜅𝑛2𝜋2𝑡)/𝐿2}, (21a) we designate the product (ση) as the “bathymetry constant” 𝛿, and since 𝐿+𝑒(𝑥)≅𝐿, we have 𝐼₂=4𝐿𝑇_𝑎𝛿𝜋^2∞𝑛=01(2𝑛+1)²𝑒^{−(𝜅(2𝑛+1)2𝜋2𝑡)/𝐿2}. (21b) Substituting (19) and (21b) in (15) and after some obvious simplifications we have(𝑒(𝑥)−𝑑_𝑟)=𝜌₀𝛽𝑇_𝑎𝐿2(𝜌₀−𝜌_𝑤)(1−8𝛿𝜋^2∞𝑛=01(2𝑛+1)²𝑒^{−(𝜅(2𝑛+1)2𝜋2𝑡)/𝐿2}) (22) Designating 𝑑_𝑠=(𝑇_𝑎𝐿𝜌₀𝛽)/2(𝜌₀−𝜌_𝑤) we have the final solution:𝑒(𝑥)=(𝑑_𝑟+𝑑_𝑠)−8𝑑_𝑠𝛿𝜋²(^∞_𝑛=01(2𝑛+1)²𝑒^{−(𝜅(2𝑛+1)2𝜋2𝑡)/𝐿2}). (23)

In comparing the bathymetry results of VBA model with the observational datasets we make use of the same data sets [55–57] and procedure as that employed in the previous study of C. A. Stein and S. Stein [36], 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 8. Note that VBA model curve (in blue color) provides a remarkably good fit to bathymetry data (indicated by the asterisk marks) for the entire age range of the oceanic lithosphere. Also shown in this figure are the GDH model curves for bathymetry. The GDH model requires two separate curves (red and green curves), for arbitrarily selected age intervals, in achieving an equivalent fit.

Figure 8 

Comparison of the fits of VBA and GDH reference models to bathymetry data. Note that VBA model curve (in blue color) of the present work provides a remarkably good fit for the entire age range of the oceanic lithosphere. The GDH model [36] requires two separate curves (red and green curves), for arbitrarily selected age intervals, in achieving an equivalent fit.

Apart from the above mentioned restriction, both VBA and GDH models provide equally good accounts of ocean floor bathymetry. The vertical temperature field of the lithosphere, derived from the VBA model, is similar to the example illustrated in Figure 3. On the other hand, unlike those derived from GDH reference models there are no discontinuities in the temperature field of the lithosphere (Hamza et al. [37]).

At this point it is convenient to consider the sensitivity of VBA model response to the values of the parameters listed in Table 1. For young ocean crust (with ages less than 55 Ma) the main parameter that controls bathymetry is 𝛿, the best fit value of which is 0.6. The dashed and dotted curves in Figure 9 are model curves for 𝛿 values of 0.7, and 0.5, respectively. These model curves bracket the observational bathymetry data for ocean crust with ages less than 55 Ma.

Figure 9 

Results of numerical simulations illustrating the response of VBA model to change in the value assumed for the bathymetry constant 𝛿. The square symbols indicate mean oceanic bathymetry values.

For old ocean crust (with ages >55 Ma) the main parameter that controls bathymetry is the basal temperature, the best fit value of which is 1300 K. The dashed and dotted lines in Figure 10 are model curves for basal temperatures of 1400 K and 1200 K, respectively. These model curves bracket most of the observational bathymetry data for ocean crust with ages greater than 55 Ma.

Figure 10 

Results of numerical simulations illustrating the response of VBA model results to change in the value assumed for the basal temperatures (𝑇𝑖). The square symbols indicate mean oceanic bathymetry values.

5. 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 regional heat flux in the VBA model are lower than those obtained in previous thermal models of the lithosphere, including the recently proposed Plate model with variable thermal conductivity (McKenzie et al. [58]). Model heat flow values calculated on the basis of (9a) and (9b) may be used along with digital isochron data of Muller et al., [59] in mapping heat flow in oceanic regions, following a procedure similar to that suggested recently by Wei and Sandwell, [60]. In addition, theoretical values derived in this manner may be appended with experimental data for the continental regions in deriving global heat flow maps, an example of which is presented in Figure 11. Note that the global map of Figure 11 displays regional features in heat flow similar to those reported in recent 36 degree harmonic representation of IHFC data set (Hamza et al., [61]). Thus, while ridge areas have heat flow in excess of 80 mW/m², the remaining parts of oceanic regions and continental areas are characterized by heat flow less than 70 mW/m². Also, the global mean heat flow is 61 mW/m² while the maximum binned value is no more than 150 mW/m².

Figure 11 

Global heat flow map derived from mixed data sets. For oceanic regions with ages less than 120 Ma heat flow values are calculated using (14a) and (14b) with age values derived from a global grid at 0.1∘ spacing (Muller et al. [59]). For continental regions and remaining areas heat flow values are derived from IHFC data set, following the procedure outlined for higher-degree harmonic representation (Hamza et al. [61]).

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 [8, 37].

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., [61]). In the previous study of Pollack et al., [3] the global heat loss was estimated at about 44 TW. However, this value is based on the use of mixed data sets, which include both experimental data as well as theoretical values; those latter ones are based on the hybrid model of C. A. Stein and S. Stein [36]. Recent work by Hamza et al. [61], based on a reappraisal of global heat flow database and with due emphasis on observational data, has concluded that global conductive heat loss falls in the range from 29 to 34 TW. This is nearly 23% to 34% less than the previous estimates. We recall that geochemical constraints discussed recently by Hofmeister and Criss [8] point to the need for downsizing the current estimates of global heat loss.

The main conclusions of the present work may be summarized as follows.(1)The new VBA model of the oceanic lithosphere allows incorporation of the thermal effects of variable heat input into its basal parts, whereas previous models (HSC and Plate) take into account only effects produced by surface heat loss.(2)The width of magma injection zone in the spreading center in VBA model is relatively narrower, and the transition to stable nonmagmatic configuration takes place on time-scales much shorter than those predicted by the conventional boundary layer theory.(3)The constant temperature Plate model envisaged by Mckenzie [2] is a particular case of the more general class of VBA models.(4)The VBA Model provides a vastly improved fit to experimental heat flow data for both the younger as well as the older segments of the oceanic crust.(5)The relation for ocean floor bathymetry derived from VBA model provides an equally good fit to observational data as that provided by the hybrid model curves of C. A. Stein and S. Stein [36].(6)The VBA model fit to bathymetry data is valid for the entire age range of the oceanic lithosphere. There is no need to introduce ad hoc adjustments (artificial changes in heat flow) in model fits for ocean floor bathymetry.(7)The fits of VBA model for sea floor heat flow and bathymetry have been achieved without introducing artificial discontinuities in the temperature field of the lithosphere.(8)Agreement of the VBA model with the observables exists for any reasonable choice of input parameters. The best agreement is obtained for values closest to those believed representative of the lithosphere, particularly the lowermost extent of the lithosphere.(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., [61]).(11)The current estimates of global heat loss need to be downsized by at least 25%, in support of recent assessments [8, 37, 61] and the view that the Earth is in quasisteady thermal state.

Appendix

A. Thermal Model of Lithosphere with Variable Magma Accretion in Its Basal Parts

A.1. Solution by the Integral Transform Method

In addressing the problem described by (3), (4a), (4b), (4c) and (4d) of the main body of this paper we introduce the dimensionless variables 𝜃,𝑋,𝑍,𝐺 and 𝜏:𝑇𝜃=𝑇_𝑆𝑥;𝑋=𝐿𝑧;𝑍=𝐿𝐿;𝐺=²𝑔𝜆𝑇_𝑆,𝜏=𝜆𝑡𝜌𝐶_𝑃𝐿². (A.1) This allows us to rewrite the differential equation (3) for the system of n discretized elements as𝜕𝜃𝜕𝜏+Pe(𝑋)𝜕𝜃=𝜕𝜕𝑋²𝜃𝜕𝑋²+𝜕²𝜃𝜕𝑍²𝑋+𝐺,_𝑗<𝑋<𝑋_𝑗+1,𝑡>0,0<𝑍<1, (A.2) where 𝑗 is the discretization index, 𝑋₀>0 and Pe is the Peclet number given byPe[𝐿(𝑋)]=𝑣𝜌𝐶_𝑃𝐿(𝑋)𝜆. (A.3) The interval [𝑋_𝑗,𝑋_𝑗+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 𝐺 is negligible. In this case the differential equation and the boundary conditions become:Pe𝜕𝜃=𝜕𝜕𝑋²𝜃𝜕𝑍²+𝜕²𝜃𝜕𝑋²,𝑋_𝑗<𝑋<𝑋_𝑗+1,0<𝑍<1,𝜃(𝑋_𝑗>0,𝑍)=𝜃(𝑍),𝜃(𝑋,𝑍=0)=1,𝜃(𝑋,𝑍=1)=0,𝜃(𝑋=0,𝑍)=1. (A.4a)(A.4b)(A.4c)(A.4d)(A.4e)

The purpose of condition (A.4b) is to avoid the well-known problem of singularity at the position 𝑋=0.

We assume that 𝜃(𝑋,𝑍) may be expressed as𝜃(𝑋,𝑍)=𝑢(𝑋,𝑍)+𝑘(𝑍). (A.5)

The problem in 𝐾(𝑍) is𝑑²𝑘𝜕𝑍²=0,0<𝑍<1,𝑘(𝑍=0)=1,𝑘(𝑍=1)=0. (A.6a)(A.6b)(A.6c)

Hence the solution of problem in 𝐾(𝑍) is𝑘(𝑍)=(1−𝑍). (A.6d)

The problem in 𝑢(𝑋,𝑍) isPe𝜕𝑢=𝜕𝜕𝑋²𝑢𝜕𝑋²+𝜕²𝑢𝜕𝑍²,𝑋_𝑗<𝑋<𝑋_𝑗+1,0<𝑍<1,𝑢(𝑋=0,𝑍)=𝑍. (A.7a)(A.7b)

The condition (A.7b) is necessary for the solution 𝑢(𝑋,𝑍) to be compatible with the condition (A.4e).

We now admit that the solution of the problem in 𝑢(𝑋,𝑍) may be expressed as𝑢(𝑋,𝑍)=^∞_𝑖=0𝜓_𝑖(𝜇_𝑖,𝑍)⋅𝐶_𝑖(𝑋). (A.8) The eigen functions 𝜓_𝑖(𝜇_𝑖,𝑍) of (A.8) are associated with the following eigen value problem:𝑑²𝜓_𝑖𝑑𝑍²+𝜇²𝜓_𝑖(𝜇_𝑖𝜓,𝑍)=0,_𝑖(𝜇_𝑖𝜓,𝑍=0)=0,_𝑖(𝜇_𝑖,𝑍=1)=0. (A.9a)(A.9b)(A.9c) The auxiliary problem presented is the eigen value problem, typical of Sturm-Liouville, which has the following properties:(a)the eigen values 𝜇_𝑖 are real and positive and the order of values is such that 𝜇₀<𝜇₁<𝜇₂<𝜇₃⋯<𝜇_𝑖<𝜇_𝑖+1, where 𝑖=0,1,2,3….;(b)the Eigen functions 𝜓_𝑖(𝜇_𝑖,𝑋) associated with the Eigen values 𝜇_i are orthogonal.

Solving the auxiliary problem we have𝜓_𝑖(𝜇_𝑖,𝑍)=sin(𝜇_𝑖𝜇,𝑍),_𝑖=𝑖𝜋. (A.9d)(A.9e)

The coefficients of the expansion 𝐶_𝑖(𝑋) are obtained by multiplying (A.8) by the following operator:¹₀𝜓_𝑗(𝜇_𝑗,𝑍)𝑑𝑍, (A.10)

and we have¹₀𝑢(𝑋,𝑍)𝜓_𝑗(𝜇_𝑗,𝑍)𝑑𝑍=𝐶_𝑖(𝑋)^∞_𝑖=1¹₀𝜓_𝑗(𝜇_𝑗,𝑍)𝜓_𝑖(𝜇_𝑖,𝑍)𝑑𝑍. (A.11)

The integral on the RHS of (11) is zero for 𝑖≠𝑗 (eigen functions are orthogonal, see [40]). For 𝑖=𝑗, the integral leads to the norm 𝑁_𝑖, associated with the problem:𝑁_𝑖(𝜇_𝑖)=¹₀𝜓_𝑖²(𝜇_𝑖,𝑍)𝑑𝑍. (A.12) Consequently, the unknown 𝐶_𝑖(𝑋) is given by𝐶_𝑖1(𝑋)=𝑁_𝑖¹₀𝑢(𝑋,𝑍)𝜓_𝑖(𝜇_𝑖𝐶,𝑍)𝑑𝑍,_𝑖(𝑋)=𝑢_𝑖(𝑋)𝑁_𝑖(𝜇_𝑖). (A.13)(A.14)

It is obvious that the transform and its inverse in (A.14) are𝑢_𝑖(𝑋)=¹₀𝑢(𝑋,𝑍)𝜓_𝑖(𝜇_𝑖,𝑍)𝑑𝑍,𝑢(𝑋,𝑍)=^∞_𝑖=1𝑢_𝑖(𝑋)𝜓_𝑖(𝜇_𝑖,𝑍)𝑁_𝑖(𝜇_𝑖). (A.15)(A.16) This procedure transforms the partial differential equation into a system of ordinary differential equations. We adopt the following strategy(a)Multiply the original equation by the operator: ∫¹₀𝜓_𝑖(𝜇_𝑖,𝑍)⋅𝑑𝑍,Pe¹₀𝜕𝑢𝜓𝜕𝑋^𝑖(𝜇_𝑖,𝑍)𝑑𝑍=¹₀𝜕²𝑢𝜕𝑍²𝜓_𝑖(𝜇_𝑖+,𝑍)𝑑𝑍₁₀𝜕²𝑢𝜕𝑋²𝜓_𝑖(𝜇_𝑖,𝑍)𝑑𝑍. (A.17) By the Leibniz rule,𝑑Pe𝑑𝑋¹₀𝑢(𝑋,𝑍)𝜓_𝑖(𝜇_𝑖=𝑑,𝑍)𝑑𝑍₂𝑑𝑋²¹₀𝑢(𝑋,𝑍)𝜓_𝑖(𝜇_𝑖,𝑍)𝑑𝑍+¹₀𝜕²𝑢𝜕𝑍²𝜓_𝑖(𝜇_𝑖,𝑍)𝑑𝑍. (A.18) Introducing (A.15) in (A.18),𝑑Pe𝑢_𝑖=𝑑𝑑𝑋²𝑢_𝑖𝑑𝑋²+¹₀𝜕²𝑢𝜕𝑍²𝜓_𝑖(𝜇_𝑖,𝑍)𝑑𝑍. (A.19)(b)Multiplying the equation of the auxiliary problem by the operator ∫¹₀𝑢(𝑋,𝑍)𝑑𝑍 leads to¹₀𝑑𝑢(𝑋,𝑍)²𝜓_𝑖𝑑𝑋²𝑑𝑍=−𝜇_𝑖²¹₀𝑢(𝑋,𝑍)⋅𝜓_𝑖(𝑍)𝑑𝑍. (A.20) Equation (A.20) becomes¹₀𝑑𝑢(𝑋,𝑍)²𝜓_𝑖𝑑𝑋²𝑑𝑍=−𝜇_𝑖²𝑢_𝑖(𝑋). (A.21)(c)Subtracting (A.21) from (A.19),𝑑²𝑢_𝑖𝑑𝑋²𝑑−Pe𝑢_𝑖𝑑𝑋−𝜇_𝑖²𝑢_𝑖(𝑋)=¹₀𝑑[𝑢(𝑋,𝑍)²𝜓_𝑖𝑑𝑍²−𝜕²𝑢𝜕𝑍²𝜓_𝑖(𝜇_𝑖,𝑍)]𝑑𝑍. (A.22)(d)Developing the integrals on the RHS,𝑑²𝑢_𝑖𝑑𝑋²𝑑−Pe𝑢_𝑖𝑑𝑋−𝜇_𝑖²𝑢_𝑖(𝑋)=[𝑢(𝑋,1)𝑑𝜓_𝑖|𝑑𝑍_𝑍=1−𝑢(𝑋,0)𝑑𝜓_𝑖|𝑑𝑍_𝑍=0−𝜓_𝑖(𝜇_𝑖,1)𝜕𝑢|𝜕𝑍_𝑍=1+𝜓_𝑖(𝜇_𝑖,0)𝜕𝑢|𝜕𝑍_𝑍=0]. (A.23)(e)To conclude we use the boundary conditions in (A.23) which lead to𝑑²𝑢_𝑖𝑑𝑋²𝑑−Pe𝑢_𝑖𝑑𝑋−𝜇_𝑖²𝑢_𝑖(𝑋)=0. (A.24)

In obtaining the solution of (A.24) it is necessary to make use of the transformed boundary conditions:𝑢_𝑖(𝑋=𝑋_𝑛)=¹₀𝑢(𝑋=𝑋_𝑗,𝑍)𝜓_𝑖(𝑍)𝑑𝑍,𝑢_𝑖(𝑋=𝑋_𝑗)=¹₀𝑍𝜓_𝑖(𝑍)𝑑𝑍,𝑢_𝑖(𝑋=𝑋_𝑗)=𝑓_𝑖, (A.25a)(A.25b)(A.25c) where𝑓_𝑖=¹₀𝑍⋅𝜓_𝑖(𝑍)𝑑. (A.25d) The solution of (A.24) is𝑢_𝑖(𝑋)=𝐴_𝑖exp(𝑎₁𝑋)+𝐵_𝑖exp(𝑎₂𝑋), (A.26a) where:𝑎₁=(Pe−Pe²−4𝜇²_𝑖)2,𝑎₂=(Pe+Pe²−4𝜇²_𝑖)2. (A.26b)(A.26c) Obviously, only the solution (A.26b) has physical meaning; hence𝑓_𝑖=𝐴_𝑖. (A.26d)Consequently the problem in 𝑢_𝑖 is:𝑢_𝑖(𝑋)=𝑓_𝑖exp(𝑎₁𝑋), (A.27) and using the inversion formula𝑢(𝑋,𝑍)=^∞_𝑖=1𝑓_𝑖exp(𝑎₁𝑋)𝜓_𝑖(𝜇_𝑖,𝑍)𝑁_𝑖(𝜇_𝑖). (A.28)

In terms of the dimensionless variables the equations for temperature (𝜃), temperature gradient (𝜕𝜃/𝜕𝑧), and heat flux (𝑞) may be expressed as𝜃(𝑋,𝑍)=(1−𝑍)+^∞_𝑖=1𝑓_𝑖exp(𝑎₁𝑋)𝜓_𝑖(𝜇_𝑖,𝑍)𝑁_𝑖(𝜇_𝑖),𝜕𝜃(𝑋,𝑍)𝜕𝑍=(−1)+^∞_𝑖=1𝑓_𝑖exp(𝑎₁𝑋)𝑁_𝑖(𝜇_𝑖)𝑑𝜓_𝑖𝑑𝑍,𝑞(𝑋,𝑍)=−𝜆𝜕𝜃(𝑋,𝑍)𝜕𝑍=−𝜆[(−1)+^∞_𝑖=1𝑓_𝑖exp(𝑎₁𝑋)𝑁_𝑖(𝜇_𝑖)𝑑𝜓_𝑖𝑑𝑍]. (A.29)(A.30)(A.31)

Acknowledgments

The authors thank Professor Anne Hofmeister (Washington University, MO, USA) for critical comments and suggestions which have contributed to significant improvements in the manuscript. The present work was carried out as part of a research project for investigating thermal isostasy in southeast Brazil, with funding from Research Foundation of the State of Rio Janeiro—FAPERJ (Grant no. E-26/100.623/2007), under the program “Scientist of the State of Rio de Janeiro”. The second author has been a recipient of a Ph.D. scholarship granted by Coordenadoria de Aperfeiçoamento de Pessoal de Nível Superior—CAPES, Brazil. The third author contributed to the development of integral transform solution for variable basal heat input, as part of his Ph.D. thesis work on thermal models of continental lithosphere.