PHYS 3070

Prof Louis Moresi

Australian National University

Continuum Mechanics and Mantle Dynamics

In which we derive the continuum equations for mantle / lithospheric flow, look at various questions of scaling, and develop some intuition for the dominant physics of mantle flow

Mantle Convection

Convection in Earth’s interior is (a little bit) like a boiling pot (as we saw in a previous lecture)

 

The hot soup rises to the surface, spreads and begins to cool, and then sinks back to the bottom of the pot where it is reheated and rises again. Why does hot soup rise and cold soup sink ?

General Observations on Convection

Without being particularly quantitative:

• Hot liquid is more buoyant than cold material and it tends to rise
• Cooler liquid is less buoyant and therefore tends to sink
• This can only happen if the two can move past each other
• Convection produces a self-stirring

Buoyancy forces are at work and viscous forces counteract these forces once the fluid is moving

\[ \textrm{buoyancy} \propto g \rho_0 \alpha(1-\Delta T) \]

Convection like this will only work when the soup is heated from below or, in the case of the Earth, if it is heated from within by radioactivity. (Can you see why) ?

General Observations on Convection

Definition: Convection is the transfer of heat by the self-organised movement of a fluid. Free convection is when the fluid is stirred entirely by rising buoyant material and sinking negatively-buoyant material. (Forced convection is produced when the fluid is stirred mechanically).

Convection is one of the ways we transfer heat from hot regions deep in the Earth to cooler, shallow regions.

Heat Transfer

Other ways that heat can be transferred include radiation, advection, and conduction.

Advection of heat is where some object carries an excess (or defecit of) heat energy from place to place.

Conduction is the transfer of heat (or electric current) from one substance to another by direct contact (lattice vibrations / electrons).

Radiation is the transfer of heat energy by photons that pass between two materials. It does not require any physical contact or indermediary material (e.g. it is not a problem to transfer heat across a vacuum this way).

“Kmecfiunit, cmglee, CC BY-SA 4.0 <https://creativecommons.org/licenses/by-sa/4.0>, via Wikimedia Commons”

Conduction & Fourier’s Law - i

Fourier’s Law tells us how much heat, \(Q\) is conducted through a sample in a unit time

\[ Q=-\frac{k A \Delta T}{L} \]

\(A\) is the cross sectional area of the sample, \(L\) is its length, \(k\) is a thermal conductivity “constant” dependent on the nature of the material and often also on the temperature, and \(\Delta T\) is the temperature difference across the sample. (Note the minus sign because heat will always flow from a higher temperature to a lower temperature.)

Conduction & Fourier’s Law - ii

Fourier’s Law is explained using a simple thought experiment such as this one but it actually refers to the fact that heat flows down a temperature gradient.

\[ q = - k \frac{dT}{dx} \]

\(q, (W/m^2)\) is called the heat flux and measures the flow of heat per unit area perpendicular to the temperature gradient

Thermal diffusion & Fourier’s Law

If the heat flux into a region is not balanced by the heat flux out, the region must change temperature. For example, if more heat flow into a region than is flowing out, it will become warmer.

The heat flow in the diagram is away (in both directions) from the high temperature to the cooler areas. The loss of heat from this region is balanced by a reduced temperature.

 

This process reduces the amplitude of temperature variations over time and makes the temperature smoother. Temperature extremes become smaller, gradients become lower, the heat flux everywhere becomes smaller: the temperature differences decay through time.

Newton’s Law of Cooling

Another example of “thermal decay” occurs when we have a hot body in a cool environment loosing heat through its surface according to the Fourier Law.

The heat flux is proportional to the difference in temperature between the hot object and the surroundings, so

\[ \frac{dT}{dt} \propto T - T_0 \quad \rightarrow \quad T = T_1 e^{-\lambda t} \]

This is exactly the same equation for the anomalous temperature as for the quantity of a radioactive isotope remaining after a certain time and has the same solution with \(\lambda\) a constant for a particular experiment.

Thermal diffusion

If a region has a different heat flux in v. out, then there must be a gradient in the heat flux. A change in temperature is driven by gradients in the heat flux like this:

\[ \frac{\partial T}{\partial t} = -\frac{1}{\rho C_p} \frac{\partial q}{\partial x} + H \]

\(H\) is internal heat generation which is often important in the Earth because of the decay of radioactive isotopes. Now we substitute for \(q\) from Fourier’s Law and have this expression:

\[ \frac{\partial T}{\partial t} = \frac{k}{\rho C_p} \frac{\partial^2 T}{\partial x^2} + H \]

Fourier was one of the first to solve this mathematically and developed his method of summing harmonics (series solutions) as part of this study.

Thermal Diffusion and Planetary Heat flow

The rate of heat loss is a primary controller of how a planet evolves.

• The heat loss controls the internal temperature of the planet and its mechanical properties are expected to be dependent on temperature

• The rate of heat loss from the interior is related to the energy available to supply mechanical work to deform the planet

• The rate of heat loss governs the amount of time the planet can remain “geologically” active before the interior stops exchanging with the surface / atmosphere.

Surface temperature and Heat loss

The Stefan-Boltzmann law tells us the rate at which a hot object loses heat to its surroundings.

\[ P = \varepsilon \sigma A \left( T^4_\textrm{hot} - T^4_\textrm{ambient} \right) \]

The average surface temperature of the Earth is about 15°C or 288K. This implies a loss of heat at a rate of 2 x 10¹⁷ W globally (to the chill of space) We can make some assumptions about the energy content of the Earth to see how long this can be sustained.

\[ \require{cancel} \require{color} \frac{dE}{dt} = \frac{dE}{dT} { \color{blue} \frac{dT}{dt} } = \varepsilon \sigma A \left( T^4_\textrm{hot} - \bcancel{T^4 _\textrm{ambient}} \right) \]

\[ \require{cancel} t_\textrm{cooling} = \frac{\rho C_p V_\oplus}{\varepsilon \sigma A_\oplus} \int^{T_\mathrm{final}}_{T_{initial}} \frac{1}{T^4}dT = \frac{\rho C_p R_\oplus}{3 \varepsilon \sigma} \left( \frac{1}{T^3_\mathrm{final}} - \bcancel{\frac{1}{T^3_\mathrm{initial}}} \right) \]

Online Calculator Tool (interactive version)

Surface temperature and Heat loss (2)

Assuming an initial temperature of 1500K (molten rock observed at the Earth’s surface) gives a cooling time of 250,000 years.

That’s quite young !

What is missing / incorrect in this analysis ?

• The incoming radiation from the sun (estimated 1.75 x 10¹⁷ W) implies the Earth’s surface is near equilibrium with the solar radiation and therefore we should not assume cooling to the cold depths of outer space (how long has this near-equilibrium been established ?)

• The assumption that the Earth’s interior is isothermal

• The assumption that heat loss from the interior to the surface is perfectly efficient

This analysis might be quite useful for a molten planet (magma ocean) where the heat transport in a vigorously convecting interior is very efficient. Once the planet begins to crystallise, conductive cooling through the chilled surface layers becomes the rate-controlling process.

Solution for the uniform conduction of heat out of a sphere

What is the equilibrium heat flow across a thin spherical shell ?

• Heat flux in balances heat flux out
• Account for heat generation within the shell

\[ Q_r = 4\pi r^2 \cdot q_r(r) \]

\[ Q_{r+\delta r} = 4\pi (r+\delta_r)^2 \cdot q_r(r+\delta_r) \]

The steady solution has this form

\[ T = -\frac{\rho H}{6 k} r^2 + \frac{C_1}{r} + C_2 \]

Which is only interesting if there is internal heating (\(H\)) and otherwise implies constant temperature / no additional heat at depth. Steady-state, of course, tells us nothing about the age of the Earth. \(C_1\) has to be zero, for the solution to be meaningful at \(r=0\).

Transient solution for a sphere cooling from a uniform T

The transient (i.e. interesting) equivalent to this problem is qualitatively different from the steady one. There are now variations with time and these trade-off with variations in the radial direction.

Fourier’s approach was to find solutions in the form of a series of terms, each of which individually satisfies the equations, but none of which, in isolation, looks much like the actual solution.

The equation we want to solve looks like this:

\[ \frac{\partial T}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r} \left( \frac{1}{r^2}\frac{\partial T}{\partial r}\right) \]

It says that the temperature changes due to radial temperature gradients and, equally, that radial temperature gradients are established by the changes in temperature.

NOTE: I have left out the material “constants” like thermal conductivity, radius of the sphere, density which means I am assuming these are all 1 (and constant)

Transient solution for a sphere cooling from a uniform T (2)

\[ \frac{\partial T}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r} \left( \frac{1}{r^2}\frac{\partial T}{\partial r}\right) \]

We can look for solutions that take the form …

\[ T(r,t) = \frac{2}{r} \sum_{n=1}^{\infty} \exp\left( n^2\pi^2 t \right) \sin(n\pi r) \cdot \int_0^1 {\color{red}{f_T(r^*)}} \cdot r^* \sin(n\pi r^*) dr^* \]

The term in red is the initial temperature distribution. If we assume that this is uniform everywhere, then this is simplified to the point where we can easily evaluate the solution:

\[ T(r,t) = \frac{2}{r} \sum_{n=1}^{\infty} \exp\left( n^2\pi^2 t \right) \sin(n\pi r) \cdot \left( \frac{\sin(n\pi) - n \pi \cos(n\pi)}{n^2 \pi^2} \right) \]

Transient solution for a sphere cooling from a uniform T (3)

\[ T(r,t) = \frac{2}{r} \sum_{n=1}^{\infty} \exp\left( n^2\pi^2 t \right) \sin(n\pi r) \cdot \left( \frac{\sin(n\pi) - n \pi \cos(n\pi)}{n^2 \pi^2} \right) \]

Ignoring how we come by this (complicated-looking) solution for the time being, what do these terms actually look like ?

A fundamental objection to Fourier’s original approach to solving equations by series solution is that the harmonics can violate the boundary conditions and other physical constraints.

Transient solution for a sphere cooling from a uniform T (4)

\[ T(r,t) = \frac{2}{r} \sum_{n=1}^{\infty} \exp\left( n^2\pi^2 t \right) \sin(n\pi r) \cdot \left( \frac{\sin(n\pi) - n \pi \cos(n\pi)}{n^2 \pi^2} \right) \]

How do the harmonics look when added together to give an approximation to the solution ?

The solution does satisfy the boundary conditions and the physical constraints. Note that the smaller time solution has much sharper gradients and many more terms are required to accurately approximate the solution.

Transient solution for a sphere cooling from a uniform T (4)

\[ T(r,t) = \frac{2}{r} \sum_{n=1}^{\infty} \exp\left( n^2\pi^2 t \right) \sin(n\pi r) \cdot \left( \frac{\sin(n\pi) - n \pi \cos(n\pi)}{n^2 \pi^2} \right) \]

Transient solution for a sphere cooling from a uniform T (5)

\[ T(r,t) = \frac{2}{r} \sum_{n=1}^{\infty} \exp\left( n^2\pi^2 t \right) \sin(n\pi r) \cdot \left( \frac{\sin(n\pi) - n \pi \cos(n\pi)}{n^2 \pi^2} \right) \]

The temperature gradient at the surface is found by differentiating this expression w.r.t. \(r\) and evaluating at \(r=1\)

\[ \frac{\partial T}{\partial r}(1,t) = -\sum_{n=1}^{\infty} \frac{2 e^{-n^2\pi^2 t} \left( n\pi \cos(n\pi) -\sin(n\pi) \right)^2}{n^2 \pi^2} \]

The surface heat flux is then given by

\[ q_0 = -k\frac{\partial T}{\partial r}(1,t) \]

Transient solution for a sphere cooling from a uniform T (5)

\[ \frac{\partial T}{\partial r}(1,t) = -\sum_{n=1}^{\infty} \frac{2 e^{-n^2\pi^2 t} \left( n\pi \cos(n\pi) -\sin(n\pi) \right)^2}{n^2 \pi^2} \]

This is, essentially, the simplest version of the Kelvin estimate for the age of the Earth based on a uniform internal temperature.

The cooling of a semi-infinite domain

This is a problem not unlike the cooling sphere, but with no inherent length scale that comes from the boundary conditions. The equations are similar but we are working in a Cartesian geometry (and also ignoring all the constants !)

\[ \frac{\partial T}{\partial t} = \frac{\partial^2 T}{\partial z^2} \]

This equation simplifies to an ODE if we make the following substitution,

\[ T^* = 1-T; \quad \eta = \frac{z}{2\sqrt{t}} \quad \rightarrow -\eta\frac{d T^*}{d \eta} = \frac{1}{2}\frac{d^2 T^*}{d \eta^2} \]

Which has the following solution

\[ \newcommand{\erf}{\mathop{\rm erf}\nolimits} \newcommand{\erfc}{\mathop{\rm erfc}\nolimits} T^* = 1-\frac{2}{\sqrt{\pi}} \int_0^\eta e^{-\xi^2} d\xi = 1 - \erf(\eta) = \erfc\left(\frac{z}{2\sqrt{t}}\right) \]

The cooling of a semi-infinite domain (2)

\[ \newcommand{\erfc}{\mathop{\rm erfc}\nolimits} T = \erfc\left(\frac{z}{2\sqrt{t}}\right) \]

This is a very different solution strategy to the cooling sphere (and yet they should, in a mathematical sense, approach each other as the radius becomes large).

The solution looks very similar in character: a cooling front which propagates in from the surface. The curves are self-similar, they just differ by a scaling factor.

Continuum mechanics - A general approach

\[ \frac{d}{dt} \int_{\Omega} \rho \phi d\Omega = - \int_{\Gamma} \mathbf{F} \cdot \mathbf{d}\boldsymbol{\Gamma} + \int_{\Omega} \rho H d\Omega - \int _{\Gamma} \rho \phi \mathbf{u} \cdot \mathbf{d}\boldsymbol{\Gamma} \]

Rate of change of a quantity, \(\phi\), is the balance between fluid transport across the boundary, \(\Gamma\), of an arbitrary volume, \(\Omega\), the diffusive flux, and sources within the volume.

The volume is arbitrary, so the equation must be valid for any choice of \(\Omega\) only guaranteed if the integrand is zero.

\[ \color{Gray}{\int_{\Omega} \Biggl[} \color{Black}{\frac{d \rho \phi}{dt} + \nabla \cdot \rho \phi \mathbf{u} + \nabla \cdot \mathbf{F} -\rho H} \color{Gray}{\Biggr] d\Omega} \color{Black}{= 0} \]

 

How it works

The general equation is

\[ {\color{Black} \frac{d \rho \phi}{dt} } + {\color{OliveGreen} \nabla \cdot \rho \phi \mathbf{u} } - {\color{Blue} \nabla \cdot \mathbf{F}} - {\color{Maroon} \rho H } = 0 \]

The colours here are:

• Black: Rate of change of \(\phi\) with time
• Green: Transport due to fluid movement
• Blue: Fluxes occuring in the absence of fluid motion
• Red: Sources of \(\phi\)

or, equivalently (a matter of notation)

\[ {\color{DarkGreen} \frac{D \rho \phi}{Dt}} - {\color{DarkBlue} \nabla \cdot \mathbf{F}} - {\color{DarkRed} \rho H } = 0 \]

• Green: Rate of change of \(\phi\) viewed in the fluid reference frame

Examples - 1. Mass

\[ {\color{Black} \frac{d \rho \phi}{dt} } + {\color{OliveGreen} \nabla \cdot \rho \phi \mathbf{u} } - {\color{Blue} \nabla \cdot \mathbf{F}} - {\color{Maroon} \rho H } = 0 \]

Let \(\phi \leftarrow 1\)

• There are no source terms for mass in classical mechanics, \(H=0\)
• There is no flux into / out of the elementary volume if the fluid is stationary, \(\mathbf{F}=0\)

Hence

\[ {\frac{\partial \rho}{\partial t} + \nabla \cdot \rho \mathbf{v} = 0} \]

Examples - 2. Heat

\[ {\color{Black} \frac{d \rho \phi}{dt} } + {\color{OliveGreen} \nabla \cdot \rho \phi \mathbf{u} } - {\color{Blue} \nabla \cdot \mathbf{F}} - {\color{Maroon} \rho H } = 0 \]

Let \(\phi \leftarrow C_p T\)

• Source terms possible, including “pure source” such as internal radioactive decay.
• Flux across the boundary is heat-diffusion given by Fourier’s Law

\[ {\frac{\partial \rho C_p T}{\partial t} + \nabla \cdot \left( -k\nabla T + \rho C_p T \cdot \mathbf{v} \right) - \rho H = 0} \]

Now, using mass conservation, we can simplify this to

\[ \left( \frac{\partial T}{\partial t} + \mathbf{v} \cdot \nabla T \right)= \kappa \nabla^2 T + \frac{H}{C_p} \]

Note the appearance of the directional derivative operator, \(\mathbf{v} \cdot \nabla\).

Examples - 3. Momentum

This is a vector equation so the conservation equation starts like this

Let \(\phi \leftarrow \mathbf{v}\)

\[ \frac{d}{dt}\int_{\Omega} \rho \mathbf{v} d\Omega = - \int_{\Omega} \rho g \hat{\mathbf{z}} d\Omega + \int_{\Gamma} \boldsymbol{\sigma} \cdot \mathbf{d}\boldsymbol{\Gamma}- \int_{\Gamma} \rho \mathbf{v} (\mathbf{v} \cdot \mathbf{d}\boldsymbol{\Gamma}) \]

• Source terms are body forces (here gravitational, buoyancy forces)
• Flux across the boundary is due to tractions applies on the boundary from stresses in the fluid:  \(t_i = \sigma_{ij} n_j\) where \(\hat{\mathbf{n}} = \left\{ n_j \right\}\) is the normal vector at a point on the surface

\[ \frac{\partial}{\partial t}(\rho \mathbf{v}) +\rho g \hat{\mathbf{z}} - \nabla \cdot \sigma + \nabla \cdot (\rho \mathbf{v} \mathbf{v}) = 0 \quad \textrm{or} \quad \frac{\partial \rho v_i}{\partial t} + \rho g \delta_{i3} - \frac{\partial \sigma_{ij}}{\partial x_j} + \frac{\partial(\rho v_i v_j)}{\partial x_j} = 0 \]

Now, using mass conservation, we can simplify this to

\[ \rho \left( \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} \right) = \nabla \cdot \boldsymbol{\sigma} - g\rho\hat{\mathbf{z}} \]

Note the appearance of the directional derivative acting on velocity, \(\left( \mathbf{v} \cdot \nabla\right) \mathbf{v}\) which is quite complicated in the vector version.

Examples - 3. Momentum (2)

The operator in the momentum equation \(\left( \mathbf{v} \cdot \nabla\right) \mathbf{v}\) stretches the vector notation quite hard, and demands explanation.

Let us look at a similar operator, the directional derivative of a vector, \(\mathbf{u}\):

\[ \begin{split} (\mathbf{v} \cdot \nabla) \mathbf{u} = & \mathbf{i} \left( v_1 \frac{\partial u_1}{\partial x_1} + v_2 \frac{\partial u_1}{\partial x_2} + v_3 \frac{\partial u_1}{\partial x_3} \right) \\ & \mathbf{j} \left( v_1 \frac{\partial u_2}{\partial x_1} + v_2 \frac{\partial u_2}{\partial x_2} + v_3 \frac{\partial u_2}{\partial x_3} \right) \\ & \mathbf{k} \left( v_1 \frac{\partial u_3}{\partial x_1} + v_2 \frac{\partial u_3}{\partial x_2} + v_3 \frac{\partial u_3}{\partial x_3} \right) \\ \end{split} \]

Continuum Mechanics: Summary

\[ {\color{OliveGreen}{\frac{d \rho \phi}{dt} + \nabla \cdot \rho \phi \mathbf{u}}} - {\color{Blue} \nabla \cdot \mathbf{F}(\phi, \nabla \phi)} {\color{Maroon} -\rho H(\phi, \nabla \phi) } = 0 \quad \textrm{[ General Template ]} \]

\[ \color{OliveGreen}{\underbrace{ \Bigl[ \frac{\partial T}{\partial t} + \left( \mathbf{v} \cdot \nabla \right) T \Bigr]}_{\dot{\mathbf{f}}}} - \color{Blue}{ \nabla \cdot \underbrace{\Bigl[ \boldsymbol\kappa \nabla T \Bigr]}_{\mathbf{F}}} - \color{Maroon}{\underbrace{\Bigl[ H \Bigl] }_{\mathbf{f}}} {\color{Black} = 0} \quad \textrm{[ Thermal Advection / Diffusion ]} \]

\[ \color{OliveGreen}{ \underbrace{ \Bigl[ S_s \frac{\partial h}{\partial t} \Bigr] }_{\dot{\mathbf{f}}}} - \color{Blue}{ \nabla \cdot \underbrace{\Bigl[ \boldsymbol\kappa \nabla h - \boldsymbol{s}\Bigr]}_{\mathbf{F}}} - \color{Maroon}{\underbrace{\Bigl[ W \Bigl] }_{\mathbf{f}}} {\color{Black} = 0} \quad \textrm{[ Darcy Flow (Groundwater) ]} \]

\[ \color{Blue}{ -\nabla \cdot \underbrace{\Bigl[ \frac{\eta}{2}\left( \nabla\mathbf{u} + \nabla\mathbf{u}^T \right)- p \mathbf{I} \Bigr]}_{\mathbf{F}}} - \color{Maroon}{\underbrace{\Bigl[ \rho \mathbf{g} \Bigl] }_{\mathbf{f}}} {\color{Black} = 0} \quad \textrm{[ Stokes Equation ]} \]

\[ { \color{OliveGreen}{ \underbrace{ \Bigl[ \frac{\partial \mathbf{u}}{\partial t} + \left(\mathbf{u}\cdot\nabla\right) \mathbf{u} \Bigr] }_{\dot{\mathbf{f}}}}} - \color{Blue}{ \nabla \cdot \underbrace{\Bigl[ \frac{\eta}{2}\left( \nabla\mathbf{u} + \nabla\mathbf{u}^T \right)- p \mathbf{I} \Bigr]}_{\mathbf{F}}} - \color{Maroon}{\underbrace{\Bigl[ \rho \mathbf{g} \Bigl] }_{\mathbf{f}}} {\color{Black} = 0} \quad \textrm{[ Navier-Stokes Equation ]} \]

\[ \color{OliveGreen}{ \underbrace{ \Bigl[ \frac{\partial h}{\partial t} \Bigr] }_{\dot{\mathbf{f}}}} - \color{Blue}{ \nabla \cdot \underbrace{\Bigl[ \boldsymbol\kappa(|\nabla h|) \nabla h \Bigr]}_{\mathbf{F}}} - \color{Maroon}{\underbrace{\Bigl[ k A^m |\nabla h|^n + \dot{h}_{b} \Bigl] }_{\mathbf{f}}} {\color{Black} = 0} \quad \textrm{[ Landscape evolution ]} \]

Navier-Stokes Equation: Constitutive Law

\[ \rho \left( \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} \right) = \nabla \cdot \boldsymbol{\sigma} - g\rho\hat{\mathbf{z}} \]

but to close this system of equations, we need to know how \(\boldsymbol{\sigma}\) relates to the momentum or velocity. In a Newtonian, viscous fluid:

\[ \tau_{ij} = \eta \left[ \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}\right] \quad \textrm{and} \quad \sigma_{ij} = \eta \left[ \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}\right] - \frac{\partial p}{\partial x_i} \]

are the deviatoric and total stress tensors, respectively. This definition of the stress in terms of a material property (\(\eta\), the dynamic viscosity), allows us to write:

\[ \rho \left( \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} \right) = \nabla \cdot \left[ \eta \left( \nabla \mathbf{v} + \nabla \mathbf{v}^T \right)\right ] - \nabla p - g\rho {\mathbf{z}} \]

\[ \rho \left( \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} \right) = \eta \nabla^2 \mathbf{v} - \nabla p - g \rho {\mathbf{z}} \] The second form assumes that the viscosity is a constant and produces the familiar form of the Navier-Stokes equation. The first form is required if the viscosity varies.

Navier-Stokes Equation / Convection

\[ \rho \left( \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} \right) = \eta \nabla^2 \mathbf{v} - \nabla p - g \rho_0 \alpha \left( T - T_0 \right) {\mathbf{z}} \]

\[ \left( \frac{\partial T}{\partial t} + \mathbf{v} \cdot \nabla T \right)= \kappa \nabla^2 T + \frac{H}{C_p} \]

We combine Navier-Stokes equation driven by thermal buoyancy, and the expression for the time-evolution of the temperature field. The two equations are coupled through \(T\). The various physical properties and their units are listed below.

\[ \begin{array} {l|l|l} \hline \textrm{Viscosity} & \eta & Pa \cdot s \\ \hline \textrm{Temperature scale} & \Delta T = T-T_0 & K \\ \hline \textrm{Lengths} & x & m \\ \hline \textrm{Thermal expansivity} & \alpha & K^{-1} \\ \hline \textrm{Gravity} & g & m.s^{-2} \\ \hline \textrm{Density} & \rho & kg.m^{-3} \\ \hline \textrm{Thermal diffusivity} & \kappa & m^2.s^{-1} \\ \hline \end{array} \]

Non-dimensional Navier-Stokes

It should not matter if we work in SI units cgs units or Imperial units — the physics does not change and the system does not care about our choice. That means there should be a way to express the equations without having to define specific units at all (or some sort of natural units for the problem that make everything close to 1).^[*]

We can rescale all our numbers by typical values of length, time, temperature etc. We have to be consistent (remember dimensional analysis) so these choices work, for example:

\[ \begin{array} {lll} x \leftarrow d . x' & \partial / \partial x \leftarrow (1/d) . \partial / \partial x' & \nabla \leftarrow (1/d) \nabla' \\ t \leftarrow d^2 / \kappa . t' & \partial / \partial t \leftarrow (\kappa/d^2) . \partial / \partial t' & \mbox{} \\ T \leftarrow \Delta T . T' & \mbox{} & \mbox{} \\ \mathbf{v} \leftarrow (\kappa / d) \mathbf{v}' & \mbox{} & \mbox{} \\ p \leftarrow p_0 - (\eta \kappa / d^2 ) p' & \textrm{Note:} \; \nabla p_0 = - g \rho_0 & \mbox{} \\ \end{array} \]

[*] this is what people used to do — feet / inches for carpenters, furlongs for farmers / horses.

Non-dimensional Navier-Stokes (2)

\[ \frac{\rho_0 \kappa}{d^2} \frac{D}{Dt'} \left( \frac{\kappa}{d} \mathbf{v}' \right) = \frac{\eta}{d^2} \acute{\nabla}^2 \left( \frac{\kappa}{d} \mathbf{v}' \right) - \frac{\eta \kappa}{d^3} \acute{\nabla} p' + g \rho_0 \alpha \Delta T T' {\mathbf{z}} \]

Tidying up as much as we can:

\[ \frac{\rho_0 \kappa^2}{d^3} \frac{D\mathbf{v}'}{Dt'} = \frac{\eta \kappa}{d^3} \acute{\nabla}^2 \mathbf{v}' - \frac{\eta \kappa}{d^3} \acute{\nabla} p' + g \rho_0 \alpha \Delta T T' {\mathbf{z}} \]

Divide throughout by \(\eta \kappa / d^3\) gives

\[ \begin{aligned} \frac{\rho \kappa}{\eta} \frac{D\mathbf{v}'}{Dt'} & = \acute{\nabla}^2 \mathbf{v}' - \acute{\nabla} p' + \frac{g \rho_0 \alpha \Delta T d^3}{\kappa \eta} T' {\mathbf{z}} \nonumber \\ & \textrm{or} \\ {\color{DarkRed}\frac{1}{{\rm Pr}}} \frac{D\mathbf{v}'}{Dt'} & = \acute{\nabla}^2 \mathbf{v}' - \acute{\nabla} p' + {\color{DarkRed} {\rm Ra}} T' {\mathbf{z}} \end{aligned} \]

We have separated out the form of the equation from its scaling by going through this process. The primed variables can be re-inflated to their normal size using the scalings from the previous slide. The new parameters in the equation are dimensionless ratios (they have to be since all the primed variables are ratios).

Non-dimensional Navier-Stokes (3)

Rayleigh number - measures the intensity of the buoyancy forces and their ability to counter the effects of thermal diffusion and viscosity in damping motion driven by buoyancy

\[ \mathrm{Ra} = \frac{g \rho \alpha \Delta T d^3}{\eta \kappa} \]

Prandtl number - measures the relative importance of stress diffusion and thermal diffusion and is purely a material property (no reference to the geometry, length- or time-scales).

\[ \mathrm{Pr} = \frac{\eta}{\rho \kappa} = \frac{\nu}{\kappa} \]

Reynolds number - measures the relative importance of inertial to viscous forces. The Reynolds number is commonly used in fluid dynamics where it is an indicator of the transition from steady flow, unsteady laminar flow, through to turbulent flow.

\[ \mathrm{Re}=\frac{\rho V_0 d}{\eta} = \frac{V_0 d}{\nu} \]

The Reynolds number does not contain any reference to the thermal state but this means that there is no natural length/time scaling and so the definition can be subjective.

Non-dimensional Navier-Stokes (4)

Questions

• What is the Prandtl number for (1) Water, (2) A planet ?
• What is the Reynolds number for plate-tectonic motions ?
• Estimate the Rayleigh number for the Earth.

Discussion points:

• If the Prandtl number is very large, how long does it take for the Navier-Stokes equation to adjust to changes in buoyancy forces or boundary conditions ?

• What is the Reynolds number for mantle flow / plate tectonics ? How about for flow in the Earth’s outer core ?

• A similar scaling question: what (roughly) is the kinetic energy of the Pacific plate ?

Incompressible flow and the stream function

In all the derivations so far we assume two incompatible things:

• \(\rho\)is constant (i.e. \(\nabla\cdot \rho \mathbf v \equiv \rho_0 \nabla \cdot \mathbf{v}\))
• Flow is driven by buoyancy forces that result from changes in density

It is not unusual to make the approximation that fluids are incompressible (even for air if the Mach number is “small”) and we do not lose a lot of detail if we do this. We assume that the density variations are small so they are only important when they appear in a term that does not depend on the actual density but only the density gradient (e.g. buoyancy)

• This is called the Boussinesq approximation.

Pressure becomes a constraint parameter rather than a physical pressure (cf a Lagrange multiplier) and the constraint equation is: \[ \nabla \cdot \mathbf{v} = 0 \]  Constrained PDEs can be very hard to solve, but there are some tricks that allow us to directly incorporate the pressure constraint.

Incompressible flow and the stream function (2)

The stream-function is a scalar quantity which defines a 2D incompressible flow everywhere. It satisfies: \[ v_1 = -\frac{\partial \psi}{\partial x_2}; \quad v_2 = \frac{\partial \psi}{\partial x_1} \quad {\color{DarkOrange} \textrm{compare this to } \mathbf{v} = \nabla \times (\psi {\mathbf{k}}) \textrm{: the curl of an out-of-plane vector} } \]

and this automatically implies \[ \frac{\partial v_1}{\partial x_1} + \frac{\partial v_2}{\partial x_2} = 0 \quad \textrm{(i.e. incompressibility)} \] and \[ (\mathbf{v} \cdot \nabla) \psi = v_1 \frac{\partial \psi}{\partial x_1} + v_2 \frac{\partial \psi}{\partial x_2} = \frac{\partial \psi}{\partial x_2} \frac{\partial \psi}{\partial x_1} - \frac{\partial \psi}{\partial x_1} \frac{\partial \psi}{\partial x_2} = 0 \] i.e \(\psi\) does not change as a result of advection

Note: a more general, 3D version of this derivation is the Helmholtz decomposition of a vector field into irrotational and a solenoidal (divergence free) parts:

\[ \mathbf{u} = {\color{DarkGreen} \nabla \phi} + {\color{DarkRed} \nabla \times \mathbf{A}} \]

Incompressible flow and the stream function (3)

Vorticity is defined as \(\boldsymbol{\omega} = \nabla \times \mathbf{v}\) but in 2D, there is just one component of this equation and we treat it as a scalar (\(\omega \leftarrow \omega_3\))

\[ \omega = \frac{\partial v_2}{\partial x_1} - \frac{\partial v_1}{\partial x_2} \]

If we take the curl of the Navier-Stokes equation we find we can eliminate the pressure gradients as \(\nabla \times \nabla P \equiv 0\). For infinite Prandtl number, the L.H.S vanishes.

\[ \frac{1}{\rm Pr} \left( \frac{D \boldsymbol{\omega}}{D t} - (\boldsymbol{\omega} \cdot \nabla) \mathbf{v} \right) = \eta \nabla ^2 \boldsymbol{\omega} + {\rm Ra} \frac{\partial T}{\partial x_1} \qquad(1)\]

Now, \(\mathbf{v} = (-\partial \psi / \partial x_2, \partial \psi / \partial x_1, 0)\) implies that \(\nabla \times \mathbf{v} = \nabla ^2 \psi = \omega\)  and, if we substitute this into \[ \color{DarkRed} \nabla^4 \psi = -{\rm Ra} \frac{\partial T}{\partial x_1} \textrm{(the biharmonic equation)} \]

where \(\nabla^4\)is known as the biharmonic operator and is defined to be

\[ \nabla^4 \equiv \nabla^2 ( \nabla ^2) \equiv \left( \frac{\partial ^4}{\partial x_1^4} + \frac{\partial ^2}{\partial x_1^2} \frac{\partial ^2}{\partial x_2^2} + \frac{\partial ^4}{\partial x_2^4} \right) \]

Does the rotation of the Earth matter ?

The Navier-Stokes equation describes the momentum balance in a viscous fluid in response to buoyancy forces and pressure gradients. In a rotating system, we may also need to consider the influence of Coriolis, and centrifugal forces

\[ \rho \frac{\partial \mathbf{v}}{\partial t} + \rho\left( \mathbf{v}\cdot \nabla \right) \mathbf{v} + 2 \rho \boldsymbol{\Omega}\times \mathbf{v} + \rho \boldsymbol{\Omega} \times \left( \boldsymbol{\Omega} \times \mathbf{x}\right) = \nabla \eta \left( \nabla \mathbf{u} + \nabla\mathbf{u}^T \right) - \nabla p - g \rho_0 \alpha \left( T - T_0 \right) {\mathbf{z}} % \frac{1}{\mu_0}\nabla\times\mathbf{B}\times\mathbf{B} \]

\(\boldsymbol{\Omega}\) is the rotation rate vector which we can also write as \(\;\boldsymbol{\Omega} \equiv \Omega_0 \mathbf{z}\)

\[ \frac{\rho _0 \kappa}{\eta _0}\left[ \frac{\partial \mathbf{v'}}{\partial t'} + \rho\left( \mathbf{v'}\cdot \nabla' \right) \mathbf{v'} + 2 {\Omega_0} \mathbf{\hat{z}} \times \mathbf{v'} \right] - \nabla' \eta \left( \nabla' \mathbf{v'} + \nabla'\mathbf{v'}^T \right) + \nabla' p^{*^{'}} = \frac{g \rho_0 \alpha \Delta T d^3}{\kappa \eta} T' {\mathbf{z}} \]

Pressure-like terms are all gathered into \(p^*\) and include contributions from the centrifugal force term

Convection and the Mantle: A Toy Model

Imagine a tank of thick, viscous fluid (like syrup, for example). This turns out to be a reasonable model for the interior of a solid (rocky / icy) planet if we consider the evolution on a geological timescale. A cold / dense blob of material sinks and stirs a tank of fluid just like this computer model shows.

What is viscosity ?

   

Viscosity is a measure of the resistance of a fluid to deform under shear stress. It is commonly perceived as “thickness”, or resistance to flow. Viscosity describes a fluid’s internal resistance to flow and may be thought of as a measure of fluid friction. Water is runny, having a lower viscosity, while honey is “thick” having a higher viscosity.

The right image is from the University of Queensland Pitch Drop Experiment

What is viscosity (ii)?

Viscous deformation is an irreversible flow that occurs in response to an applied shear stress. The stress is found to depend on the strength of the shearing velocity gradient.

\[ \tau = \eta \frac{d v_1}{d x_2} \]

Think of this as the stress that resists the shear deformation, i.e. how hard it is to stir the fluid. This is much harder if the fluid is more viscous. Or think of it as how fast the fluid responds to a given force (e.g. gravity) so a viscous gravity current will spread more slowly if the viscosity is high.

Viscosity only opposes the formation of velocity gradients; not a driving force, only a resistance.

The Toy Model Examined

In this toy model, do hot blobs and cold blobs behave in the same way (other than moving up v. moving down ?) - is this what you expect from the NS equation ?

The Toy Model - Examined

Now let us think about how the physical properties of the fluid and the boundary conditions change how this blob behaves.

\[ \textrm{Activity} = \frac{\textrm{Whatever tends to enhance the flow}}{\textrm{Whatever tends to retard the flow}} \]

We will consider how the behaviour changes with various parameters. When we change a parameter does it tend to enhance flow or retard it ?

Buoyancy forces

\[ \mathrm{buoyancy} \propto g\rho_0\alpha\left(1-\Delta T\right) \]

Lower strength of buoyancy forces on the left, increasing towards the right.

Does buoyancy drive the flow harder when it is a stronger effect ?

Viscous drag

Lower viscous drag stresses on the left, increasing towards the right.

Does viscous drag drive the flow harder when it is a stronger effect ?

Thermal conductivity (or diffusivity)

\[ \frac{\partial T}{\partial t} = \color{blue}{\frac{k}{\rho C_p}} \frac{\partial^2 T}{\partial x^2} + H = \color{blue}{\kappa} \frac{\partial^2 T}{\partial x^2} + H \]

Slower diffusion rates (lower diffusivity) on the left, faster towards the right.

Does thermal diffusivity drive the flow harder when it is a stronger effect ?

The Size of the Box

Is the flow moving faster when the box is larger ?

The Size of the Box (Scaled to fit)

Is the flow moving faster when the box is larger ?

Convective Activity

Intuition from these simple toy models produces the same dominant effects as we obtained from the non-dimensionalisation.

\[ \mathrm{Activity} \propto \frac{\textrm{buoyancy} \times d}{\kappa \eta} \]

cf. Rayleigh number that we found earlier.

\[ \mathrm{Activity} \propto \frac{\color{blue}{g\rho_0\alpha\Delta T} d^3}{\kappa \eta} \]

The advantage of the detailed analysis / re-scaling, is that we obtain the true dimensionless ratios.

Thermal Convection

In very-viscous convection, only the Rayleigh number plays a role (\(\mathrm{Re}\rightarrow 0\) and \(\mathrm{Pr} \rightarrow \infty\))

The structure of the flow and the pattern of the temperature field is predictable once we know the Rayleigh number but it is helpful first to develop some feeling for the patterns by seeing them.

Blue is cool, red is warm. The bottom of the box is kept hot, the top cold and the difference in temperature is \(\Delta T\), the depth of the box is \(d\) and so on. These numbers are set up so that the Rayleigh number is the required value. In a lab experiment, the simplest parameter to change is \(\Delta T\) assuming that the material has a constant viscosity with \(T\) and does not freeze, melt, burn or boil.

Thermal Convection (ii)

Rayleigh number: \(Ra=100\)

Thermal Convection (iii)

Rayleigh number: \(Ra=1,000\)

Thermal Convection (iv)

Rayleigh number: \(Ra=10,000\)

Thermal Convection (v)

Rayleigh number: \(Ra=100,000\)

Thermal Convection (vi)

Rayleigh number: \(Ra=1,000,000\)

Thermal Convection (vii)

Rayleigh number: \(Ra=10,000,000\)

Flat Earth, Round Earth

Rayleigh number: \(Ra=1,000,000\)

The pattern is not changed much when we make the box a little wider

Flat Earth, Round Earth (ii)

Rayleigh number: \(Ra=1,000,000\)

Curvature does not qualitatively alter the picture

Flat Earth, Round Earth (iii)

Rayleigh number: \(Ra=1,000,000\)

Some sensitivity to initial conditions

Flat Earth, Round Earth, 3D

In the sphere, the patterns are harder to see by eye, but it turns out that we can draw general conclusions from the 2D cases and they still hold in 3D.

Steady State Convection

In many of the examples, we reach a steady state in which the system stops changing and the movies appear to be frozen.

However, this steady state actually represents a dynamic balance between:

• heat conduction into and out of the model at the bottom and top respectively
• heat advection by continual rotation in the fluid that drags hot material between the two boundaries

Again, the best thing is to watch the system evolve and get a feel for what happens.

Steady State Convection (ii)

Rayleigh number: \(Ra=100,000\)

Steady State Convection (iii)

Rayleigh number: \(Ra=1,000,000\)

Steady State Convection

Rayleigh number: \(Ra=10,000,000\)

Understanding Thermal Convection

A relatively simple pattern emerges from the competition between advection of heat by fluid flow and diffusion:

• Convection cells that become independent
  • Boundary layers where horizontal advection balances vertical diffusion
  • Rotation in the middle of the box without much deformation

Understanding Thermal Convection (ii)

Convection length and timescales are determined by the Rayleigh number.

• High \(\mathrm{Ra}\) narrower boundary layers, higher temperature gradients

  • High \(\mathrm{Ra}\) more likely to be time-dependent, unsteady

  • Low \(\mathrm{Ra}\), convection suddenly shuts off

We will see this picture again when we consider the patterns of heat flow in the Earth’s ocean floor.

Average temperature

The temperature (velocity, stress etc) profiles vary systematically with \(\mathrm{Ra}\).

Nusselt number

Another, independent, dimensionless quantity can be derived for the convecting system and it is called the Nusselt number.

\[ \mathrm{Nu} = \frac{Q d}{k \Delta T} \]

Which is the ratio of the heat transported in the convection model (\(Q\)) compared to what would be transported if only conduction were at work. It turns out that

\[ \mathrm{Nu} \propto \mathrm{Ra}^{1/3} \]

The results (right) are from an experiment using liquid He which could not be much more different that solid rock but the principals still hold and the agreement with theory is excellent

Boundary layer approximation

•  

  • Thin horizontal boundary layers where heat diffuses into the layer and is carried away by horizontal advection. These are connected by vertical boundary layers of similar dimension.

  • Gradients along boundary layers are typically on a scale of the layer depth (large so leads to small gradients). Gradients across the boundary layers are typically large because the boundary layers are thin

Boundary layer approximation (2)

For the vertical thermal boundary layer, assume that:

\[ \frac{\partial v_2}{\partial x_2} \sim \frac{V}{\delta} ; \quad \frac{\partial v_2}{\partial x_1} \sim \frac{V}{d} ; \quad v_1 = 0 ; \quad \frac{\partial T}{\partial x_1} \sim \frac{\Delta T}{\delta} ; \quad \]

Then the curl of the Stokes equation (or the biharmonic equation, substituting for the stream function),

\[ \frac{\partial^3 v_2}{\partial x_1^3} - \frac{\partial^3 v_1}{\partial x_2^3} + \frac{\partial^3 v_2}{\partial x_1 \partial x_2^2} - \frac{\partial^3 v_1}{\partial x_2 \partial x_1^2} = \frac{g\rho\alpha}{\eta} \frac{\partial T}{\partial x_1} \]

becomes

\[ \frac{V}{d^3} + \frac{V}{d\delta} \sim \frac{g\rho\alpha}{\eta} \frac{\Delta T}{\delta} \quad \quad {\color{DarkRed} \frac{Vd^3}{d\delta^2} \sim \frac{g\rho\alpha d^3}{\eta} \frac{\Delta T}{\delta}} \]

Boundary layer approximation (3)

For the energy equation, assume steady state:

\[ \kappa \nabla^2 T = \left( \mathbf{v}\cdot\nabla \right) T \]

For the horizontal thermal boundary layer, assume that:

\[ \frac{\partial T}{\partial x_1} \sim \frac{\Delta T}{d} ; \quad \frac{\partial T}{\partial x_2} \sim \frac{\Delta T}{\delta} ; \quad v_1 \sim V ; \quad v_2 \sim 0 \]

Then the curl of the Stokes equation (or the biharmonic equation, substituting for the stream function),

\[ \frac{\partial^3 v_2}{\partial x_1^3} - \frac{\partial^3 v_1}{\partial x_2^3} + \frac{\partial^3 v_2}{\partial x_1 \partial x_2^2} - \frac{\partial^3 v_1}{\partial x_2 \partial x_1^2} = \frac{g\rho\alpha}{\eta} \frac{\partial T}{\partial x_1} \]

becomes

\[ \kappa \nabla^2 T = \left( \mathbf{v}\cdot\nabla \right) T \quad \quad {\color{DarkRed} \frac{\kappa \Delta T}{\delta^2} \sim \frac{V \Delta T}{d}} \]

Boundary layer approximation (4)

Now we have \(V/d \sim \kappa / \delta^2\) and some rearrangement of the previous approximations gives us

\[ V \sim \frac{\kappa}{d} \textrm{Ra}^{2/3} \]

and

\[ \textrm{Nu}^{-1} \sim \frac{\delta}{d} \sim \textrm{Ra}^{-1/3} \]

This approach balances the transfer of momentum and thermal energy vertically into the horizontal boundary layer (diffusion) with advection of each quantity along the boundary layer in a steady state.

Critical Rayleigh number

It can be a little difficult to see from a small set of movies, but if we reduce the Rayleigh number enough, conduction completely wins over advection and no motion takes place at all.

This can be approached theoretically, and we find that there is a critical value of the Rayleigh number and a critical wavelength where this transition occurs.

\[ \mathrm{Ra}_c = \frac{27}{4} \pi^4 = 657.51 \quad \textrm{and} \quad k = \frac{\pi}{\sqrt{2}} = 2.22 \]

If \(\mathrm{Ra} \gg \mathrm{Ra}_c\) then convective motion is inevitable

Discussion: Ra for the Earth’s Mantle

:::

Use the following values:

• Gravity, \(g \approx 10 \mathrm{m/s}^2\) through the mantle
• Density, \(\rho_0 \approx 3000 \mathrm{kg/m}^3\)
• Thermal Expansivity, \(\alpha \approx 10^{-5} \mathrm{K}^{-1}\)
• Temperature drop, \(\Delta T \approx 2000 \mathrm{K}\)
• Thickness of the mantle, \(d \approx 3000 \mathrm{km}\)
• Thermal diffusivity, \(\kappa \approx 10^{-7} \mathrm{m}^2 / \mathrm{s}\)
• Viscosity, \(\eta \approx 10^{22} \mathrm{Pa.s}\)

\[ \mathrm{Ra}_{\textrm{Earth}} = ? \]

::::

Is it above the critical value ?

What does that tell us ?

::::

Summary

Convection is a heat engine (i.e. it converts heat energy into mechanical work)

Convection is a balance between heat transported by fluid motion and diffusion. The fluid self-organises to create large scale patterns “out of nowhere”

In tanks of viscous fluids like syrup or honey, convection depends on just the one free parameter which is a combination of fluid properties, geometry and boundary conditions — this is called the Rayleigh number (Ra).

If we know Ra, we can predict the heat flow and typical velocity of the system

If Ra is below a critical value, convection dies away but if it is more than (about) ten times this value then convection cannot be suppressed. The Earth’s mantle is very much super-critical and so it must be convecting.

None of these simple models actually produce plate tectonics but they do still tell us about the heat flow in the Earth.

End