4.2.4. Isoviscous rheology in cylindrical geometry¶

4.2.4.1. Defining a coordinate system¶

In any convection model, gravity defines the natural down direction and gives us our first most important scale: the depth \(z\) from the surface, or its complement, the height from the model base \(h=1-z\).

If the domain is allowed to curve around a certain locus, a cylindrical or annular geometry is obtained which is more appropriate for planetary mantles. While we retain \(h\) and \(z\) as terms relevant to any action within the domain, we must also introduce a concept of radial height \(r\), understood here to represent the distance from the planetary centre of gravity. The cylindrical domain, for us representing the mantle, is thus bounded by the inner radius \(r_{i}\) and the outer radius \(r_{o}\), defining an area of \(\pi(r_o^2 - r_i^2)\).

Our choice of radii implies a degree of curvature \(f\):

\[ f \equiv \frac{r_o}{r_i} \]

Where \(f=1\) is equivalent to an infinitely wide Cartesian box, \(f\to0\) represents a full disk, and the values \(\sim 0.5\) and \(\sim 0.9\) would be appropriate for the whole mantle and upper mantle respectively. The ratio of radii \(f\) is identical to the ratio of circumferences, so that \(f=0.5\) represents a system where the arc length of the base is half that of the surface. (Note that this would imply infinite planetary radii at \(f=1\) - hence the planar-like endmember \(f=1\) is not strictly reachable under an assumption of curvature, though arbitrarily high values can be set to reproduce that behaviour [Jarvis, 1993].)

If we further stipulate:

\[\begin{split} \Delta r = r_{o} - r_{i} = 1 \\ r_{o} \to 1 \quad as \quad f \to 0 \end{split}\]

Then:

\[\begin{split} r_i = \frac{f}{1 - f}, \quad r_o = \frac{1}{1 - f} \\ r(h) = r_i + h \end{split}\]

Which suggests non-dimensionalising the planetary radial height \(r\) as \(r^{*}\) (r-star) such that \({r^{*}}_o = 1\):

\[\begin{split} {r^{*}}_i = f, \quad {r^{*}}_o = 1 \\ \Delta r^{*} = 1 - f \\ r^{*}(h) = \frac{h + r_i}{r_0} = \frac{\Delta r^{*}}{r} = h(1-f) + f \end{split}\]

This leaves us with four different terms to describe radial position: \(h\), the dimensionless height from the mantle base; \(z\), its complement; \(r\), the radial scale such that \(r_o - r_i = 1\); and \(r^{*}\), the radial scale such that \(r_0 = 1\). Each of these scales will prove natural in some contexts and less so in others, and all find use in our analysis.

To complete our coordinate system, we require an angular coordinate as well: the angle \(\theta\) in radians anticlockwise from an arbitrary origin. Often we will only want to reproduce a small wedge of the annulus, as time-dependence and numerical workload scale exponentially with aspect ratio; we can define our wedge selection \(\Theta\) in radians:

\[ \theta: 0 \to \Theta, \quad \Theta \leq 2\pi \]

If the simulation is to represent periodic flow around the annulus, values of \(\Theta\) must satsify \(\pi / l\), where \(l\) is a positive integer representing the number of convection cell pairs it would take to populate the full annulus at equivalent curvature (i.e. the number of upwellings or downwellings).

The selection of \(\Theta\) also selects an aspect ratio for the domain, but only once we choose a depth at which to measure it. It is most convenient to use the arc length through the mid-depth, which then relates to \(\Theta\) and \(f\) via a new term: the radial distance from the planetary centre of gravity to the mantle mid-depth, \(r_m\):

\[\begin{split} r_m \equiv \frac{r_{i} + r_{o}}{2} = \frac{1 + f}{2 \left( 1 - f \right)} \\ A = r_m \Theta \end{split}\]

If radius has been non-dimensionalised as above, this is equivalent to stating that the aspect ratio is equal to the area of the domain - regardless of \(f\) - just as it is in the Cartesian case.

Such a scheme leaves us with two competing claims for a ‘natural’ denominator of the angular coordinate - \(\Theta\) and \(r_m\). While authors have sometimes preferred to keep \(\Theta\) and \(r_m\) constant and allow \(A\) to vary [Jarvis, 1994], we have for the most part chosen to fix \(A\) and \(r_m\) with \(\Theta\) as the free parameter, as in [Jarvis, 1993]: this simplifies comparisons with plane-layer simulations at the cost of producing planforms which would be unstable if scaled to the full annulus.

Some system forcings, like internal heat, scale with area. The proportion of the annulus lying below a particular height \(h\) - which we shall call \(D\) for ‘disk’ - is a function of the inner and mid-radii:

\[ D(h) = \frac{r^2 - {r_i}^2}{2 r_m}, \quad r = h + r_i \]

Having declared that the total area shall always equal the aspect ratio \(A\), the true area under any depth \(h\) is simply \(D \cdot A\).

From \(r_m\) we are also able to derive \(s\), a statement of the angular length in ‘real’ terms (i.e. in the same metric as \(r\)). It is convenient to non-dimensionalise \(s\) as \(s^{*} = s / A\), such that the dimensionless length through the mid-depth \({s^{*}}_m = 1\). We can then write \(s^{*}\) very simply as a function of \(r^{*}\), and the inner and outer lengths accordingly:

\[\begin{split} s^{*} = 2 \frac{r^{*}}{1+f} \\ {s^{*}}_i = 2 \frac{f}{1+f}, \quad {s^{*}}_o = 2 \frac{1}{1+f} \end{split}\]

The length \(s\) is, among other things, the factor by which an average measurement of some variable taken across a layer can be converted into a total value for that layer. It is vital to account for varying \(s\) whenever comparing between different layers in a given system, or between equivalent layers in systems of differing \(f\).

4.2.4.2. Basal heating¶

4.2.4.2.1. Conductive solution¶

Fig. 4.3 Summary of the scaling behaviours of isoviscous conduction for varying curvature parameter \(f\). We obtain a natural scaling for \(f\) versus \(T_{av}\) with an \(R^2\) of 0.998.¶

Fig. 4.4 The analytical scaling of conductive temperature with \(\ln{r^{*}}/\ln{f}\) holds empirically with extreme precision.¶

It is a requirement of a conductive steady state (\(Nu=1\)) that the thermal flux must be the same through every layer. In the planar case this results in a linear geotherm which, in a model with fixed and unitless boundary temperatures, results in a simple function of \(T = z\) where \(z\) is dimensionless depth from the top of the model. The average temperature is then trivially \(T_{av}=0.5\). (For any system in pure conduction the Nusselt number is by definition \(1\).)

In a cylindrical domain, however, the length of each layer \(s\) is a function of depth and curvature as we have shown; consequently, shallower layers may transmit the same flux with a smaller temperature drop:

\[ \phi_q \propto s \cdot \frac{dT}{dh} \]

To define the flux, we need the geothermal gradient. The conductive geotherm can be elegantly stated in terms of \(r^{*}\) Fig. 4.3 Fig. 4.4:

\[ T(h) = \frac{\ln{r^{*}}}{\ln{f}} \]

And so the geothermal gradient:

\[ \frac{dT}{dh} = \frac{f-1}{r^{*}\ln{f}} \]

And finally the flux itself can be written as:

\[\begin{split} \begin{align*} \phi_q &\propto \frac{s^{*}(f-1)}{r^{*}\ln{f}} \\ &= \frac{2(1-f)}{(f+1)\ln{f}} \end{align*} \end{split}\]

\[\begin{split} \begin{align*} &\to -1 &as \quad f \to 1 \\ &\to 0 &as \quad f \to 0 \end{align*} \end{split}\]

Or very succinctly in terms of the ‘true’ radius of the mid-depth:

\[ \phi_q \propto \frac{1}{r_m \ln{f}} \]

To facilitate comparison between systems of different curvature, we can then use the above to define a dimensionless planetary flux \({\phi_q}^{*}\) - which is really just another name for the Nusselt number \(Nu\):

\[\begin{split} \begin{align*} {\phi_q}^{*} &= \frac{ {\phi_q} }{ {\phi_q}_c } \\ &\equiv Nu \end{align*} \end{split}\]

Where the subscript \(c\), here and elsewhere, denotes a purely conductive endmember. Because \(Nu\) now inherits a dependency on \(f\), it is no longer equivalent to the dimensionless surface temperature gradient, and so it is important always to present and discuss it in its proper terms as a ratio of fluxes.

Just as the flux now scales with \(f\), so must the average mantle temperature. In the planar case, the average temperature of the system is always half the temperature drop. In the cylindrical case, however:

\[\begin{split} \begin{align*} T_{av} &= \dfrac{1}{2} \large{\root \huge{e} \of f} \\ &\equiv T_{c} \end{align*} \end{split}\]

The relationship is apparent in the numerical results Fig. 4.3.

4.2.4.2.2. Instability and convection¶

An implication of \(Nu\)’s dependency on curvature is that the upper and lower boundaries must no longer be symmetrical. This invalidates many of the assumptions that made the planar case amenable to analysis. The additional space at the top of the model now allows more room for downwellings relative to basal upwellings, tending to promote instability [Jarvis, 1991]; on the other hand, the curved geotherm and the increased surface for radiating heat would tend to permit a comparatively thicker upper boundary layer. The effect of these countervailing forcings on the fundamental scalings of \(Nu\), \(Ra\), \(Ra_{cr}\), and the all-important relation \(Nu \propto R^{\beta}\) is not obvious.

To begin to unpack the complexities of convection in the annulus, we can start with the assumption that - as in the planar case - the convective steady state will eventually result in a broad intracellular region of uniform temperature \(T_{cell}\). Assuming a unit temperature drop \(\Delta T = 1\), we can write:

\[\begin{split} \begin{align*} {\Delta T}_o &= T_{cell} \\ {\Delta T}_i &= 1 - T_{cell} \end{align*} \end{split}\]

Knowing that the inner and outer fluxes \({\phi_q}_i\) and \({\phi_q}_o\) must be equal at steady state, and that the outer boundary - due to its greater length - can sustain that flux with a gradient shallower by a factor of \(f\), we can deduce a relation between the outer and inner thermal gradients, and thence between \(T_{cell}\) and the inner and outer boundary layer thicknesses \({\Delta r}_i\) and \({\Delta r}_o\):

\[\begin{split} \begin{align*} f \frac{{\Delta T}_i}{{\Delta r}_i} &= \frac{{\Delta T}_o}{{\Delta r}_o} \\ \frac{{\Delta r}_i}{{\Delta r}_o} &= f \frac{1 - T_{cell}}{T_{cell}} \end{align*} \end{split}\]

For each of the two layers, we can prescribe a layer-specific Rayleigh number accordingly:

\[\begin{split} \begin{align*} Ra_o &\propto T_{cell} {{\Delta r}_o}^3 \\ Ra_i &\propto (1 - T_{cell}) {{\Delta r}_i}^3 \end{align*} \end{split}\]

Having maintained non-dimensionality throughout, it is simple relate these two boundary Rayleigh numbers to the bulk \(Ra\) value:

\[ Ra_{layer} = Ra \cdot {\Delta T}_{layer} \cdot {{\Delta r}_{layer}}^3 \]

At this point, however, we have exhausted the insight we can obtain without making further assumptions. If we provide that the inner and outer boundary thicknesses must be the same, as they are in the planar case, we can see that:

\[ T_{cell} = \frac{f}{f + 1} \quad \leftarrow {\Delta r}_i = {\Delta r}_o \]

This, however, would imply that the inner and outer Rayleigh numbers are divergent. If we instead choose to conserve \(Ra\), then: ([Jarvis, 1993])

\[ T_{cell} = \frac{1}{1 + f^{-3/4}} \quad \leftarrow Ra_i = Ra_o \]

Both possibilities converge on \(0.5\) when \(f\to1\) and \(0\) when \(f\to0\), as we would expect.

However it is estimated, it is clear that, as \(Ra\) increases and boundaries thin, more of the mantle will fall in the intracellular region and global temperatures as a whole will approach \(T_{cell}\). Conversely, if \(Ra\) slips below its critical value, the boundary layers will disapper and the entire domain will enter the conductive regime: \(T^{av} = T_{c}\). These two temperatures therefore make up respectively the lower and upper endmembers of global temperature:

\[\begin{split} \begin{align*} T_{av} &\approx T_{c}, \quad Ra < Ra_{cr} \\ &\to T_{cell}, \quad Ra \to \infty \end{align*} \end{split}\]

It makes intuitive sense that the effect of increasing \(Ra\) should be to decrease global temperatures, since that is exactly why convection is preferred wherever possible - though this intuition may not hold for all rheologies.

Of course, what we desire most of all is a cylindrical scaling for the mantle convection power law \(Nu \propto R^{\beta}\). Following [Jarvis, 1993] and mandating equality of inner and outer \(Ra_{layer}\), it is possible to construct a ‘geometric correction’ \(g(f)\) that functions as a coefficient of the beta scaling:

\[\begin{split} g(f) = \frac{Nu_{c}}{{T_{cell}}^{4/3}} \quad \leftarrow Ra_i = Ra_o \\ Nu = g(f) \cdot R^{\frac{1}{3}} \end{split}\]

Using this scaling, Jarvis was able to obtain a beta exponent of \(0.321 \pm 0.001\) across four values of \(f\) from \((1.0 - 0.1)\) [Jarvis, 1993].

4.2.4.3. Internal heating in the annulus¶

4.2.4.3.1. Conductive solution¶

Fig. 4.5 Summary of the scaling behaviours of isoviscous conduction under internal heating \(H\) for varying curvature parameter \(f\) (colours as in previous charts). While samples of varying heat have been plotted, they do not appear in these charts due to the intentional factoring out of \(H\), demonstrating that this parameter is a simple coefficient.¶

It was established previously that, for an internally heated system, the geotherm and geothermal gradient are represented by:

\[\begin{split} \begin{align*} {T(h)}_{c(internal)} &= \frac{H}{2} \left( 1 - h^2 \right) \\ {{T(h)}_{c(internal)}}^{'} &= H\cdot h \end{align*} \end{split}\]

This is intuitive because the source flux visible to each layer is proportional to the area below that layer, which goes linearly with height \(h\) in a planar domain.

In the annulus, though, the proportion of the domain beneath a given height \(h\) is instead represented by \(D\), as we have shown. If we further assume that \(H\) is non-dimensionalised so as to represent the total flux of the model (i.e. it equals \(1\) for all geometries), then the flux through each layer height \(h\) of the annulus must simply be:

\[ {\phi_q}(h) = -H \cdot D(h) \]

We show that this holds exactly Fig. 4.5.

As before, the geothermal gradient required to transmit this flux must account for the varying layer length \(s^{*}\) - a function of \(h\) and the \(f\) parameter. Thus:

\[ \frac{dT}{dh} \propto \frac{\phi_q}{s^{*}} = -\frac{HD}{s^{*}} \]

The integral with respect to \(h\) yields the geotherm:

What kind of thing is a planet?