Chapter 2 Session 2

2.1 Part 1: One-dimensional diffusion

Lecture slides for this session can are here.

Next, we will familiarize ourselves with the case in which with a variable (in this case temperature) that changes its value not only through time, but also with location. The 1-D heat diffusion equation is given by ∂T∂t=κ∂²T / ∂z² , where the partial derivative symbol ∂∂ indicates that in the equation TT has derivatives to more than one variable (in this case to time tt and depth zz). A differential equation such as this one is therefore called a partial differential equation (PDE), and, in contrast, the DEs that we saw before (with derivatives to only one variable) are usually referred to as ‘ordinary differential equations, ODEs).

In this exercise, we will solve the heat diffusion equation for a situation of a cooling oceanic lithosphere: a hot column of mantle material (with temperature T=Tm=1350^oT=Tm=1350^o C) is emplaced at a mid-ocean ridge, where it suddenly will be in contact with cold ocean water. This column will start cooling from above.

2.1) Sketch (on paper) this temperature distribution TT as a function of depth zz for time t=0t=0 . This temperature distribution will be a good initial condition for our modelling exercise.

2.2) Our modelling domain will be the column of mantle material extending to, say, 300 km depth. We already have a governing equation (the heat diffusion equation described above) and an initial condition. To obtain a unique solution, we need to provide boundary conditions at the boundaries (i.e. end points) of the one-dimensional domain, e.g. a fixed temperature, that apply at every time step at the ends of our domain. What would be good boundary conditions for the temperature TT at z=0z=0 and z=300z=300 km?

2.3) What does the forward-Euler finite difference form of the diffusion equation look like?

2.4) We will first work out a timestep calculation by hand. Suppose we have a system, in which the mantle column is described by just 5 grid points (i.e. i=[0,1,2,3,4]i=[0,1,2,3,4]). At the boundaries (i=0i=0 and i=4i=4) boundary conditions apply, and at the remaining, interior points, the above finite difference equation applies. Assume the first solution T0iTi0 to be the initial condition you obtained in 2.1, and furthermore assume κ=10−⁶m² s, Δt=5Δt=5 Myr, Δz=75Δz=75 km. * Add your 5-point initial solution to your sketch. * Write out (again on paper, not in Python) for each point ii the equation that calculates the temperature at the new time step T1iTi1, where ii refers to the spatial dimension and the 11 to the time step. * Add this first timestep solution to the schematic plot of the initial condition. * (If time permits), apply a second time step and again add the solution to the plot.

2.5) Now copy the the following Python template heat1d.py into your Python editor. Read through, and try to understand how the code is set up. Then, add the missing heat diffusion function oneDdiff. Run the code, and check how well the numerical solution corresponds to the analytical one:

# heat1D.py
#
# Subitop modelling course, Edinburgh 2017
# Jeroen van Hunen, March 2017
# purpose: calculates 1D heat diffusion
# method: explicit time integration

import numpy as np
import pylab as plt
from scipy.special import erfc

# Subfuctions: 
### your oneDdiff function goes here ###
def halfspacecooling (Tm, z, kappa, t):
    # Calculates Halfspace cooling solution:
    fout = np.array(Tm-Tm*erfc(z/(2*np.sqrt(kappa*t))))
    return fout

# Main code: 
# Initialisation:
# Time variables:
dt       = 0.15                # timestep in Myrs
tmax     = 100
nt       = int(tmax/dt)+1      # number of tsteps to reach tmax Myrs
secinmyr = 1e6*365*24*3600     # amount of seconds in 1 Myr
dt       = dt*secinmyr         # unit conversion to SI: time in sec
time     = np.zeros(nt)
t        = 0                   # set initial time to zero
nplot    = 1                   # plotting interval: plot every nplot timesteps

# Mesh setup:
h        = 3e5                 # height of box: 3x10^5 m = 300 km
dz       = 1e4                 # discretization step in meters
nz       = h/dz+1         
dz       = h/(nz-1)            # Adjust reqested dz to fit in equidistant grid space
z        = np.linspace(0,h,nz) # array for the finite difference mesh

# Heat diffusion variables:
kappa    = 1e-6                # thermal diffusivity
Tm       = 1350                # mantle temperature in degC
Ttop     = 0                   # surface T
Tlith    = 1200                # T at base of lithosphere
Told     = Tm*np.ones(nz)      # initial T=Tm everywhere ...
Told[0]  = Ttop;               # ... except at surface, where T=0

# Timestepping
for it in range(1,nt):
    #update time
    t=t+dt 
    time[it]=t

    # numerical solution
    Tnew = oneDdiff(Told, nz, dz, kappa, dt)
                   
    # analytical solution
    Tana = halfspacecooling (Tm, z, kappa, t)

    if (it%nplot==0):
        # plot solution:
        plt.clf()
        plt.figure(1)
        plt.plot (Tnew,-z)
        plt.plot (Tana,-z)
        plt.xlabel('T [^oC]')
        plt.ylabel('z [m]')
        tmyrs=round(t*10/secinmyr)/10
        plt.title(' T after '+str(tmyrs)+' Myrs')
        plt.pause(0.0005)

    # prepare for next time step:
    Told = Tnew

A model solution is available here.

2.6) (EXTRA) If time permits, continue with the following exercise. We will define the thermal lithosphere as mantle material cooler than 1200^oC, and want to calculate how thick is the lithosphere as a function of its age. Add to your code a function that calculates the depth for which T=1200^oT=1200^oC, store this information in a time array, and plot it after the time looping has finished. The easiest method is probably to first work out the z-coordinate of the shallowest grid point for which T>1200^oT>1200^oC. For a smoother (and more accurate) solution, linearly interpolate between the grid points on either side of the lithosphere lower boundary to find the exact depth for which T=1200^oT=1200^oC.

2.7) (EXTRA) If time permits, replace your Forward Euler method with a Backward Euler and Crank-Nicholson methods: – Start from the basics and write out the finite difference discretisation for this model; – Put all known (i.e. old timestep) quantities to the right, all unknown (new timestep) quantities to the left of the ==-sign. – To solve this, you’ll need a system of equations solved as a matrix-vector problem. This is quite beyond the learning objectives of this course, and no model answers are given, but some more information on this is given in session 4. – Throughout this course, you can continue to add the features discussed for the Forward Euler method to this Backward Euler and/or Crank-Nicholson model setup.

2.2 Part 2: Numerical stability

2.7) Test the time step stability criterion for the radioactive decay problem. Open your radioactive decay code again, increase the total model time to 100 Gyrs, and empirically test which timesteps are stable or instable for the FE, BE, and CN timestepping methods. Calculate the maximum stable timestep for FE using the lecture notes, and compare it to your emperically-found maximum stable timestep.

2.8) Run the heat diffusion code multiple times with larger and smaller time steps, and larger and smaller number of grid points. Do you encounter unstable or unphysical solutions? Also in this case, BE and CN would result in unconditionally stable solutions, but solving a spatially one-dimensional (or more generally more-dimensional) problems with BE or CN is more complex (it needs to be solved using a matrix-vector system) than the one without any spatial variation (such as our radioactive decay problem).

2.3 Part 3: Natural boundary conditions

A bottom-boundary for vertical temperature models is often poorly defined. For example, the temperature at the bottom of the lithosphere is not well known. More often a heat flux from the mantle into the lithosphere is taken to calculate continental geotherms, because this avoids assigning a certain temperature to a certain depth. In the next set of exercises we will replace the bottom temperature boundary condition T=TmT=Tm by a bottom heat flux boundary condition −kdTdz=qm−kdTdz=qm.

Earlier, we derived a forward Euler code for heat diffusion, and applied it to lithospheric cooling. We used a function to calculate dT/dt for all interior points of the mesh, and explicitly set dTdt=0dTdt=0 for the boundary points. To implement a bottom heat flux boundary condition, we need to calculate dTdtdTdt for the bottom boundary point. For the calculation of the temperature change in the end point of the mesh, temporarily extend your mesh and temperature array with the imaginary point (or ‘ghost point’).

2.9. Adapt your 1-D thermal diffusion code to implement this heat flux boundary condition. Take a model box height of 150 km, q_m=0 mW/m², and k=3.3k=3.3 W/m,K. Replace the fixed-T bottom boundary condition by an insulating one: heat will disappear through the top boundary, but is not replenished from below, and ultimately the whole model domain should cool down.

A model solution is available here.