NNESMO Session 1 – The basics: 0-D and 1-D models
Prof Jeroen van Hunen
Department of Earth Sciences
Durham University
Contents
- Practical 1, Part 1: Modelling an Emptying Bathtub
- Practical 1, Part 2: Radioactive Decay
- Practical 1, Part 3: One-dimensional diffusion
- Practical 1, Part 4: Numerical stability
- Practical 1, Part 5: Cooling of the Earth (ADVANCED)
Practical 1, Part 1: Modelling an Emptying Bathtub
1.1) Download the file bath.m, and open it in Matlab. Study it briefly, and run the code to see what it does. Make sure you understand what the code is doing.
1.2) Add the analytical solution from the lecture notes to the code and plot.
1.3) How does the numerical model perform for the emptying bath?
1.4) How does the size of the timestep affect your results?
Practical 1, Part 2: Radioactive Decay

Uranium 235 ( ) decays to Lead 207 ( ), with a half-life
Myrs. In other words, in 704 Myrs, half of the available
decays to
. More generally, each unit of time (say each million years) a constant fraction of
decays to
. The removal of
(to produce
) is therefore directly proportional to the available
itself (check this for yourself!). Let’s simplify notation and call
simply
. Mathematically, this radioactive decay is then written as:

, with being the ‘decay constant’.

NOTE: An equation like this one, containing a variable and its derivative, is called a differential equation (henceforth referred to as DE) It nicely describes the physics of (in this case) radioactive decay, but it doesn’t provide an immediate solution for the unknown variable ; we have to solve it first. In some cases, including this one, it is possible to solve the DE analytically — in other cases, a numerical approximation must be used.

Quantities (or rather concentrations) of and are usually measured with respect to the concentration of a stable isotope of the daughter element, in this case
. At the start of Earth’s life at 4.567 Ga, the concentrations in the mantle of
and
relative to
were
and
.

1.5) The analytical solution to the DE for radioactive decay is given by . What is the physical meaning of ?

1.6) The relationship between the decay constant and the half-life can be calculated using the analytical solution. Given that the half-life is the time when half of the original radionuclide remains, find
as a function of
. What are the units of
?
1.7) Before programming a numerical model, we first need to translate a (geological?) model into a numerical model. Try to answer the following questions:
- What is/are the key mathematical equation(s) we need to solve?
- Are there any constants you will need to know? Can you look them up or estimate them?
- What will the model look like at time zero? This is the initial condition.
- Roughly for how much (model) time will you need to run our model for?

1.8) Rewrite the differential equation in a Forward Euler finite difference format. Rearrange the equation such that the unknown new quantity of appears on the lefthand side and all known quantities appear on the righthand side of the equation.

1.9) Finally, write the Matlab script to calculate the radioactive decay of and the growth of in a closed
system from the origin of the Earth to the present day using a Forward Euler method. Build your code in a logical and structured way: set up the sections you think you’ll need, then fill them in one at a time, testing as you go. Plot your results and check that they behave the way you expect. Hint: Rather than referring to ‘old’ and ‘new’ solutions, it is easier to use a timestep counter it, so that the ‘new’ solution refers to the value at timestep it and the ‘old’ solution refers to the value at timestep it-1.

1.10) Add the analytical solutions for both and to your code, and modify your plots to show the numerical and analytical solutions for both
and
. How do the numerical solutions compare to the analytical ones?
1.11) Try different timesteps to check that the numerical solution approaches the analytical one for decreasing timesteps. This is called a ‘resolution test’.
1.12) Repeat steps 1.8 to 1.11 for the Backward Euler method and add it to your existing code. How do the numerical and analytical solutions compare now?
1.13) Repeat steps 1.8 to 1.11 for the Crank-Nicholson method and add it to your existing code.

1.14) Calculate the error between the analytical and numerical solutions for for each of the three methods. Add a third section to your plot showing the error through time. Which of the three methods has the smallest error for a given size timestep and why?
Practical 1, Part 3: One-dimensional diffusion

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 , where the partial derivative symbol indicates that in the equation
has derivatives to more than one variable (in this case to time
and depth
). 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 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.

1.15) Sketch (on paper) this temperature distribution as a function of depth for time
. This temperature distribution will be a good initial condition for our modelling exercise.

1.16) 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 at and
km?
1.17) What does the forward-Euler finite difference form of the diffusion equation look like?
![$i =[1,2,3,4,5]$](https://community.dur.ac.uk/jeroen.van-hunen/NNESMO2019/Practicals/Practical1_eq12684994882133834298.png)
1.18) Suppose we have a system, in which the mantle column is described by just 5 grid points (i.e. ). At the boundaries ( and
) boundary conditions apply, and at the remaining, interior points, the above finite difference equation applies. Write out (again on paper, not in Matlab) for each point
the equation that calculates the temperature at the new time step
, where
refers to the spatial dimension and
to the time step. If the old solution
is known for all
, (and
,
, and
are all known as well), we then have all the necessary information to calculate all
. Perform this calculation manually for this 5-point system for one time step. Assume the first solution to be the initial condition you obtained in 1.15, and assume
s,
Myr,
km. Note that different units of time and space are used here. Draw this solution onto the schematic plot of the initial condition. Now apply a second time step and again draw the solution to the plot.
1.19) Now copy the Matlab template heat1d.m, and add the missing functions. Run the code, and check how well the numerical solution corresponds to the analytical one.
If time permits, continue with the following exercise.

1.20) We will define it thermal lithosphere as mantle material cooler than C, 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 C, 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 shallowest grid point for which
C.
Practical 1, Part 4: Numerical stability
1.21) 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 emperical-found maximum stable timestep.
1.22) 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 that the one without any spatial variation (such as our radioactive decay problem), and will be discussed later.
Practical 1, Part 5: Cooling of the Earth (ADVANCED)

The (potential) temperature T of the Earth’s upper mantle today is about , but heat released from accretion of the Earth and differentiation into a core-mantle layered system probably heated the early internal Earth to several hotter than that. Since then, radio-active decay added (and still adds) another significant amount of heat, while the only way for the Earth to lose its heat is through its surface. A simple heat balance for the Earth states that the difference between the total radioactive heat source for the whole Earth,
, and the total heat sink for the Earth,
, (surface heat flow) leads to cooling (or heating) of the Earth. Mathematically, this can be expressed as:
, with
the Earth’s total heat capacity.
,
, and
generally are all functions of time
.

1.23) The following table lists the 4 most important decay systems for the Earth’s thermal history, using , with related to the half-life of a decay system as
:


Give an analytical expression for the total H(t) for the whole Earth through time. depends on past surface tectonics, which is less well known: today, plate tectonics with plate speeds of 5-10 cm/yr is responsible for TW, but whether this value is valid for the Earth’s whole history is unclear. Maybe plate tectonics was slower or faster (or even non-existent) in the past. To start with, we will assume that
was always 36 TW. Note that since
is not dependent on
, this not a d.e., and can be easily solved analytically, provided
and
are simple functions of time.
1.24) Construct a finite difference code to calculate the thermal history of the Earth. Go through all the necessary steps to arrive at a numerical model:
- What is the governing equation for this problem?
- Discuss the initial condition.
- How to discretise your problem?
- Write the numerical code.
- Hint:
can be chosen to be either today (in which case time runs backward) or at the time of Earth’s accretion (so that time runs forward). Note that this forward or reverse sense of time has nothing to do with Forward or Backward Euler timestepping.

1.25) Plot the resulting temperature as a function of time. According to this (simple) model, did the Earth always cool down over time?
1.26) Test the code with a convergence test, and compare results with those of your colleagues.
If time permits, implement the following adaptations.

1.27) The Urey ratio is defined as the present-day ratio of , and is thought to be anywhere between
and
. Vary
by multiplying the total heat production with a constant prefactor, and explore how that affects the Earth’s thermal evolution. Plot multiple
results in a single plot.

1.28) is probably not constant, but mantle temperature dependent (e.g. a hotter, weaker mantle might lead to faster mantle convection and therefore faster cooling). Adjust to be linearly proportional to
. Note that this will make the governing equation a d.e. (since now
is dependent on
). Examine the different solution for Forward Euler, Backward Euler, and Crank-Nicholson timestepping.