Small-scale convection beneath oceans and continents

Small-scale convection supplies heat flow of ∼17 mW m−2 to the base of stable continents where xenolith studies resolve the geotherm. However, effects of small-scale convection are difficult to resolve in ocean basins. On first pass, most seafloor appears to subside to an asymptote compatible with ∼40 mW m−2 convective heat flow. These common regions are tracked by hotspots so uplift associated with ponded mantle material is an attractive alternative. Unaffected seafloor in the North and South Atlantic continues to subside with the square root of age as expected from pure conduction. The theory of stagnant-lid convection provides good scaling relationships for heat flow. For linear viscosity, heat flow is proportional to the underlying “half-space” viscosity to the −1/3 power and the temperature to change viscosity by a factor of e to the 4/3 power. The formalism is easily modified to represent convection beneath a lid of highly viscous and buoyant cratonal lithosphere and to represent transient convection beneath thickening oceanic lithosphere. Asthenospheric mantle with linear, strongly temperature-dependent, and weakly depth-dependent viscosity is compatible with both oceanic and continental data. More complicated rheology may allow vigorous small-scale convection under most but not all old ocean basins. Still viable hypotheses require poorly understood global features, including lateral variations of asthenospheric temperature. Seismological studies have the potential to resolve the lithosphere-asthenosphere boundary, including local variations of its depth associated with small-scale convection.


Introduction
The construct of lithospheric plates provided a kinematic basis for understanding global tectonic processes in the 1960s. Geodynamic work quickly followed. Early modelers confirmed the reality of seafloor spreading from the relationship between ocean depth and seafloor age [1][2][3][4][5]. These scientists concentrated on young and middle aged (<100 Ma) seafloor where subsidence over plate age is obvious. The physics and kinematics are quite simple. Hot material rises at the ridge axis. The lithospheric material spreads laterally cooling from top down. The cool material thermally contracts causing it to gradually subside as it moves away from the ridge axis. Kinematic analytical thermal models of oceanic crust provide straightforward predictions.
The bases of plates received initial little attention because of their inaccessibility and the lack of obvious implications to surface phenomena. There was a tendency to associate the pre-plate (and modern) concept of a low velocity (S-wave) seismic zone beneath the lithosphere with a thin low-viscosity asthenosphere. For example, Figure 1 of Isacks, Oliver, and Sykes [6] shows the return flow from trenches to ridges occurring within this thin layer and slabs confined to the upper mantle. The need to consider the dynamics of the basal lithosphere first arose with respect to effective bottom boundary condition for thermal models of the oceanic lithosphere. This topic and the nature of small-scale convection near the lithosphere-asthenosphere boundary in general are the subjects of this paper.
The extent of which small-scale convection supplies heat to the base of the lithosphere was not evident to the early modelers and is still not resolved. Two classes of bottom boundary conditions were initially imposed mainly for numerical convenience. (1) Classical plate models had an isothermal boundary representing the base of the lithosphere at a given depth on the inference that small-scale convection is vigorous. The model temperature evolves to a linear thermal gradient at large distances and ages away from the ridge axis. Parsons and Sclater [7] recognized that this simple boundary condition has no obvious physical significance. Modern "CHABLIS" versions of plate models use a heat flux boundary that moves with the thickening the base of the lithosphere to represent the net effect of small-scale convection [8][9][10][11]. However, Niu et al. [12] noted that the hydrous mineral amphibole is unstable at asthenospheric temperatures deeper than ~90 km depth and that the melt-bearing mantle below 90 km depth may convect vigorously. The gross effect is similar to a classic plate model. (2) Small-scale convection supplies little heat to the lithosphere. There was hence no need for a bottom boundary condition and the plate was represented as a half space. The lithosphere continued to thicken by conduction over time with the square root of age [5]. Extant seafloor is all less than ~180 Ma so that conductive cooling penetrates to ~150 km depth, which is small enough compared to the radius of the Earth that spherical geometry and the finite thickness of the mantle can be ignored.
Early modelers also recognized that most of the extension between lithospheric plates occurs in a narrow (few kilometer) axial zone rather than the broad zone implied by the original sense of the term seafloor spreading. It is hence straightforward to mathematically account for the horizontal conduction of heat by imposing an energy conserving boundary condition at the plane of the ridge axis (Appendix A). Half-space and plate boundary conditions then give indistinguishable results at young and intermediate (<70 Ma age) lithosphere. The Fourier series formulation for a classical plate model is simpler to evaluate numerically than the Fourier integral in a half-space model. A simple way to evaluate a two-dimensional half-space model on a computer is thus to use the plate model with a very deep base [4,13].
As discussed in section 2, the predictions of plate and half-space models differ significantly for old lithosphere. The seafloor depth approaches an asymptote in the plate model but continues to increase with the square root of plate age in the half space model. It is still unresolved whether plate or half-space models better represent the actual physical processes that control the depth to the base of oceanic lithosphere. Physically, does small-scale convection supply significant heat at the base of oceanic lithosphere?
The base of the continental lithosphere is a related issue that was poorly resolved at the advent of plate tectonics and also received little initial attention. Early researchers applied the classical plate model [14,15] to explain the subsidence of passive margin and continental interior basins. These efforts obtained a thickness similar to that of plate models for the oceanic crust ~100 km. In contrast, the tectosphere model of Jordan [16] included chemically buoyant lithosphere as thick as 400 km. Little progress was made until petrologists studied xenoliths from diamond pipes (kimberlites in the general sense) (Figure 1). The general consensus is that the continental lithosphere in stable regions is ~200 km thick. Ancient (often >3 Ga) chemically buoyant lithosphere underlies cratons. Ordinary mantle forms the lithosphere of platforms. It is modestly thinner than cratonal lithosphere. Modern continental basin models involve ~200 km thick platform lithosphere [17,18].
Notably, the persistence of ~200 km thick continental lithosphere over billion-year time scales requires the addition of heat from below. With conduction alone, the lithosphere would thicken over billions of years to well over 400 km, like the lithosphere of the Moon [19]. The tendency of many regions of older seafloor to subside less rapidly than predicted by the half-space model also requires the addition of heat flow below. Small-scale convection and mantle plumes are obvious candidates. The net heat flow added by plumes is small because ponded plume material has a stably stratified base that suppresses subsequent stagnant-lid convection [17,20]. Drag from the lateral movement of plates  [60]. The thick shallow parts of the arrays are constrained by data. The curved geotherm between the arrays and the MORB adiabat is not resolved. The chemical base of the lithosphere in the Wajrakarur region of India [62] provides an upper limit for the rheological boundary layer. Its geotherm is similar to the cratonal one. Modified after Sleep [102]. that removes the basal boundary layer is an inadequate source of heat flow [19].
This paper compares small-scale convection at the bases of continental and oceanic plates with the goal of obtaining more constraints than accrue by modeling the processes in isolation. I include lines of evidence that with improved data potentially bear on this issue. I begin with a review of kinematic half-space and plate models for ocean lithosphere. I then review small-scale convection with the intent of constraining the rheology of the uppermost mantle with oceanic and continental data. The main new contribution is to apply simple relevant scaling relationships.

Kinematic lithosphere theory
Oceanic crust subsides as the lithosphere cools ( Figure 2). This process depends mostly on conduction of heat through the surface of the cooling lithosphere. Convection at its base replenishes only a fraction of the heat lost from its top until at least the lithosphere is ~100 Ma. It is useful to approximate heat transfer with conductive cooling of an initially hot half-space in a kinematic model, as the horizontal scale of plates is much greater than their thickness. This exercise provides purely conductive predictions to compare with observations. The deviations of data from these predictions quantify the effects of small-scale convection and plumes.
These conductive simplifications apply as shown in Appendix A for lithosphere younger than ~70 Ma that is a few kilometers away from the ridge axis. I begin with this case Ocean are plotted as a function of plate age [42]. The data from the South American plate are systematically deeper than the data from the African plate. The classic subsidence line with an axial depth of 2500 m and a subsidence coefficient of 350 m provides a good fit to the South American data (purple dotted line). The green (thick dotted) line is parallel to this line but shallower represents the older African data. Solid lines are from statistical fitting by Hayes [42] and differ little from my eyeball lines. Note that data younger than 1 Ma are not utilized.
to put the main subject of this paper in context. A rheological boundary layer with some convection and some conduction exists between the rigid overlying lithosphere and the adiabatic underlying asthenosphere. This layer supplies a finite but non-obvious amount of heat to the overlying oceanic lithosphere. I discuss relevant observations in sections 2.3, 2.4, and 2.5 and dynamic issues in sections 3, 4, 5, and 6.
Kinematic models for cooling of the lithosphere combine the construct of isostasy with simple one-dimensional or two-dimensional heat transfer. The formalism is intended to facilitate explanation and provide quick approximations, rather than being exact mathematically or physically. In this regard, Hofmeister and Criss [21][22][23] claimed that traditional depth-age and heat flow formalism is in gross (factor of >2) error. I give some attention to mathematical approximations, as the works of Hofmeister and Criss [21][22][23] have been widely cited. It is often productive to reexamine whether putative second order terms are really second order. Conversely, it is inappropriate to sow seeds of doubt by qualitatively claiming that the existence of ignored second order terms renders a construct invalid. Overall, one can always include small second-order terms in numerical models, but they obscure the implications of scaling relationships, making it difficult to see the effect of parameters and to generalize from numerical results. Moreover, scaling relationships allow one to omit negligible terms in calculations. Without scaling and exclusion of negligible effects, one has to construct a vast number of models to blindly explore parameter space. I alert when second order effects are ignored in derivations. Productively, the detailed nature of isostasy is important to the subsidence of old oceanic lithosphere as discussed in Section 5.4.

Lithospheric isostasy
I begin with isostasy as it provides first order predictions and hence a reference with which to compare observations. It suffices for discussion to present the equation of motion or equivalently the momentum equation in one vertical ( z , depth below sea level positive downward) and one horizontal x Cartesian dimension. Inertial forces are negligible for plate tectonics yielding and 0 xz zz where  is the deviatoric stress tensor, P is scalar pressure,  is density, and g is the acceleration of gravity. Isostasy in its simplest form arises from (1b) and is a restatement of Archimedes principle. The asthenosphere beneath the lithosphere is assumed to have laterally invariant density and no strength so that the deviatoric stress is zero and the pressure is lateral invariant (strictly laterally invariant on a gravitational equipotential surface). The lithosphere is assumed to have no strength on vertical planes so that the first term in (1b) is 0 or equivalently (1b) is laterally averaged over a scale greater than the thickness of the lithosphere to greatly diminish that term. Equation (1b) is then integrated with depth. The second term in (1b) integrates to zero since the vertical stress is zero at the free surface and by assumption deviatoric stresses are zero in the asthenosphere. These assumptions yield the well-known expression for constant hydrostatic (here lithostatic) pressure the com- where E is elevation above sea level. This equation needs to be integrated from the sea surface and include the density of water in marine settings. The approximate statement of isostasy that there is equal mass in each unit area column is obtained by making g constant. This simplification ignores the effects that arise from the spherical shape of the Earth, which are negligible for terrestrial lithospheric isostasy [24]. These include the direct geometrical effect, self-gravity of topography and its compensating root, the variation of global gravity with depth, and the tendency of the curved lithosphere to support vertical stresses as an arch. More caveats are necessary as causal use of approximate equations for isostasy can lead to confusion and even the misstatements by Hofmeister and Criss [22] that thermal isostasy is grossly incorrect. The depth to an equipotential and equal pressure model surface Z I is taken beneath lateral density variations and deep enough that the asthenosphere behaves as a fluid. It is immaterial if the model compensation depth is taken to be somewhat below the maximum depth of lithospheric density variations, as equal masses per area are then added to each column. Asthenospheric material moves vertically across this imaginary surface when shallow density changes, so Z I is an Eulerian coordinate. That is as expected; the mathematical surface at the somewhat arbitrary compensation depth Z I is not a physical barrier. For example, the rock surface moves upward when an ice sheet melts and material moves upward across Z I . There is no global nonconservation of mass as the melted water displaces the ocean bottom and its underlying mantle downward elsewhere. There is certainly no implied conversion of water to rock as claimed by Hofmeister and Criss [22,23].

Oceanic lithosphere formalism
The depth-age relation for oceanic crust is a salient feature of the Earth's surface. For this application, it is easiest to apply (2) in a differential sense by comparing different columns of material through the lithosphere, as one does not then have to precisely know the absolute density. One takes the density of material that ascends adiabatically at the ridge axis (MORB abiabat) as a convenient reference. Fast ridge axes are close to this limit. The linearized density is then where  0 is the density on the ridge axis,  is the volume thermal expansion coefficient, and T is the temperature deficient below the MORB adiabat. My derivation retains terms of the order of T in (1b) and ignores (T) 2 and higher terms. The volume thermal contraction coefficient of the upper mantle is ~310 5 K 1 and the temperature at the base of the lithosphere is ~1350°C so the maximum T is ~4% and (T) 2 is 0.16%. The linear version of (2) with these assumptions is where the depth D is referenced to the ridge axis. Solving yields that d .
This expression assumes that the effective density of the topography is that of the mantle or equivalently the thickness and density of the oceanic crust does not change away from axis. There is thus a slight error in the axial region where the temperature of the crust changes significantly. At this point, the load of water is ignored, but will be included later in this section.
The predicted subsidence is kinematically the same as would occur if all the thermal contraction occurred in the vertical direction. This situation would occur if the seafloor surface was perfectly rigid and all the underlying material behaved as a fluid. On this point, Hofmeister and Criss [23] state that horizontal thermal contraction cannot occur because the plate is one-dimensional, so  is (5) should be replaced by /3. This grossly incorrect statement ignores both elastic and inelastic deformation. Appendix B addresses the details.
Classical treatments of depth-age relationship of oceanic lithosphere use the linearized integral (5). Further simplification accrues if one assumes that physical properties including the thermal expansion coefficient are constant with temperature and depth. A natural temperature boundary condition of ~0°C exists at the seafloor and provides the origin the depth coordinate z f . A Lagrangian coordinate system that moves with the plate is convenient off axis. Such models include only the vertical conduction of heat. If one wishes to model the temperature and lateral conduction of heat in the axial region, an Eulerian coordinate system is conveniently centered at the ridge axis (Appendix A).
I begin with the well-known half-space conduction formulas. The half-space temperature as a function of depth is where T is relative to the surface temperature, T L is the half space temperature of material at the ridge axis, z f is depth from the seafloor,  is thermal diffusivity, and t L is lithospheric age [25]. The half-space heat flow as a function of depth is where k is thermal conductivity, and q S is the surface heat flow [25]. This treatment is partly Lagrangian in the sense that it follows a place on the seafloor and Eulerian at depth as some geometrical distortion occurs as vertical thermal contraction moves material at depth toward the seafloor. The surface area of a column of material is also changed by horizontal thermal contraction but the relative error is bounded by the order of T.
Continuing with thermal isostasy, (5) is evaluated by integrating the rapidly convergent (6) with respect to depth to infinity m m w 4 , π where the term in the bracket corrects for the weight of the seawater,  m  0 is the density of the compensating asthenosphere, and  w is the density of water [25]. This treatment ignores the small direct effect of the water pressure on the seafloor on the density of the underlying rock. The seafloor subsides as heat and its associated buoyancy escapes from the top of the lithosphere. Differentiating (8) and noting that  k/C where C is specific heat per mass yields the O'Connell relationship when no heat is supplied from below where the heat flow to the base of the plate is q B . This relationship is convenient in that only heat flow into and from the lithosphere are involved. In particular, one does not need to know the thermal conductivity within the lithosphere [26]. Conversely, the ratio of the thermal expansion coefficient to specific heat is critical and involves depth and temperature averaged quantities in the lithosphere. The relationship thus does not apply simply when phase changes drive subsidence. Modern CHABLIS plate models are based on (9b) [11]. To compute temperature within the lithosphere, the basal heat flow is mathematically supplied at the base of the lithosphere at a geotherm that represents the top of the rheological boundary layer. (Applying a heat flux at a constant depth would unphysically increase model temperatures above the mantle adiabat beneath young lithosphere that has not yet thickened to that depth.) Stagnant-lid convection cannot increase the boundary layer temperature above the adiabat of the upwelling material.

Summary of oceanic lithosphere data
On face value, it appears that measurement of both surface heat flow and subsidence rate would constrain the heat flow from below from (9b). The situation on the Earth is complicated as accurate data on the heat flow from the mantle are unavailable. Rather, it became evident early in the history of plate tectonics [7] that: (1) The central rift at slow ridges and the axial high at fast ridges obscure thermal subsidence within tens of kilometers of the ridge axis. It is also rarely possible to measure conductive heat flow close to the axis. Traditional heat flow probes penetrate into sediment. In most localities, there has not been enough time for sediment to accumulate on very young crust. The details of the ridge axis do not significantly affect the behavior of old lithosphere as shown in Appendix A. (2) Equation (8) reasonably predicts subsidence out to at least ~70 Ma crust.
Heat flow values are highly scattered with a mean systematically below that predicted by (7). The scatter and deviation are greatest near the ridge axis and decrease to a small value for crust older than ~40 Ma. (3) Old oceanic crust is quite rare because of subduction. Most old crust resides along the passive margins of the Atlantic and in the northwest Pacific. The mean seafloor depth (and probably the mean heat flow) trend to asymptotes rather than continuing to increase (decrease) with the square root of age. Subsequent compilations have confirmed these inferences [27][28][29][30].
I concentrate discussion in this paper on bathymetry. Seismic velocity at the asthenosphere-lithosphere transition [31] and closely spaced parallel hotspot tracks [32,33] provide constraints that I discuss in subsection 5.3 and section 6 after discussing physics of convection. Gravity data confirm that the systematic relationship between depth and plate age as well as regional variations in depth result from density differences at lithospheric and asthenospheric depths [29,34]. Such data, however, provide little resolution on the lithosphere-asthenosphere transition. Moreover, lateral density heterogeneities below the asthenosphere dominate the gravity field and dynamic topography at wavelengths comparable to plate dimensions making it very hard to resolve shallow effects [35]. Interpretation of gravity data is beyond the scope of this review for brevity.
The scatter and deviation of measured conductive heat flow from conductive predictions on young seafloor are traditionally attributed to hydrothermal circulation through outcrops and thin sediments. I digress as Hofmeister and Criss [21,22] objected to this interpretation. The hydrothermal circulation is shallow (upper few kilometers) so it does not significantly affect the temperature at lithospheric depths and hence thermal contraction away from the ridge axis. It is certain that off-axis hydrothermal circulation occurs but it is, however, difficult to obtain its precise magnitude independently of subsidence. The argument for its importance is hence phenomenological. Following Stein and Stein [28], the measured heat flow is closer to predictions when it is measured above thick impermeable sedimentary cover that suppresses venting hydrothermal fluids directly into the ocean. In that case, heat flows through the sediments by conduction. Hydrothermal circulation in the igneous oceanic crust causes the temperature at the base of the sediments and the heat flow through the sediments to vary laterally. This effect causes scatter but not systematic bias in the measured heat flow. Hydrothermal precipitation within the igneous crust suppresses hydrothermal circulation by reducing permeability as the oceanic crust ages, greatly reducing the bias and scatter in measured heat flow. The measured and predicted heat flows are similar in the 50-70 Ma age range, where plate and half-space models yield similar predictions. Applying (7) and (8) thus calibrates (9a) well enough to exclude gross errors in the formalism but not well enough to precisely estimate basal heat flow.

Classes of data on seafloor subsidence
The deviation of subsidence toward an asymptote is of relevance to small-scale convection. Its existence or nonexistence of an asymptote has been widely argued since the works of Davis and Lister [5] and Parsons and Sclater [7]; Stein and Stein [27] and Zhong et al. [31] review the discussions and provide additional references. Goutorbe [36] presented a combined inversion of bathymetric, seismic, and heat flow data which support a CHABLIS plate model. My conclusions are closest to those of Korenaga and Korenaga [35].
Relevant data are already in hand in the sense that plate age [37] and sediment-corrected depths are known for virtually all of the Earth's seafloor. It is highly unlikely that oceanographers will find unknown in-place seafloor significantly older than 180 Ma. The subsidence rate from (8) is far too slow to monitor over even a century.
Paleodepth data from paleontological studies of sediment cores are currently woefully inadequate for resolving sub-sidence of old oceanic crust [38]. For example, Cliff [39] analyzed data provided by deep sea drilling. These compiled data partition depths to 0-150 m, 150-500 m, 500-2000 m, 2000-4000 m, and deeper. This irksome compilation method renders deep ocean data totally worthless right where one needs information on asymptotic subsidence. It also gives the spurious impression that the ocean depth was exactly 4000 m, when the compiled depth moved from 2000-4000 m into the deepest bin. Near-sea-level erosion surfaces with the current depths provide reliable information, but only for hotspot ridges. However, it is then critical to know the seafloor age at the time of beveling. For example, not knowing whether the rock surface was beveled at zero crustal age or 4 Ma would cause an uncertainty of 700 m in predicted subsidence. (I use a round number version of (8) where the subsidence in the first million years is 350 m [7] and perfect square ages in my quick examples to aid in following calculations.) Careful paleodepth studies are certainly warranted.
Thus, most modern treatments of depth-age relations use current crustal ages and seafloor depths corrected for sediment load. One then obtains apparent subsidence as in Figure 2. Hotspots have several effects. Plumes directly cause transient uplift with thermal buoyancy. The oceanic crust of on-axis hotspots is thicker than normal crust, which results in permanent shallow depth. The chemistry of seafloor basalts, however, does provide a reliable constraint on the zeroage depth of ancient oceanic crust [40] that I apply at the end of the next subsection. Off-axis hotspots create volcanic edifices and hence shallower depths.

Representative bathymetry data and some philosophy
I illustrate issues by viewing actual data on a scatter plot rather than depending on statistics of data suites ( Figure 2). In practice, one obtains a slope as in (9a) and (9b) over a significant age range, avoiding the complicated axial region. The reference depth for zero age in (5) and (8) is thus extrapolated inward from somewhat older seafloor and the roughness of the seafloor averaged out. With regard to magnitude of the tail of the subsidence curve, asymptotic seafloor with an apparent age of 81 Ma is 1050 m shallower than half-space seafloor with an age of 144 Ma.
To explain my philosophy and semantics issues, I consider an analogy, a "friendly society", the equivalent of a modern life insurance and annuity company in medieval Europe. A quick look at death records indicates that few people survive beyond age 35 so it is attractive to consider this to be the "normal" lifespan. Further looking at records shows that most people die from barbarism, famine, and pestilence. It is a matter of personal worldview whether this situation is the normal state of affairs. One might well consider smallpox scars to be more normal than skin wrinkles. The question of lifespan in the absence of calamities might well rise. Finding a few very fortunate septuagenarians would provide much more information on maximum lifespan, than massive statistics on famine victims.
With regard to the ocean floor, asymptotic subsidence (or the lack thereof) would be obvious for 300-400 Ma crust, but in-place oceanic lithosphere older than ~180 Ma does not exist on the Earth. Thus, the situation is analogous to a medieval savant speculating on the effects of aging by studying 45 year olds. Like my savant, it is productive to consider specific mechanisms, here for adding heat to the oceanic lithosphere: mantle plumes and small-scale convection. The former like medieval smallpox is common on the Earth but can be recognized by being local and episodic. It is significant, as old oceanic lithosphere has had a considerable time to be underplated by mantle plumes and old hotspotfree oceanic crust is quite rare [31]. Small-scale convection conversely is an inevitable consequence of the aging of the oceanic lithosphere.
Returning to data, I consider the South Atlantic ( Figure  3). There are several hotspot tracks. The Rio Grande and Walvis tracks are oblique to flow lines parallel to fracture zones. A depth profile along a flow line hence crosses on and off the tracks, producing apparent uplift and subsidence. Similarly, the Southeast Indian Ridge south of Australia taps mantle that is ~50 K lower that the typical MORB value. Profiles from the axis cross from seafloor above cool mantle to seafloor above normal mantle [41].
There are numerous volcanic edifices east of the ridge axis in the South Atlantic and no obvious ones west of the ridge axis in the swath in Figure 3. The depths on the west are as expected systematically deeper than those on the east [42] (Figure 2). The subsidence slope of the older seafloor in compatible with that predicted by (8). Notably, the slope on the shallower seafloor to the east of the axis is similar, although the seafloor is overall shallower.
The North-Central Atlantic illustrates further issues ( Figure 4). A volcanic margin exists from the Bahamas into southern Canada [43]. The New England, Corner, and Great Meteor seamounts are a long-lived track. The Bahaman chain and Cape Verde may be segments of another track separated by subsequent seafloor spreading. The Azores and the Canaries are additional hotspots. Only the Sohm plain is not tracked by obvious hotspots [44]. Its northern part lies offshore of a nonvolcanic margin [43]. Its depth is that predicted by the half-space model.
Topographic raises associated with the Bermuda, Cape Verde, and Canary tracks occur on ~115 Ma crust. The northern end of Bermuda rise approaches the New England track. These bands of seafloor are shallower than both the younger and older seafloor on their flanks. This seafloor when included in compilations causes ~90-130 Ma crust to be shallower than both older and younger crust [27,29,30]. I return to the physics of this feature in subsection 5.3.
Fortunately, studies of the chemistry of oceanic crust in a few boreholes allow inference of actual as compared with apparent subsidence [40]. These workers culled data where off-axis edifices are present and obvious hotspot tracks. They then estimated zero-age depth from the chemistry of the basaltic basement. The total subsidence with respect to the computed zero-age depths is in agreement with the halfspace model at more than half the sites. Despite the culling, there is some remaining effect of near-axis hotspots. The Figure 3 Gravity map of the South Atlantic Ocean after Sandwell and Smith [103]. The data in Figure 2 come from the latitude band in the box. Most of the South American plate in this band is well away from hotspot tracks. There are several tracks on the African plate but the locations of active plumes are not obvious.

Figure 4
Gravity map of the North Atlantic Ocean after Sandwell and Smith [103]. Shallow seafloor of ~90-130 Ma is associated with the Bermuda, New England, Cape Verde and Canary rises. The Sohm Plain has depths expected from a half-space model. The rises are evident on a gravity map, which indicates that the compensating density variations are deep, that is, within the lower lithosphere or asthenosphere.
computed mantle source temperature of the ascending material was on average 50 K hotter than modern MORB material. As correctly stated by Humler et al. [40], the variation is likely a distal effect of hotspots as the cooling rate of the Earth's mantle is much slower, only ~50-100 K Ga 1 [45][46][47][48].
Petrological studies of mid-plate oceanic lavas confirm that the lithosphere thickens with age until at least 70 Ma [12]. These data, which provide snapshots at the time of eruption, are compatible with a plate model with a final thickness of ~90 km. A sampling bias, however, exists, as mantle plumes that thin the lithosphere are associated with much off-axis volcanism. One does not have lavas to sample in hotspot-free regions that continue to subside with the square root of age.

Implications to small-scale convection
Overall, most old oceanic crust is shallower than expected from a half-space subsidence model. This is evident in all statistical compilations [7,27,29,30]. The effects of hotspot tracks explain much of the discrepancy, as is evident when statistical efforts are made to cull that part of the data [35,40]. There is clearly a preservation bias as the oldest crust formed near hotspots. Specifically the Atlantic crust is associated with volcanic break up margins most likely an effect of plumes. The depth in regions with no obvious hotspots is often close to that predicted from the half-space model. That is, there is no strong evidence that lithosphere that was always distal to hotspots received significant heat flow from small-scale convection.
For quantification, it is useful to compute the total amount of heat that needs to be added to the base of the lithosphere over its evolution to cause the depth to deviate from a half-space cooling model. In my round numbers, 1 W m 2 over 1 Ma causes 350 m of elevation change. I give results for the mean deviation of heat flow over 100 Ma and 144 Ma as examples. Keeping extra digits so differences are apparent, the half-space heat flows at these crustal ages are 50 and 41.7 W m 2 respectively. The basal heat flow averaged over plate age needed to cause a 1000 m deviation in elevation is 28.6 and 19.8 W m 2 .
Conversely, observations including Figure 2 do not preclude a deviation of few hundred meters from the half-space model even in the regions that seem to follow it. This deviation implies mean basal heat flow of ~1/3 of those given in the previous paragraph. Moreover, the depth coefficient implying 350 m for 1 Ma crust, is not precisely enough constrained to predict depth better than a few hundred meters. For example, a 10% error gives a predicted 350 m difference for 100 Ma in the range suggested by Hillier [30]. For reference, the constant basal heat flow in the models of Doin and Fleitout [8] and Wei and Sandwell [26] is 38 W m 2 . Wei and Sandwell [26] retained the traditional rounded value 350 m for the depth coefficient from Parsons and Sclater [7]. Korenaga and Korenaga [35] preferred 323-336 m and resolved no basal heat flow. Hillier [30] agreed that these values best represent the data, but did not attempt to resolve mechanisms.
Finally, heat flow itself is too strongly affected by hydrothermal circulation to provide more than qualitative agreement with the inferences from depth-age. Moreover, it takes a long time for heat that is currently supplied to the base of the lithosphere to reach the surface by conduction. The current basal heat flow associated with modern ancient crust is thus unconstrained by current surface heat flow [49].

Dynamics of small-scale convection
Small-scale convection likely supplies heat to the base of both continental and oceanic lithospheres. Parameterized convection theory seeks to obtain laterally averaged properties including temperature and heat flow that one needs to compare with observations. Like the germ theory of disease and epidemics, it provides insight, here with regard to heat transfer into the base of the lithosphere. In particular, the vigor of convection depends on the properties of the actively convecting rheological boundary layer rather than the properties of the deep fluid and the rigid lithosphere. This eliminates numerous potentially free parameters from discussion. I present dimensional relationships and make limited use of two-dimensional numerical models to appraise and illustrate their implications.
The shallow part of the lithosphere acts as a lid to convection. Three cases pose the issues related to oceanic and continental lithosphere. (1) In stagnant-lid convection, the shallow lithosphere does not partake in convection because it is cooler and much more viscous than the underlying rheological boundary layer. (2) In chemical-lid convection, the lithosphere above the boundary layer does not downwell because it is chemically less dense than the underlying mantle. It is not also entrained significantly into downwellings because it more viscous than normal mantle material at the same temperature. (3) Rigid-lid convection is a useful approximation of infinite chemical-lid viscosity for the Earth and a convenient laboratory and computer situation.
A well-known formalism exists for stagnant-lid convection. I present results compiled from the work of Solomatov [50], Solomatov and Moresi [51], and Sleep [52] with modifications for chemical-lid by Sleep and Jellinek [53] in order to illustrate how simple formulas arise. A Rayleigh number approach yields the same dimensional results as a force balance approach. In both cases, the heat flow depends only on the properties of the boundary layer, not on the total thickness of the convecting region. I use physical properties rather than dimensionless variables to aid in following calculations and so that the convective heat flow through moderate age seafloor is apparent.

Rheology of the mantle
Traditional thermodynamic expressions for the rheology involve activation energy and activation volume. They apply to a single well-defined solid phase. Mantle mineralogy is complicated in that the phase assemblage and the compo-sition of each phase depend on pressure and temperature. Trace amounts of partial melts or grain-boundary phases may significantly affect rheology. In addition, there is little agreement in the literature on the thermodynamic viscosity parameters of the mantle.
Until better insight into mantle rheology is available, I treat the Earth as a natural experiment and attempt to resolve rheological parameters. To begin this treatment, I express scalar viscosity, stress, and strain rate from simple shear with a change of velocity in the x-direction of V x between two parallel planes separated by distance Z in the Z-direction. The shear stress (formally traction) on the planes is , where the spatial derivative of velocity is the strain rate ′.
For small scale convection, it is convenient to take the temperature of the asthenospheric half space T half as a reference. It thus is helpful to compact notation by making the viscosity explicitly temperature and depth dependent where  half is the viscosity of the adiabatic half-space, T  T half T is the temperature below that in the half space, T  is the temperature scale for viscosity and D  is the depth scale for viscosity. This relationship is intended to apply over a limited range of temperatures and depths centered at the base of the rheological boundary layer. In terms of activation energy, the temperature scale is approximately RT 2 abs /E*, where R is the gas constant, T abs , and E* is the effective activation energy. I use values of 40-100 K for T  in this paper. For T abs = 1600 K, the activation energy is 532.5-213.0 J mol 1 .
The strain rate in the Earth's mantle may well depend nonlinearly on stress. The tensor expression for an isotropic material is where ′ is the strain rate tensor,  is the deviatoric stress tensor, ij are tensor indices, |  | is the second invariant ij ij   of the deviatoric stress tensor here normalized so that it yields the shear stress in simple shear,  ref is a reference stress. The factor of 2 is necessay for (12) to reduce to the scalar expression (10) for linear viscosity n=1. The term apparent viscosity is in widespread use in geodynamical literature; it is here I include nonlinear rheology in the deviations because it has been associated with the depth-age curve [54,55]. I concentrate discussion on linear n=1 Newtonian viscosity for simplicity. The rule of thumb that a linear rheology with nT  is approximately equivalent to a nonlinear rheology with T  is in wide use in the geophysical literature. It often gives a tolerable representation as in the work of van Hunen et al. [54]. This approximation is, however, not strictly correct, as can be seen by substitution into the dimensional formulas derived in this in this paper.

Stagnant-lid formalism
For now, I consider purely temperature-dependent viscosity to compact notation. Intuitively and with forethought to section 4.1, the quantity T  is a scale for the temperature contrast across the rheological boundary layer. Following Solomatov and Moresi [51], I let the rheological temperature contrast be where a n is a dimensionless constant of order 1 that depends on n, the power of the rheology. The thickness of the rheological boundary layer is instantaneously Z rheo and for now free to vary in the derivation. I consider convection started from strong perturbations as opposed to spin-up from a tiny perturbation and constrain the vigor of convection by balancing body forces with viscous forces. The stress from body forces resulting from density variations is obtained by dimensionally integrating (1b) across the rheological boundary layer to obtain rheo , where  is unperturbed density, g is the acceleration of gravity, and  is the volume thermal expansion coefficient.
This expression implicitly includes the effect of the slope of the boundary layer between upwellings and dowwellings, which I consider explicitly in section 4.3. Velocity changes most rapidly in the z-direction from a finite value at the base of the boundary layer to zero at its top. Components of the strain rate tensor in (12) where the effective viscosity in the boundary layer is that at the temperature T=a n T  below the half-space adiabat. The strain rate occurs across the boundary in a region that is Z rheo thick so the (vertical) velocity in the boundary layer is dimensionally The lateral temperature contrast in the boundary layer between upwellings and downwellings is dimensionally T rheo =a n T  . The convective heat flow is thus, where C is volume specific heat.
The linear form n=1 of this equation illustrates the relationship of this approach to a Rayleigh number approach The conductive heat flow across the boundary layer is The Nusselt number, the ratio of convective to conductive heat flow, is a measure of the vigor of convection: where =k/C is the thermal diffusivity. The dimensionless term in the bracket in (20) has the form of a Rayleigh number local to the rheological boundary layer. The vigor of convection depends on this quantity.

Steady state heat flow
Continental lithosphere is to the first order in equilibrium with the heat flow into its base. Simple scaling relationships arise because convection and conduction carry comparable heat flows in the middle of the rheological boundary layer. The thermal gradient in the boundary layer is thus ~1/2 that in the overlying lithosphere. Dropping this factor yields the dimensional expression for thermal gradients: Equating the conductive heat flow through the rigid lithosphere with convective heat flow into the bottom of the lithosphere from (17) and (21) gives The thermal gradient in the lithosphere T L /Z L appears on both sides of the approximate equality in (22). Algebraically solving for this quantity and multiplying by k yields that steady state heat flow For the special case of a Newtonian rheology n = 1, eq. (23) simplifies to where the final calibrated equality comes from the works of Davaille and Jaupart [56][57][58].

Rheological boundary layer and the continental xenolith geotherm
The temperature contrast across the rheological boundary layer can potentially be inferred from mantle xenolith data ( Figure 1). As shown in this section, such data provided coupled constraints on the temperature scale T  and the power-law n.

Rheological boundary layer thickness for stagnantlid convection
I obtain the temperature contrast across the rheological boundary layer following Solomatov and Moresi [51] and Sleep and Jellinek [53]. Solomatov and Moresi [51] defined this quantity as the lateral temperature contrast across a convection cell from the apex of an upwelling to a downwelling. With this definition, the rheological boundary layer thickness is the depth range for a isotherm in the boundary layer; that is how far the rheological boundary extending up into the lithosphere. This quantity is conveniently obtained from numerical models and physical experiments ( Figure 5). We first derive the dimensionless temperature contrast parameter a n in (13) that represents rheological temperature contrast. I then discuss related temperature contrasts in the rheological boundary layer.  [20]. The viscosity in the models is calibrated so that the conductive geotherm extrapolates to the MORB adiabat (here 1300°C) at 200 km depth. Parameters are in Table 1.
Solomatov and Moresi [51] suggested that that the temperature contrast of the rheological boundary layer selforganizes so the heat flow is a maximum. This occurs when / 0    v n q a in the transient equation (17) giving, The proportionality in (25b) applies to any definition of T rheo . Solomatov and Moresi [51] obtained 1.2 in the calibrated equality in (25b) from numerical and laboratory experiments for the specific definition of T rheo as the lateral temperature contrast. These results apply for both transient and steady state convection. Note that T rheo ln(10)T  for a linear material so the viscosity within the rheological boundary layer changes by round factor of 10. Parameterized convection equations are sometimes expressed in base 10 so the temperature scale is ln(10)T  and the multiplicative constants modified accordingly.
An important but subtle point arises from this formalism. The convection pattern self-organizes toward an extremal in heat flow. The general property of extremals applies: the heat flow changes slightly if the flow pattern differs somewhat from the optimal one. In particular, wavelengths imposed by lithospheric relief in Earth or artificial side boundaries in numerical models will only slightly affect laterally averaged heat flow. Moreover, it does not matter if  the extremal concept is only a good approximation.
The dimensional derivation of (25a) does not depend on the precise definition of a n . However, attention to the definition in (25b) of T rheo is needed for application to the Earth. Its definition in terms of lateral temperature contrast is convenient for numerical and laboratory experiments. I suggest two alternative definitions of T rheo that are applicable respectively to xenolith geotherms and convection beneath a rigid chemical lid.
Xenolith geotherms provide an instantaneous estimate at the time of kimberlite eruption. They are most likely to sample upwelling regions where partial melting occurs. I present numerical models to illustrate the estimate of the rheological boundary layer thickness obtained in this way. The models retain well-constrained parameters from Nyblade and Sleep [20] and Sleep and Jellinek [53]: thermal  Figures 5 and 6 illustrate features of the boundary layer in vigorous nearly steady state stagnant-lid convection. First, the deep isotherms are "scalloped", not flat. The lateral temperature between upwelling and downwellings is a 1 T   2.4 T  as expected from (25b). However, the geotherm in upwellings (from where kimberlites are most likely to ascend) diverges from the linear conductive gradient at a temperature contrast of ~T  , not 2.4 T  as suggested by (25b). I define another rheological temperature contrast T up = c 1 T  T  that applies within sampled upwellings to distinguish it from T rheo . The next two subsections apply this effect to the transition between stagnant-lid and rigidlid convection.

Rigid-lid convection
Chemical lithosphere beneath cratons likely forms a lid on the rheological boundary layer. Niu et al. [12] discuss an analogous possibility for the oceanic lithosphere that leads to a plate thermal model. The hydrous mineral amphibole is unstable below ~90 km depth at temperatures near the MORB abiabat. The region below 90 km depth contains hydrous melt and has a much lower viscosity than shallower material. The overlying material acts a chemical lid. I do not specifically discuss this case as scaling relationships for cratonal lids apply.
For simplicity, I consider rigid-lid convection, an endmember case of chemical-lid convection. I let the base of the lid where the velocity is 0 be flat. Returning to first order physics, the stagnant-lid equations apply to a rigid lid when one replaces the rheological temperature contrast a 1 T  with the laterally averaged temperature contrast AT  beneath the rigid lid [53]. This yields a steady state formula where the exponential term in (24)

Transitional rigid-lid convection
The thickness of platform lithosphere that is presumably not chemically buoyant is modestly less than cratonal lithosphere that is chemically buoyant (Figure 1). Stable continents are thus in the transition between stagnant-lid and chemical-lid convection, so (26) does not give a straightforward prediction. That is, the attention to the limit where heat flow is modestly less stagnant-lid value is required. Sleep and Jellinek [53] noted that this transitional regime is complicated but did not attempt to explain its physics and stated only that there is a critical value of the normalized temperature contrast A crit above which the rigid lid does not affect convection and the process is effectively stagnant-lid convection.
The transitional regime involves a subtle feature of rigidlid convection: the deep isotherms are much flatter ( Figure  7) than the scalloped isotherms below a stagnant lid ( Figure  5). The dynamic formalism needs to be subtly modified. Temperature contrasts within the rheological boundary indicated by geotherm curvature in upwellings (Figure 8) are similar to those for stagnant-lid convection ( Figure 6). Thus the relevant temperature scale is still the stagnant-lid limit a 1 T  . Strain locally occurs across the boundary layer and The derivation of (24) carries through so the ratio of transitional rigid-lid heat flow to stagnant-lid heat flow at steady state for a linear material is which quantifies the statement of Sleep [59] that convection modulated by lithospheric relief is slightly more vigorous than convection beneath flat lithosphere. Sleep and Jellinek [53] used normalized coordinates from (26) A and (q/q SL ) 3/4 to illustrate time histories (Figure 9).
The normalized temperature contrast increases with time when calculations are started with a linear temperature profile to the base of the rigid lid. Evolution in the model takes over 1 Ga, so real continental lithosphere does not reach full steady state. Models with T  of 100 K approach the stagnant-lid limit at lower values of A than do models with T  of 60 K. This is expected for two reasons. (1) At steady state, the heat flow through the rigid lid is  Table 2.
where the Z lid is the depth to the base of the chemical lid. Increasing A in the transition region modestly increases the convective heat flow as in (26) The ratio of scallop height to rheological boundary layer thickness is thus reduced at high T  . That is, the coefficient  in (28) is larger at low T  .
This reasoning defines two regimes for normalized heat flow. At low values of this parameter, the lack of scalloping beneath a rigid lid as a small effect on convective heat flow. At higher normalized heat flow, the lack of scalloping causes normalized heat flow to increase slowly with normalized temperature contrasts. Equation (28) provides some indication of the transit between the regimes: the transitional normalized heat flow coordinate (to the ¾ power is  1/4 . Values of  of ~1.5 and ~2 apply for T  of 100 and 60 K, respectively. The normalized heat flow coordinate at the transition is ~0.90 and ~0.84, respectively, in crude agreement with the calculations in Figure 9.
There is a degree of unpredictability in extrapolating the rigid-lid results to real chemical lithosphere that can flow to some extent and thus have a scalloped base. The depths to the bases of the chemical lithosphere and the thermal lithosphere vary within cratons [61,62]. The trend of scalloped relief on the base of the chemical lithosphere may sometimes align with the direction of plate motion and hence the strike of the axes of convection cells. In this case, the geotherms then would be similar to those beneath a stagnant lid and the heat flow would be enhanced above that for a flat chemical lid. This geometry would develop when the direction of plate motion and hence the axes of convective rolls remained constant over a long period of time. Changes in plate motion would leave the chemical scalloping poorly oriented to aid convection.

Application to continental lithosphere
Xenolith studies provide instantaneous geotherms. These data alone do not provide an indication of the position of the base of the chemical lithosphere ( Figure 8) or even directly whether chemical lithosphere exists at a given site. Still, the curvature of the geotherm (or the lack thereof) does constrain T  and n through (25b). For this application, one needs to use geotherms underlain by normal adiabatic mantle not plume material. Ideally xenolith data should be analyzed specifically for these purposes. I select analyses from the literature that present results that I can use in a straightforward manner, rather than attempting to reanalyze a massive number of results. Figure 1 shows typical results of xenolith studies compiled by Francis and Patterson [60]. No curvature is resolved. Assuming that the kimberlites arose from MORB adiabat asthenosphere an upper limit of T up . Here this difference between the hottest xenoliths and the MORB adiabat is ~100 K for both the platform and craton geotherm. Assuming a linear rheology n=1 and that the geotherms represent upwellings as in Figure 6 and 8, T  100 K. This exercise requires that the relative calibration of the MORB adiabat and the xenolith geotherm is accurate and applicable but does not depend critically on the absolute values of the inferred temperatures.
Composition of xenoliths indicates the base of chemically buoyant lithosphere. The absolute radiometric ages of xenoliths indicate whether material resided within the lithosphere for a long period of time. I discuss examples of the former approach, which requires only the temperature contrast between the base of the chemical layer and the MORB adiabat, which for purpose of calculation is 1350°C plus 0.3°C times the depth in kilometers. This method is less sensitive to the presence of plume material as it involves a long-lasting geological boundary.
Kopylova and Caro [61] discussed the xenolith geotherm in the southern Slave Province of Canada. Depleted samples extend to ~250 km depth and ~1300°C. Karmalkar et al. [62] compiled results for Wajrakarur district of India. Depleted samples extend to ~185 km and ~1200°C. Saltzer et al. [63] studied the Kaapvaal craton of South Africa. Depleted samples extend to ~150 km depth and ~1150°C. The temperature differences between the base of the chemical lithosphere and MORB adiabat are 75 K, 156 K, and 195 K, respectively. From the discussion in subsection 4.3, geotherms are expected to be relatively flat in upwelling under a flat chemical lid so that the normalized temperature contrast measured by geotherms is ~A crit >3. So T  < 1/3 of the measured temperatures. This result applies even if the sampled regions were in the stagnant-lid regime, rather than the chemical-lid regime.
The analysis based on sections 4.1 and 4.3 of the geotherm curvature data and the chemical lithosphere data indicate that the value of T  is relatively small. Allowing for considerable imprecision, the results are marginally compatible with 100 K for the rheological temperature scale; 60 K or lower is more likely. The permitted value of the temperature scale becomes smaller if the rheology is non-linear n > 1 from (25b).
Another potential calibration arises because the Earth's mantle adiabat has slowly decreased geological time [45][46][47][48]. The heat flow beneath platforms underlain by chemically normal mantle has thus waned over time and the thickness of platform lithosphere has increased [53]. However, this process provides only a loose constraint on n in this case where T rheo in (25b) is obtained as in this subsection. We assume that we have used xenolith geotherms to measure the temperature contrast in upwellings and (25) to obtain T up =T  (n+1)/2 (calibrated by the n=1 result) and that the asthenosphere temperature decreased by T ad over some period of geological time. The steady state heat flow from the temperature dependence of viscosity (11) and (23) is the new staganant lid heat flow is proportional to the old heat flow before the temperature change where the adiabat decreased by T ad over a given time. The factor involving n varies from 1/3 for n=1 to 2/5 for n=3. The analysis in this subsection did not constrain T up to this precision. Neither is the secular temperature change T ad over a given time that well constrained.

Transient heat flow through oceanic lithosphere
It is straightforward to obtain scaling relationships for convective heat flow beneath ocean lithosphere. Initially convection is weak and the conductive geotherm in (6) is a good approximation of the actual geotherm. The thickness of the rheological boundary layer is proportional to the lithosphere thickness, which increases with the square root of age [64].

The transient heat flow equation (17) then becomes in normalized form
where t L is plate age and t s is a temperature scale proportional to the thermal time constant of steady state lithosphere with heat flow q SL , here the time for the half space heat flow to reach the steady state heat flow [49]. These equations hold when the predicted heat flow is significantly less than the steady state value. Sleep [17,52] discussed the final approach to steady state. Equation (32) has a simple calibrated form for linear rheology n=1: where the dimensionless term (of the order of 1)  is a modest function of the rheological scale T  . The value from numerical modeling is 0.27 and 0.50 for T  = 60 K and 120 K, respectively [52]. The constant  arises because the rheological boundary layer thickness scales inversely to thermal gradient within the rheological boundary layer. The deep thermal gradient is proportional to the surface thermal gradient as assumed in the derivation of (34), but decreases as the half-space temperature is approached. The thickness of the rheological boundary layer from a given surface heat flow thus is greater for small values of T  than large ones.
Sleep [59] recommended T  as an acceptable approximation by evaluating limits of the error function and appraising numerical models. A calibrated relationship is where T cal = 240 K. For larger values T  > T cal , the temperature drop across the rheological boundary layer is a significant fraction of the temperature drop across the lithosphere. Convection is not in the stagnant-lid regime. Some caveats apply to (34). It does not apply as steady state is approached, as the geotherm is significantly different than the half-space geotherm (7). Rather the actual heat flow is on average less than predicted by (34) and the thermal gradient in the rheological boundary layer is ~1/2 the conductive gradient. Sleep [17] gave calibrated relationships based on (17) for this situation. For comparison, Dumoulin et al. [64,66] consider the time of convection to become vigorous which is analogous to t s using Rayeligh numbers. Their work utilizes t s , but their equations cannot be manipulated into the forms of (34) and (35).
The transient heat flow in (32) can be applied to the base of the model lithosphere in a modified CHABLIS plate model. If n=1 and T  is small, the convective heat flow quickly approaches its steady state value so constant basal heat flux is a reasonable approximation. Furthermore the time to approach steady state is small if the final steady state heat flow is high. Conversely, a low-value of steady state heat flow and a large value of T  imply a transient heat flow much less than its steady state value.
With regard to subsidence, integrating q v (in (34)) yields the total heat flow supplied by basal convection and hence the predicted deviation of subsidence from the square root of age relationship in (9) where t is a dummy variable for plate age. Equation (36) has been used in models of oceanic lithosphere subsidence, for example, equation (1) of Colin and Fleitout [67]. Their intent was to simply represent the increase in the vigor of convection as the rheological boundary layer thickens.

Comparison of oceanic heat flow with continental basal heat flow
The average heat flow added to the base of the plate in (36) is half the maximum at age t L . This effect needs to be included when one constrains the basal heat flow on older oceanic crust from observed deviation (or lack there of) of subsidence from the square root of age relationship in (8).
For an example, I assume that the heat flow through old oceanic lithosphere is in balance with stagnant lid convection of ~40 mW m 2 at 144 Ma. Then the average basal heat flow of ~20 mW m 2 would cause a 1000 m deviation for 144 Ma crust. A ~350 m derivation would not be evident in regions that appear to follow the square root of age relation for depth (8). Thus, a current basal heat flow in (34) (34) for 144 Ma oceanic crust is then 12 mW m 2 . This is near the limit of 14 mW m 2 that could be detected. It is thus acceptable that the steady state heat flow is the same beneath platforms and oceanic lithosphere. With regard to my assumed continental heat flow, Lévy et al. [68] (paragraph 67) and Perry et al. [69] use heat flow and xenolith systematics to constrain the Canadian Shield Moho heat flow to be between 12-18 mW m 2 . They note that xenoliths studies tend to yield higher heat flows than heat flow studies.

Depth-dependent viscosity
To illustrate further points, I assume for this subsection that seafloor depths actually trend to an asymptote different from the square root of age relationship because of strong stagnant-lid convection. I provide a simple numerical example. The heat flow is near its steady state value 41.7 W m 2 for 144 Ma crust and the seafloor is ~1000 m shallower than in the half-space model. The conductive geotherm with this heat flow intersects the MORB abiabat at 101 km. The platform geotherm in Figure 1 extrapolates to the MORB adiabat at 188 km. I proceed for illustration, assuming that stagnant-lid heat flow has reached steady state in both places.
In this regard, Doin et al. [70] proposed that viscosity increases with depth along an adiabat so that the different steady states exist beneath thin oceanic lithosphere and thick platform lithosphere. From the steady state linear formula in (24), the viscosity needs to be a factor of (41.7/ 17.2) 3 = 14.25 higher at 188 km depth than 101 km depth. Solving (11), D  is 32.7 km.
This value needs to be refined somewhat. Dumoulin et al. [71] correctly pointed out that strongly depth-dependent viscosity suppresses the upward increase of viscosity along a conductive geotherm. They recommended that T  be replaced with an effective temperature scale to account for this effect, where Z L is the lithosphere scale thickness. Equation (37) becomes singular when viscosity is constant on the conductive geotherm, unphysically all the way to the surface. Stagnant-lid formalism is not applicable as this limit is approached and a thermodynamic representation of viscosity is required to adequately impose shallow cold lithosphere. It is convenient here use steady state basal heat flow when comparing continents and oceans to obtain The second term in the bracket involves the thermal gradient q SL /k that is directly measured by xenolith studies. The predicted steady state heat flow is obtained by replacing T  with T eff in (24) which increases predicted heat flow by a factor of (T eff /T  ) 4/3 . I present a simple example assuming T  is 60 K and that the steady state ocean and platform heat flows are 17.2 and 41.7 mW m 2 . By iterating, a depth scale D  of 29.1 km provides this difference. The effective temperature scale in (38) is 70.4 and 93.7 K for the model ocean and continent, respectively For the assumed example value of T  = 60 K, equation (38) is away from its singularity and stagnant-lid formalism provides reasonable predictions. However, a larger value of the scale temperature will bring the equation toward singularity. For example, the continent parameters yield T eff = 214 K for T  = 100 K. In addition, Sleep [19] noted that the steady state thermal gradient at the middle of the rheological boundary layer is ~1/2 the lithospheric conductive gradient and recommended (37,38) with a factor of 2 in the numerator of the second term in the bracket. The actual factor in (37,38) for evaluating heat flow in (24) should be somewhere between 1 and 2 which brings the formalism closer to singularity.
With regard to ocean lithosphere, it is straightforward to use the transient heat flow in (18) as a basal boundary condition in plate models in analogy with the work of Dumoulin et al. [71]. Then (37,38) can be used to modify (18) or a numerical method can be used to evaluate the thickness of the rheological boundary layer for the assumed viscosity function. The transient heat flow in (18) depends on the inverse of viscosity not the inverse of the cube root for steady state. Therefore if viscosity is strongly depth dependent, the initial convective heat flow beneath thin lithosphere would rapidly increase. Thereafter as the lithosphere thickens, the effect of a thicker rheological boundary layer increasing heat flow is offset by the effect of higher viscosity decreasing heat flow. Such models could be calibrated so that basal heat flow remains nearly constant over a range of plate ages.
Returning to both oceans and continents, Doin et al. [70] note the possibility of two stability depths for stagnant-lid convection. Stability occurs when the predicted convective heat flow at the base of the lithosphere equals the predicted conductive heat flow ( Figure 10). The shallow equilibrium is stable. If the model system is started with thin lithosphere like at the ridge axis, the conductive heat flow exceeds the convective heat flow and the lithosphere thickens toward equilibrium. If the lithosphere is started too thick, the convective heat flow is greater than conductive heat flow and the lithosphere thins. However, if the lithosphere is started deeper than the second equilibrium point, conductive heat flow always exceeds convective heat flow and the lithosphere thickens indefinitely.
Doin et al. [70] point out that there may be a range of lithospheric thicknesses where the conductive and convective heat flows are close as in Figure 10. The lithospheric thicknesses slowly evolve and continental lithosphere may stay near the lower equilibrium for a considerable time. This inference is mathematically correct, but may not represent actual thick chemical lithosphere. That is, a range of cratonal lithospheric thicknesses greater than platform lithosphere thickness exists [61,62]. Such lithosphere is deeper than the postulated lower equilibrium. The presence of chemical lithosphere can only suppress heat flow. Thus the predicted cratonal convective heat flows are less than the conductive heat flow and craton lithosphere should have greatly thickened over its existence.
To reiterate, mildly depth-dependent viscosity may well exist near the base of the lithosphere. Its effect would mildly enhance the convective heat flow in ocean basins relative to that beneath continents. Strongly depth-dependent viscosity is required for multiple equilibrium lithosphere depths to exist. This hypothesis does not apply simply to cratonal lithosphere with a range of depths. However, it should not be discarded at present as discussed at the end of subsections of 5.3 and 5.4.

Nonlinear rheology
van Hunen et al. [54] and van Hunen and Zhong [55] proposed that nonlinear rheology can be detected by examining the apparent thermal age of oceanic lithosphere with a crustal age of ~90-130 Ma. They used seismic data from the northwestern Pacific Ocean to confirm that the deep lithosphere is in fact hotter and more buoyant than expected from the square root of age model. The bathymetry from the North Atlantic Ocean also indicates strong deviation from square root of age at ~90-130 Ma crust of the Bermuda, New England, Cape Verde, and Canary rises. This crust is shallower than the younger and older crust on its flanks (Figure 4). These shallow depths at ~90-130 Ma are evident in global and regional compilations where culling has not specifically removed them [27,29,31].
van Hunen et al. [54] correctly showed that nonlinearly rheology will cause convection to spin up rapidly as the lithosphere approaches its steady state thickness. The transient convective heat flow equation (18) illustrates the process. Convective heat flow depends on the thickness of the rheological boundary layer to the n+1 power. Until convection becomes vigorous, the rheological boundary layer thickness is proportional to lithosphere thickness. Hence convective heat flow is proportional to plate age to the (n+1)/2 power. van Hunen et al. [54] preferred n = 3.5 where convective heat flow increases with age to the 2.25 power. Heat flow is quite weak until the steady state lithosphere thickness is approached. There is a tendency for the initial vigorous convection to cause the rheological boundary layer to avalanche and hence thin the lithosphere below its steady state thickness, causing thermal uplift. van Hunen and Zhong [55] confirmed this behavior with three-dimensional convective models.
This hypothesis, however, does not provide a simple explanation for regions like the Sohm Plain [44] (Figure 4) that continue to follow the half space subsidence relationship. Moreover, both the New England and Bermuda rises appear to lie on hotspot tracks that crossed the continent and hence plume tracks. Physically, perturbations, including fracture zones, are present to trigger convection so the lithosphere cannot thicken greatly beyond its equilibrium thickness without nonlinear convection becoming quite vigorous. That is, the simple hypothesis implies shallow depths for all, rather than just most, 110 Ma crust.
However, this avalanche hypothesis may still be applicable as discussed at the end of subsection 5.5 if asthenospheric temperatures vary significantly laterally from a single adiabat. Depth-dependent viscosity may play a part.
Avalanche instabilities may also have relevance to continental lithosphere. The stress in (14) increases as the lithosphere and hence the rheological boundary layer thickens. The increased stress may cause the effective rheology to transition from basically linear to nonlinear. The net effect is a delamination avalanche where the temperature contrast across the rheological boundary layer is greatly and suddenly increased from (25b).
Physically, linear and nonlinear creep mechanisms occur in different places on a molecular to grain scale and hence act in parallel. The creep rate in (12) of coexisting linear and nonlinear mechanisms is formally, where the subscript lin indicates property of the linear mechanism and parameters without the subscript refer to a nonlinear parameter for this subsection only. The linear depth and temperature scales in (11) are not necessarily the nonlinear ones. These rheological modifications add several nearly free parameters, but under most conditions either the linear or the nonlinear creep will dominate as viscosity varies by orders of magnitude. It is thus inappropriate to extrapolate a power-law rheology for cold shallow lithosphere derived from plate flexure to the base of the lithosphere where linear creep may well dominate. Sleep [72] noted that stresses within the plate (as during continental collision, rifting, and conceivably plume impingement) may trigger this instability. There are no examples that have been related to this instability in mid-plate regions although it is a candidate for thinning the lithosphere beneath continental interior basins that later thermally subside [17,18]. Cooper and Conrad [73] modeled a related effect where the drag stress on the base of cratons increases as they thicken over time. Increased lithosphere thickness puts its base into more viscous mantle increasing the drag. Cooper and Conrad [73] contend that nonlinear rheology will limit the thickness of cratons at a time well into the geological future.
The effect of a delamination event on a xenolith geotherm would be apparent. The upper part of the geotherm would extrapolate to the base of the lithosphere prior to delamination and the lower part to the base after delamination. The curvature of the geotherm would be evident. Thinning of the lithosphere by ponded plume material would produce similar but moderate curvature in the geotherm [74]. The two cases can be distinguished in principle because the delamination event would not heat the mantle above the ambient MORB adiabat.

Dynamic topography and the fate of downwelling material
Lithospheric isostasy is clearly an approximation as lateral density variations clearly occur beneath the lithosphere and the uppermost asthenosphere, including slabs and plumes. The downwellings associated with small-scale convection produce local density heterogeneities. It is simplest to include density variations in the lithosphere, the rheological boundary layer, and any plume material ponded beneath the lithosphere in isostasy. Then dynamic topography can be defined conveniently as any process that causes the pressure at the compensation depth Z L in (4) and (5) to vary. The dynamic topography defined in this manner is weakly dependent on the choice of Z L if lateral density variations in the asthenosphere are in fact small.
It is inescapable that dynamic topography exists [75,76], but the point here is that it is unlikely to cause local variations in depth far from slabs like the difference between the Sohm Plain and the New England and Bermuda rises [44] ( Figure 4). Physically, distant dynamic topography does not depend on the details of its density source and reciprocally distant density variations cause only long wavelength variations in dynamic topography. For the mathematically inclined, the equations of motion (shown as (1a) and (1b)) result in an elliptical boundary value problem and the previous sentence is a form of Saint-Venant's principle.
Dynamic topography caused by dense material in the downwellings of convection is a related issue. If this dense material behaved according to simple isostasy, it would cause the seafloor to subside more, not less, than predicted by the half-space model. Two related processes are involved. First the viscosity of the mantle increases with depth below the asthenosphere. Deep density variations tend to produce downward directed forces that are ultimately compensated by relief on the core-mantle boundary not by seafloor topography. Positive geoid anomalies associated with slabs calibrate this effect [77].
Second and more importantly, most of the heat flow from the mantle arises from the internal radioactive heating and from internal transient cooling. The core-mantle boundary provides only a modest fraction of the total heat flow. This leads to an effect that Huang and Zhong [78] and Zhong et al. [31] call heat trapping.
Qualitatively, the thermal gradient is subadiabatic at mid-mantle depths [79]. Viscosity increases with depth so that surface plate moments do not effectively stir the deep mantle. Slabs sink and pond near the base of the mantle. Their cool dense material spreads out slowly. Subsequent slabs sink to the bottom displacing previous slab material upward. It takes ~2 Ga for slab material to return to the shallow mantle. A batch of dead-slab material remains essentially adiabatic except for radioactivity heating. During its ascent, the material heats up from radioactivity. The longterm cooling of the mantle acts as a virtual heat source in the batch relative to average mantle temperature.
Ridge axes are basically passive. The forced upwelling against the subadiabatic thermal gradient taps material at depth that is cooler than the material on the ridge flanks. The cool material is dragged laterally beneath the lithosphere above slightly warmer asthenosphere. Conversely, slabs drag down surrounding material causing it to be warmer than average for its depth (Figure 11). Buoyancy forces from this process contribute to dynamic topography and geoid anomalies.
Lateral variations of temperature rather than a single MORB adiabat are thus expected beneath plates. The smallscale convection off-axis thus mixes the cooler ridgeentrained material with the underlying warm material. Down-wellings below the warm region sink eventually into the underlying subadiabatic region. They tend to pond at a level of neutral temperature and buoyancy where they Figure 11 Highly schematic diagram illustrates temperature variation beneath plates. The mantle is stably stratified at middle mantle depths. Ridges drag potential isotherms upward and slabs drag them downward. The asthenosphere near the ridge axis is thus somewhat cooler that the asthenosphere in mid-plate regions. Downwelling from stagnant-lid convection encounter progressive cooler material as they sink and pond at a level of neutral buoyancy. Slabs sink into deep mantle composed of previous slabs. The temperature contrast of the slab with its surroundings decreases as it sinks. produce no further lateral density contrasts that would cause dynamic topography. The temperature difference between the downwellings and the ambient mantle is a fraction of the temperature scale T  (Figures 5 and 7).

Delayed avalanche instability
If the lateral asthenospheric temperature differences are a significant fraction of T  , the predicted transient heat flow (18) and steady state heat flow (23) will vary laterally. Cooler regions might have not yet reached the avalanche instability associated with nonlinear viscosity in subsection 5.3. This provides an attractive explanation for regions like the Sohm Plain ( Figure 4) that continue to subside with the square root of age [44]. I present of simple example with the North Atlantic in mind to estimate the temperature variation at the base of the lithosphere needed to cause the onset of vigorous nonlinear convection to vary. Convection is vigorous at ~110 Ma except in the Sohm Plain where it is not yet vigorous by ~170 Ma. That is, the value of t s needs to vary by a factor of 170/110 = 1.55. From (33) and (23), the scale time is proportion to the half-space viscosity to the 2/(n+1) power. For n = 3.5, the viscosity has to change by a factor of 4.26. The half-space temperature beneath the Sohm Plain thus needs to be 1.45T  less that where vigorous convection occurred.
Even the low value T  =40 K implies a significant temperature deficit of 58 K beneath the Sohm Plain. This is comparable to the temperature variation of ~50 K between the cool mantle of the Southeast Indian ridge compared with normal ridges. See the work Phillips and Coltice [80]. There is, however, no rapid plate divergence beneath the Sohm Plain to entrain dense mantle from depth and no obvious reason why local upwelling of cool dense material should occur in a mid-plate region. Neither is it likely that slab material underlies the Sohm Plain but not the adjacent seafloor. This simple form of lateral temperature variations is thus unattractive.
However conceivably, strongly depth-dependent viscosity by between 100-200 km depth might allow small variations in temperature to have a greater effect. Then vigorous stagnant-lid convection would occur beneath most old oceanic crust and weaker convection beneath platforms. Weak dependence of viscosity on depth between 200 and 300 km depth would allow chemical-lid convection to heat the base of all cratons. Such a formulation introduces lots of poorly constrained parameters. Insight on the rheology and global flow pattern of the mantle is needed for such modeling efforts to be productive rather than ad hoc.
Overall, three-dimensional thermal models that include plates and long-term thermal evolution of the Earth that allowed a subadiabatic gradient to form would quantify these issues. The numerical grid needs to be fine enough to resolve the rheological boundary layer and the ponding of downwellings in the mid-mantle below it. It would greatly help if plumes spontaneously formed in the model. It is also important to include continents. Even the sign of the temperature variation of continental asthenosphere relative to ridge asthenosphere is not well constrained. For reference, continental asthenosphere is hotter in the models of Phillips and Coltice [80].

Further work on three-dimensional dynamical features
So far in this review, three-dimensional relief on the base of the lithosphere acts a strong perturbation that assures finite-amplitude small-scale convection. Such relief slightly enhances the vigor of convection leading to a complicated transition between chemical-lid convection and stagnant-lid convection. In this section, I examine the potential use of three-dimensional features to constrain the physical properties of the mantle and the vigor of convection.
I confine discussion to oceanic lithosphere where the geometry is simply related to surface features. The lithospheric relief associated with fracture zones modulates the planform of convection [52,64,66]. Closely spaced hotspot lines of volcanic edifices may indicate underlying smallscale convection in somewhat hotter regions of the mantle [32,33]. The sedimentation rate with sedimentary basins may fluctuate with time [81]. I begin with ambient temperature variations in the asthenosphere that is relevant to these processes.

Ambient basal temperature variations
The temperature at the base of the lithosphere changes somewhat as the plate moves away from the ridge axis as discussed in Section 5.4. The numerical models in this paper do not provide evidence of the actual temperature change. I imposed a permeable basal boundary condition where the temperature of fluid entering from the bottom remained constant, so that comparison of model results with scaling relationships is simple.
The permeable boundary condition is not typically implemented into three-dimensional models for tractability. Still, local three-dimensional calculations do not provide much information on actual asthenospheric temperature variations. Rather modelers impose easily coded isothermal and no heat flow basal boundary conditions. The isothermal boundary condition eventually reheats cool material from downwellings thereby causing warm upwellings that interact with the stagnant-lid convection [78]. The vigor of convection wanes as the asthenosphere cools with a no heat flow boundary conditions (or isothermal ones where convection from below is sluggish) [82,83]. One can also generate heat internally within the model from radioactivity. The models of Dumoulin et al. [66] have significant internal heat generation.
Three-dimensional global models have "natural" mathematical boundary conditions at the surface and the coremantle boundary. The model Earth can be started hot following the moon-forming impact. The model details do not matter, as the mantle and the core then quickly cool to the slowly evolving state for the rest of the Earth's history. It is essential to include plates with slabs that emplace cold material near the base of the mantle. The numerical grid needs to be fine enough that stagnant-lid convection and plumes can arise from the dynamics. The models could conceivably impose plate boundary boundaries over the last 200 Ma where there are constraints on plate geometry and plausible plate motions at earlier times.
As a parallel approach, lateral variations in seismic velocities indicate lateral temperature variations in the asthenosphere. The inversion method needs to resolve lateral variations in the depth to the lithosphere-asthenosphere boundary from velocity variations within the asthenosphere. Small amounts of partial melt have significant effects so the variation of seismic velocity with temperature is not obvious; Till et al. [84] reviewed this topic. The asthenosphere with partial melt is likely to be layered and hence anisotropic [85].
Studies of the mantle source temperatures of mid-plate volcanoes calibrate the relationship between asthenospheric seismic velocity and temperature. Niu et al. [12] discussed the related issue of using this petrology to constrain lithospheric thickness. To have relevance to normal asthenosphere, the petrologist needs to concentrate on edifices formed well away from plume orifices were asthenospheric temperature varies slowly with position. At present it is difficult to recognize such features. For example, the number and location of plumes in the South Atlantic ( Figure 3) is poorly constrained. I return to this topic in Section 6.3 with closely spaced hotspot tracks.

Fracture-zone modulated convection
Study of the mantle beneath fracture zones has the potential to calibrate the vigor of small-scale convection. In the presence of conduction only, the deep boundary between thick old lithosphere and thin young lithosphere directly underlies the surface trace of ridge-ridge fracture zones. The escarpment on the base of the lithosphere modulates small-scale convection so that downwelling occurs on the old side and the boundary between thick and thin lithosphere migrates toward the old side as the plate ages [52,64,66]. Lithospheric relief associated with the base of the fracture zone disappears as steady state is approached [66]. The process is fully three-dimensional with flow parallel to the fracture zone [66].
The position of the lithospheric escarpment and its demise can in principle be detected seismically. Gravity data provide some constraint. Kawakatsu et al. [85] resolved the lithosphere-asthenosphere boundary beneath the Philippine Sea and the northeast Japan subduction zone with receiver functions. This method resolves waves converted at a sharp transition of seismic velocities, as do methods that resolve reflected waves. It thus reveals lateral variations in the depth to the convertor.
It is essential to deduce the nature of the convertor. The Philippine Sea interface lies beneath moderate-age oceanic crust at the depth where partial melt is expected to freeze either in a plate or half-space model. The Japanese Pacific plate interface is too shallow for a half-space model for ~130 Ma lithosphere [85]. Surface wave studies resolve this situation in the Pacific to the east [31]. The Pacific plate in this area has significant volcanic edifices so plumes rather than vigorous stagnant-lid convection may be involved.
The amplitude of the converted waves indicates that the wave convertor separates deep warm material with horizontal layers with partial melt from cooler shallower unstratified material [85]. The convertor persists within the subducting plate and has approximate depth range expected for isotherm. However, a thermodynamic boundary that quickly responds to changing pressure conditions does not behave in this manner. A purely melt convertor would disappear with tens of kilometers of subduction depth contrary to observations, as it freezes due to increasing pressure Further studies with closely based seismometers are needed to see whether the convertor responds rapidly enough to resolve the instantaneous lithosphere-asthenosphere boundary. Fracture zones are especially promising in that the conductive (and hence young plate) position of the lithosphere-asthenosphere escarpment is known.
The local sharpness of the convertor in principle provides information on the rheology of the mantle. Convection within a nonlinear material n  1 involves avalanches of the rheological boundary layer. The interface is hence sharp where an avalanche has recently occurred and diffuse elsewhere. Kawakatsu et al. [85] resolved similar interfaces beneath moderate-age lithosphere where avalanches are unlikely and ~130 Ma crust where avalanches are an attractive mechanism to explain the deviation of seafloor depth and asthenosphere depth from the predictions of a halfspace model. These inferences provide a very tentative argument against strongly nonlinear rheology and avalanches. This method would have utility if paleo-depth studies delineated candidates from recent uplift by avalanche.
Plumes provide a natural probe. Mantle plume material at the base of the lithosphere will pond on the young side and cascade from old side to the young side causing pressure release-melting and hence volcanic edifices. Sleep [52] and Dumoulin et al. [64,66] note that volcanic edifices do occur on the old sides of fracture zones where they would not be expected from conduction alone. Dumoulin et al. [66] in their review of available data do not obtain constraints that tell how close convection has approached its steady state value and hence constraints on mantle rheology.

Parallel hotspot tracks
Ballmer et al. [32,33] note that several systems of closely spaced hotspot tracks exist within the ocean basins. One must distinguish flowline edifices from those that align with the direction of absolute plate motion. The latter features have been associated with extension of the lithospheric plate [86]. Mittelstaedt and Ito [87] provided a compilation of both types. I note that many volcanic edifices are not well dated at present. I follow Ballmer et al. [32,33] and Holmes et al. [88] attribute the plate motion edifices to a combination of smallscale convection and channelized flow of warm buoyant material back toward the ridge axis. Plume material that is somewhat hotter than normal mantle may lower the viscosity facilitating convection. I continue to obtain scaling relationships.
Hot plume material spreads rapidly along the slope of the base of the lithosphere [89]. The ponded material also convects rapidly, cooling itself and thinning the overlying lithosphere. These processes occur quickly, leaving a thinned layer of warm plume material that convects slowly over a long period of time. The conductive heat flow from in the deep lithosphere (a few tenths of the surface heat flow in (7)) balances the heat flow from the plume material when lithospheric thinning ceases.
Mathematically, a calibrated form of the steady state convection equation (23) for this heat flow then involves three unknown rheological parameters, the temperature scale T  , the viscosity of the warm ponded material  warm (which I distinguish from the viscosity of asthenosphere at the MORB adiabat  half ), and the power n. Studies of the petrology of lavas sourced by the upwellings constrain the temperature difference between the warm material and the MORB adiabat. Further, the lines of volcanic edifices do not resemble the expected effect of avalanches of the rheological boundary layer, so it is attractive to assume that n1. Analysis is still model dependent. One may make an assumption on the steady state heat flow at the MORB adiabat to apply (23) to this situation. There are then two known heat flows from (23) at a known temperature difference from petrology that can be solved for T  and  half .
The mass balance of volcanism provides further constraint. Pressure-release melting occurs in the upwellings and the flux of melt is proportional to vertical velocity until the warm material is consumed. This velocity scales with the convective heat flow in (17), but is a factor of a few greater. This occurs because downwelling cool dense material drives the flow. The temperature difference between the downwelling and the mantle adiabat T down scales with rheological temperature contrast T rheo . The width of downwelling scales with the thickness of the rheological boundary layer and hence also with the rheological temperature contrast. If the lateral spacing of downwellings does not vary, the fraction of the base of the lithosphere that is be-neath downwellings scales with T down . Vertical velocities scale inversely with this fraction. Thus where the second bracket is the inverse of the fraction of the surface occupied by downwellings and there is no adiabatic region in the middle of the circulation when the temperature contrast between the downwelling and the adiabat T down = T max . Figure 7 provides an eyeball calibration of ~300 K and that for a linear viscosity, T down T  .
I provide examples for 25 Ma lithosphere with a surface heat flow of 100 mW m 2 and a small-scale convective heat flow of 30 mW m 2 . I ignore depth-dependent viscosity for simplicity. Equation (40) is then readily evaluated as the volume specific heat is ~410 6 W m 3 K 1 . The convective heat flow implies an upwelling velocity of 7.1 and 28.4 km Ma 1 , for T  of 100 and 50 K, respectively. Thus even modest stagnant-lid convection causes significant upwelling.
The buoyancy of source regions that have lost melt that ascended to the surface needs to be included in a numerical calculation [32,33]. I use the relationship from these authors that the density decrease of residuum is  depl = 72.6 kg m 3 time the fraction of melt removed F for simple scaling examples. The buoyancy of the residuum becomes a major hindrance to convection when it approaches thermal buoyancy T  . That is, the criterion for thermal convection is For  = 3.410 3 kg m 3 and  = 310 5 K 1 , the thermal buoyancy is 5 and 10 kg m 3 for T  = 50 and 100 K, respectively. Depletion of 7% and 14% melt produces these density differences, respectively. The value of T  in the models of Ballmer et al. [32,33] is ~180 K. However, the downwelling volume includes entrained material that has partially melted but otherwise remained near the abiabat. The average temperature of the downwelling material is hence 2 This provides a criterion of melting not significantly impeding convection Thus, even 1% and 5% melt depletion will interfere with convection for T  = 50 and 100 K, respectively. This effect counteracts the tendency of material with low T  to convect rapidly in (40). The value of T  = 180 K assumed by Ballmer et al. [32,33] is large enough that 12% melting with depletion is needed to significantly impede convection. The calculations of Humler et al. [40] put these melt depletions in context. Melting at their normal ridge begins at 2 GPa pressure (~60 km depth) and produces 7 km of crust. The mean melt depletion is thus ~11.7% and the maximum melt depletion about twice that. Their 50-K warmer mantle begins melting at 2.5 GPa (~75 km depth) and produces 8.6 km of crust. The mean melt depletion is ~11.5%. Thus the Ballmer et al. [32,33] models are in a parameter range where melt depletion has only a modest effect on smallscale convection near the ridge axis where the lithosphere has not yet thinned to 80 km depth. Melt depletion has a negligible effect for small fractions of melting off axis. Calculations with strongly temperature-dependent viscosity and explicitly nonlinear rheology where avalanches can occur are warranted. Seismic studies to resolve the lithosphere-asthenosphere boundary near active volcanic edifices that are aligned with absolute plate motions would be valuable.

Time-dependent convection beneath sedimentary basins
Thermal subsidence causes sedimentary basins to subside in analogy with oceanic crust. Sediments sequences deposited near seafloor record apparent subsidence. Petersen et al. [81] presented models to demonstrate that vagaries in smallscale convection cause the buoyancy of a column of material to fluctuate. With thermal isostasy, this process causes the rate of subsidence to vary from smooth curve. The effects of this process are most evident when basin subsides in the long term. That is, the transient in heat flow in (32) has approached the steady state value so convection is vigorous and the surface heat flow is substantially greater than this value.
Petersen et al. [81] used a very complicated rheology where the temperature contrast T rheo across the thermal boundary layer is readily not predicted from theory. I thus proceed in a manner that one would from seismic and sedimentological data where it is unnecessary to know the rheology. First, the lateral temperature variations can be inferred to be ~200 K from the Petersen et al. [81] color contoured figures. This quantity and the relief on the base of the lithosphere are potentially observable below a real basin from seismology. Second, the scale thickness of the rheological boundary layer is given by (19) in terms of heat flow. The amplitude of subsidence fluctuation (beneath air) is dimensionally where q v is the convective heat flow which can be inferred from lithosphere thickness. This quantity is measurable from the sedimentary record, adjusted for the isostatic effect for sediment and water loads. As pointed out by Petersen et al. [81], the thermal time scale of the rheological boundary layer is measurable from the time scale of sediment-rate fluctuations. It is given by Equations (42) and (43), contain only non-rheological physical properties, the convective heat flow q v and the rheological temperature contrast T rheo . The latter quantities are constrained to some extent independently of the sedimentary record. Petersen et al. [81] obtained variations of 2-20 Ma beneath relatively thin ~100 km lithosphere. Fully three-dimensional models of this process are relevant. Changes in the direction of plate motion should disrupt convection. Relief on the base of the chemical lithosphere may modulate convection cell and suppress transient variation. Nonlinear rheology should result in avalanches with sudden uplift in their aftermath.

Conclusions
The discovery of plate tectonics unified the earth sciences and provided a basis for quantitative geodynamics. Smallscale convection occurs at the lithosphere-asthenosphere boundary. This fact was appreciated early in the history of plate tectonics. Still small-scale convection is poorly understood with regard to its effects on observable geology and geophysics. Its existence is reasonably established beneath old stable continental lithosphere that transmits heat by conduction from the asthenosphere to the surface.
The effect of small-scale convection in ocean basins is difficult to separate from transient effects of mantle plumes. The depth of most oceanic crust approaches an asymptote rather than continuing to subside with the square root of age as expected from pure convection. I concentrated on the regions that did not approach an asymptote on the grounds that these hotspot-free regions were affected mainly by small-scale convection.
The resolution in depth-age data in these regions is a few hundred meters. The physical parameters that predict conductive depth-age for old oceans are imprecisely known, implying a similar model uncertainty for old oceanic crust. Acknowledging data and model imprecision, cooling by conduction of oceanic plates explains the first order relief of mid-oceanic ridges and the hypothesis that oceanic smallscale convection supplies an insignificant amount of oceanic heat flow is viable.
Modeling of the effects of small-scale convection is quite feasible. There are good scaling relationships to relate convective heat flow to physical properties and numerous numerical calculations. I treat the rheology of the mantle as unknown and limit the number of independent parameters by using (11) and (12). The temperature scale T  for viscosity is essential because stagnant-lid convection and plates are effects of strongly temperature dependent rheology. A linear viscosity that does not vary with depth along an adiabat to asthenospheric depths, perhaps surprisingly provides a reasonable representation of continental and oceanic observations.
More complicated models with depth-dependent and nonlinear rheology are permitted when the deviations of these models from simple models with constant asthenospheric viscosity and adiabat are subtle. Viable models that predict localized vigorous small-scale convection beneath ocean basins involve poorly resolved features of gross mantle convection including lateral variations of temperature in the asthenosphere. They require three-dimensional global and local calculations for appraisal. Insight into mantle rheology would help limit the number of free parameters in such efforts.
It would be quite helpful if the apparent subsidence from depth-age could be supplemented with actual subsidence. This requires understanding dynamic topography and hence gross convection. It may be possible to use boreholes into edifices to calibrate regional paleontological benthic fauna provinces and hence infer relative subsidence and uplift.
Xenolith geotherms (Figure 1) resolve the heat flow in continental mantle lithosphere. It is bothersome that the xenolith studies give significantly higher heat flow at the Moho than do heat flow studies. The linear geothermal gradient typically continues to the base of the sample suite as in Figure 1. The curved geotherm within the rheological boundary layer conversely is not observed. The data thus provide only an upper limit on the temperature contrast cross the rheological boundary layer. The depth and temperature of the base of buoyant chemical lithosphere at present provides better constrains on the rheological temperature contrast. Assuming a linear viscosity, continental observations indicate that T  <100 K and probably <60 K. The viscosity of the adiabatic asthenosphere then determines the laterally averaged heat flow through (24).
The rheological boundary layer between the rigid conducting lithosphere and the nearly adiabat asthenosphere has proved difficult to observe. The region is thin so it is not well resolved by tomographic seismic methods. In addition, small amounts of partial melt and their effects on rock fabric significantly affect seismic velocity so the relation between seismic velocity and temperature is not well constrained.
Currently, studies of converted seismic waves resolve lateral variations in the depth to the base of the oceanic lithosphere [85]. The convertor is not a simple phase boundary between melt-rich and melt-poor rocks as it advects downward with the slab. It is unclear until more measurements are taken whether the convertor adjusts rapidly enough that local variations in the base of the lithosphere associated with small-scale convection and fracture zones can be resolved. The strength of the converted phase gives some indication of the thermal gradient in the rheological boundary layer, but there are too few data to calibrate this effect.
Similar issues of relating seismic data to Earth processes exist with the lithosphere-asthenosphere boundary beneath continents. Combined converted-wave and surface-wave detect the boundary at 180-240 km depth beneath stable North America [90] in agreement with xenolith studies and heat flow systematics. The methods are sensitive to the detection of anisotropy associated with ancient processes above the boundary and current plate motions below the boundary. The change in anisotropy with depth is thus associated with rheology, but cannot be simply associated with the rheological boundary layer. The resolution is too coarse to look for downwelling material or even scalloped relief of the lithosphere-asthenosphere boundary.