Skip to main content

session3

Session3, Part 1: 2-D diffusion

Lecture slides for this session can are here.

Many geodynamical problems cannot be addressed satisfactorily in 1-D, and require (at least) a 2-D solution method. Today, we’ll set up a code for that. We will apply the 2-D solution method discussed in the lecture on a 1000×1000 km vertical cross section of the top of the mantle, in which around 660 km depth (the interface between upper and lower mantle) a 600 km long (i.e. wide) piece of subducted slab is positioned. A numerical heat diffusion code can be used to examine how this cold slab will thermally equilibrate with its hot surrounding.

3.1. Copy the following Python template into into your Python editor.

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

import numpy as np
import pylab as plt

# Subfuctions: 
### twoDdiff function goes here ###

# Main code: 
# Initialisation:
# Heat diffusion variables:
kappa    = 1e-6                # thermal diffusivity

# Mesh setup:
h        = 1000e3              # height of box: 1000 km
w        = 1.0*h               # box of aspect ratio 1
dx       = 1e4                 # discretization step in meters
dz       = 1e4
nx       = w/dx+1
nz       = h/dz+1         
dx       = w/(nx-1)            # Adjust requested dx & dz to fit in equidistant grid space
dz       = h/(nz-1) 
x        = np.linspace(0,w,nx) # array for the finite difference mesh
z        = np.linspace(0,h,nz)
[xx,zz]  = np.meshgrid(x,z)

# Time variables:
nt       = 50                  # number of tsteps
secinmyr = 1e6*365*24*3600     # amount of seconds in 1 Myr
t        = 0                   # set initial time to zero
nplot    = 5                   # plotting interval: plot every nplot timesteps

# Initial condition:
Tm       = 1350                # mantle temperature in degC
Ttop     = 0                   # surface T
Told=Tm*np.ones(np.shape(xx))  # Create T-field with T=Tm everywhere

iplot=0
# timestepping
for it in range(1,nt):
    # numerical solution
    Tnew = twoDdiff(Told, dx, dz, kappa, dt)
    
    #update time
    t=t+dt                  # in sec
    
    # plot solution:
    if (it%nplot==0):
        tmyrs=int(t/secinmyr)   # time in Myrs
        plt.clf()
        plt.figure(1)
        plt.imshow(Told, 
                   extent=(0,h,0,h),
                   clim=(0,Tm),
                   interpolation='bilinear', 
                   cmap='jet')
        plt.title('T after '+str(tmyrs)+' Myrs')
        plt.pause(0.00005)
  
    # prepare for next time step:
    Told = Tnew

3.2. Add/complete the missing parts:

  • Put a 100-km thick and 600-km wide slab with T=0oC as initial condition around 650 km depth within a mantle with uniform temperature of 1350oC.
  • An automatically calculated stable time step
  • A subfunction for the calculation of the new temperature as a function of the old temperature field for each of the nodal points, assuming essential boundary conditions.

3.3. Once the code is working, try to test it:

  • A rough order-of-magnitude test is to see how quickly thermal anomalies fade away. The typical diffusion time scale δtδt for a thermal anomaly of size δxδx to fade away to a large extent is δt=δx2κ.
  • A benchmark comparison, by checking your code with those produced by your peers.
  • A resolution test.

A model solution is available here.

3.4. The solution can be made computationally faster by recognizing that the problem is symmetric across the vertical plane x=500 km. So by calculating only half of the solution, the other half is an identical mirror image. The symmetry boundary at x=500 km now becomes the right-hand side boundary of the mesh and should be assigned a different boundary condition.

  • Why is a natural boundary condition the most appropriate one on this boundary?
  • Implement this boundary condition. Implementation of this heat flux through the boundary (i.e. normal to the boundary) is fairly similar to the 1-D case. But keep in mind that, in 2-D, on top of this prescribed boundary-normal heat flux, heat should also diffuse parallel to the boundary at these boundary points.
  • Check that the solution (mirrored across x=500 km) is identical to the one for the full 1000-km wide domain.

3.5. Implement the same natural boundary conditions also on the x=0 boundary and test it in the same way. We will use natural boundary conditions on both side boundaries for Session 4 of this course.

Session3, Part 2: 2-D advection-diffusion

In this session, we will add advection to the heat equation, which will make it suitable for modelling convection. To do so, we will modify the code from the previous session:

3.5. Reduce the size of the thermal anomaly to a 200-by-200 km cold block centred around (xc,zc)=(500,300) km.

3.6. Create a solid rotation velocity field around the centre of the box: v⃗ =(vxvz)=(−vα(z−zc)vα(x−xc)), with vα=2πτvα=2πτ the magnitude of the angular velocity (in rad/sec), rr the distance to the centre of the box (in m), and ττ the convection overturn time (in sec). If you are unsure how to implement this, a possible Python implementation is given here.

3.7. Add the advection part to your twoDdiff subfunction (which you now might want to call twoDadvdiff). Follow the following ‘pseudo-code’:

# Horizontal advection: 
#    loop over all x-points (internal ones only, boundary points have vx=0)
#       loop over all z-points (including boundary points)
#           if vx>0:
#               dTdt_x = -vx / dx * (T_j,i - T_j,i-1)
#           else: 
#               dTdt_x = -vx / dx * (T_j,i+1 - T_j,i)
#           add dTdt_x to diffusion dTdt part
# Vertical advection: 
#    repeat same procedure ...

If you are unsure how to implement this, a possible implementation of this is given here.

A model solution is available here.