4.2.3. Isoviscous rheology with internal and mixed heating¶
On Earth and other planets, the geothermal flux is not solely driven by basal heating from the core: there is a substantial, variable, and transient contribution from radioactive decay as well [Daly, 1980]. In mantle convection studies this is dubbed ‘internal heating’, and a system driven by both basal and internal heating is said to be ‘mixed-heating’.
The ratio of internal heat production to surface heat flux inside a planet is called the Urey ratio \(Ur\) [Urey, 1955]. Values less than one signify that a planet on the whole is cooling, while values greater than one represent global interior warming. To avoid projecting unreasonable global temperatures into Earth’s past, it was originally calculated that, for a beta exponent of \(1/3\), the Urey number must be greater than \(\sim0.7\) [Christensen, 1985]. This assumption was forced into revision when it transpired that the radiogenic inventory of bulk mantle materials is far too low to support such a high radiogenic fraction [Jochum et al., 1983]. Direct evidence from geoneutrinos now suggests that the Urey ratio must be around \(0.5\) [Gando et al., 2011], and further argues that heat production is non-uniform throughout the mantle [Huang et al., 2013]. Thermal histories that honour these facts have proved very difficult to construct [Dye, 2012, Jaupart et al., 2015, Korenaga, 2003, Korenaga, 2008, Korenaga, 2017, Mareschal et al., 2012].
Analytically, the addition of internal heating introduces many complications. Most significantly, the temperature scale is no longer an independent variable, but a natural one emerging from the dynamic balance of heat production and heat transport. This dependency invalidates the derivation of \(Ra\) previously provided, and unfortunately, no obvious or unambiguous alternative scaling is known [Korenaga, 2017, Moore, 2008, Schubert et al., 2001].
4.2.3.1. Internal heating¶
Fig. 4.1 Summary of the scaling behaviours of isoviscous conduction for varying internal heating parameter \(H\). Temperature goes linearly as a function of dimensionless height from the base of the mantle, \(h^2\).¶
To provide some bounds on the expected behaviour of a mixed-heating model, it is helpful to review the effects of purely internal heating. Without the ability to arbitrarily specify a temperature scale, it is natural to use that which arises from pure conduction. Under a scenario of pure internal heating, the conductive geotherm is no longer linear, but scales with length squared [Turcotte and Schubert, 2014] Fig. 4.1:
Where \(H\) is the per-mass heating rate and \(h\) is the dimensionless height from the base of the mantle. The maximum mantle temperature possible for a given value of \(H\) is thus \(H/2\). Because, at the outer boundary, the thermal flux is exactly \(H\) (or, technically, \(H \cdot \rho\) if dimensionalised), \(H\) here is identical to the ‘conductive Nusselt number’, \({Nu}_{c}\).
Having obtained the conductive geotherm, we can select a nominal temperature scale:
Where \(k\) is the conductivity, \(\rho\) is the density, \(H\) is the per-mass heating rate, and \(b\) is the characteristic spatial length. The temperature scale allows an ‘internally heated Rayleigh number’ \({Ra}_H\) to be defined [Roberts, 1967]:
Where \({Ra}_B\) is the Rayleigh number derived for basal heating, introducing the parameters of thermal expansivity \(\alpha\), gravity \(g\), thermal diffusivity \(\kappa\), and momentum diffusivity \(\nu\) [Turcotte and Schubert, 2014]. As in the basally-heated case, there exists for this number a critical value \({{Ra}_{H}}_{cr}\) which divides the purely conductive regime from the convecting regime [Roberts, 1967]:
Usually we would be interested in deriving some scaling for the Nusselt number with respect to \({Ra}_H\). In this case, however, \(Nu\) must always equal one. This follows because there are no other sources or sinks in the model except for \(H\), which must be radiated in full from the upper boundary whether the interior is convecting or not. Rather than serving to augment the surface flux, convection in a purely internally-heated system only functions to smooth out interior temperatures. To describe this, we must choose a represenative internal temperature \(T_{i}\), which should reflect a typical temperature in those regions where advection dominates [Solomatov and Moresi, 2000]. Normalising \(T_{i}\) by the true (measured) temperature scale \(\Delta T\) gives \({T_{i}}^{*}\), the non-dimensional internal temperature, which then informs the definition of an ‘internal Rayleigh number’ \({Ra}_i\):
This observes, in other words, that - in an internally-heated system - the thermal contrast available to drive convection will always be less by some factor than the thermal contrast across the system as a whole.
While the flux across the upper boundary is fixed by the choice of parameter \(H\), the actual geometry of the boundary that supplies this flux must be a function of the Rayleigh number. At steady state, the interior temperature and the layer thickness adjust until a balance is reached, resulting in an outer temperature drop \({\Delta T}_o\) and outer boundary thickness \({\delta}_o\) that scale with, but are not solely determined by, the heating rate [Schubert et al., 2001]:
Wherein we observe that the implied flux (temperature drop divided by boundary thickness) remains equal to \({Nu}_c = H\), and hence \(Nu=1\), as expected: all of that inner motion occurs without leaving any thermal trace on the surface.
By taking a boundary layer approach to the heat-producing regions [Jaupart and Mareschal, 2010] and exploiting certain known requisites of the purely conductive state [Vilella and Kaminski, 2017], one can derive an alternative statement of the expected outer boundary layer properties in terms of the critical Rayleigh number:
This relationship appears to hold empirically with good confidence [Vilella and Kaminski, 2017].
Qualitatively, the addition of internal heating to the isoviscous convection problem imposes an asymmetry between upwellings and downwellings, which would otherwise be temperature-reversed mirrors of one another [Weinstein and Olson, 1990]. Because the thermal gradient of the conductive state rapidly shallows with depth, the local Rayleigh numbers of deeper layers can quickly become subcritical, such that a large portion of the domain from the base up is locked in a near-conductive state. At high \(H\), large-scale motions cease, and only a thin sub-surface ‘mixing layer’ witnesses any meaningful convection at all [Parmentier et al., 1994].
4.2.3.2. Mixed heating¶
Fig. 4.2 Summary of the scaling behaviours of the conductive solution for mixed heating. The inner boundary is set to dimensionless temperature \(T=1\), the outer to \(T=0\). On the right, average temperature and the slope of \(\delta T / \delta h\) are given as functions of \(H\). Gradient clearly varies as a function of dimensionless height \(h\) above the mantle base according to a slope given by \(-H\). The conductive geotherm for mixed heating is therefore, in fact, a parabola.¶
For systems driven by both internal and basal heating, the choice between basally-derived and internally-derived Rayleigh numbers leads to ambiguity and confusion. One path forward is to define a dimensionless heating rate \(H\) [Schubert et al., 2001]:
The conventional basally-heated \(Ra\) derivation can then be used, with \(H\) as a correcting coefficient.
While the purely internally-heated case has only one boundary layer (the outer), and the purely basally-heated case has symmetrical outer and inner boundaries, the mixed-heating case has two asymmetrical boundaries which are at least pseudo-independent of each other. While, at steady-state, the outer boundary flux is constrained to be always greater than the inner, the inner boundary flux may freely adjust itself as a function of the temperature differential between the prescribed lower boundary condition and the interior temperature. This intrinsic feedback makes quantitative analysis much more complicated.
As usual, a close consideration of the situation at the boundaries can set us on the right course. In the sorts of scenarios relevant to planetary settings, the lower boundary will remain hotter than the interior and the flux will be positive and substantial, in which case the outer boundary is required to transport the volumetric heating \(H\) in addition to the basal heating [Moore, 2008]:
Where \({q}_{o}\) and \({q}_{i}\) are the inner and outer heat fluxes.
At the limit where internal heating becomes negligible compared to the inner-outer temperature drop, internal heating can be ignored altogether and the system approximates a purely basally-heated model. At the opposite extreme, if the interior temperature due to internal heating is high enough, the inner flux may become negative, and the mantle will cool into the core: a case hard to envision in nature. Between these regimes is a transition point when the interior temperature becomes equal to the lower boundary condition. In this borderline case, the inner (core-mantle) flux is exactly zero (i.e. insulating), so that the mixed-heating system reproduces a purely internally-heated system. Thus, the mixed-heating model includes both internally-heated and basally-heated endmembers.
To sharpen this analysis, we can look at more closely at the conductive profile. While it is possible to derive this from first principles, it is more easily illustrated by a numerical approach Fig. 4.2. It is clear that the conductive geotherm forms a parabola whose second derivative is exactly equal to \(-H\); the rest follows by integration, with the constants provided by logic and observation:
If each geotherm traces a parabola, a maximum ‘natural’ temperature is implied where the first derivative is zero:
If \(h_{T_{max}}\) is less than zero for a particular value of \(H\), then that maximum will never be realised and the true maximum temperature will be that at the mantle base, which is unit in our dimensionless treatment. The condition \(h_{T_{max}} > 0\) thus represents the regime boundary between those conductive solutions that cool into the core and those that only cool into space. This occurs at exactly \(H = 2\), or more generally:
Where \(H_{crit}\) is the critical heating factor above which the mantle cools into the core and at which the lower boundary becomes effectively insulating, as previously discussed.
If we define a ‘conductive Nusselt number’ - something of an oxymoron - we can see how the true Nusselt number should be expected to scale:
I.e. the surface flux in the mixed case is simply the sum of the two heat drivers considered separately, just as we would expect from first principles.
We now have an understanding of how the system behaves below the critical Rayleigh number; now we must go above. One thing we can immediately observe from first principles is that the outer heat flux should be equal to the balance of \(H\) and whatever flux, positive or negative, is occurring over the lower boundary:
To analyse the mixed heating case, then, we need to invert our perspective and focus not on the outer boundary but on the inner boundary, whose core-ward temperature is equal to \(1\) and whose mantle-ward temperature - the ‘interior temperature’ \({T_i}^{*}\) - is to be determined. Unfortunately, no universally accepted analytical treatment of \({T_i}^{*}\) post-convection has yet been devised.
One workaround starts with the interior temperature analytically obtained for a purely internally-heated system and modify it with empirical constants to accommodate the addition of basal heating [Moore, 2008]:
This scaling, however, is back-formed from an observed and arbitrary measurement of \(T_i\) taken across the middle 60% of a finite-element numerical model, and no justification for this particular choice is ventured nor defended across the range of \(Ra\) and \(H\) values sampled.
Given a known value of \(T_i\), the associated Nusselt number should be proportional to that temperature drop divided by the stable outer boundary thickness \({\delta}_o\). Conventional boundary layer can be used to argue a scaling of the form [Schubert et al., 2001]:
Empirically, the same authors [Moore, 2008] found that this scaling holds only for a much lower exponent of \(\sim -0.303\) for \(Ra<=10^8\) and \(H<10\), while apparently approaching the theoretical value with increasing \(H\); the discrepancy was attributed to disruption by plumes from the lower boundary.
The next step is to attempt to derive the Nusselt-Rayleigh scaling, which we expect to approach the conventional power law behaviour under the beta exponent \(\beta = 1/3\). The expected relationship, however, does not manifest even approximately for the mixed heating case unless alternative definitions of \(Nu\) and \(Ra\) are used. Subtracting the basal component from \(Nu\) and using the arithmetic difference between \(Ra\) and \({Ra}_{cr}\) yields the following with some confidence [Moore, 2008]:
While the above was calculated using a planar 2D model, an equivalent model suite in the spherical shell closely reproduced this relationship after geometry was accounted for [Weller et al., 2016].
An alternative parameterisation [Vilella and Deschamps, 2018] expresses \(Nu\) in terms of the critical heating rate \(H_{crit}\) after convection:
Where \(C_N\) is an empirically-obtained constant argued to be exactly equal to \(1.5\). This definition of the convective \(H_{crit}\) further defines two independent Nusselt numbers depending on the magnitude of heating:
Where \({\Delta T}_{TBL,b}\) is the temperature jump across the inner boundary layer, recalling that the mantle must cool into the core above \(H_{crit}\).