Finite-Difference Model for the Thermal History of the Earth

• Download source code - 28.4 KB

Introduction

Assuming that the earth was a hot ball at a homogeneous temperature upon its formation, we can model the cooling of the earth using heat transfer modeling. The model allows to compute the geotherm of the earth, and make predictions about the earth’s crust thickness and the geothermal gradient for the first km of the earth, 4.5 Ga after earth’s formation. For this purpose, we apply the Fourier heat diffusion equation to the earth for the initial and boundary value problem defined by the initial temperature of the earth and the earth’s surface temperature.

The Fourier heat diffusion model expressed in spherical coordinates is as follows:

where T is the temperature, r the radius, t the time, and D the thermal diffusivity.

Consider a hot ball of radius a.

At time zero, our hot ball has a uniform temperature:

The earth's surface temperature is a constant:

The explicit discretization yields T_i+1,j = A*T_i,j+1 + B*T_i,j + C*T_i,j-1 + D, where A, B, C, and D are given in the table below.

To account for the latent heat fusion of magma, or the energy converted into heat during the crystallization of magma, we introduce the enthalpy which is expressed as follows:

where ρ_s and ρ_l the density of the solid and liquid phases, c_s and c_l are the specific heat of solid and liquid phases, h the latent heat of fusion of magma, and T_f the melting point of magma.

As the cooling of the earth starts at a fully liquid phase, we introduce the energy of the latent heat fusion of magma Q_i,j at each node of the lattice, which is initially set to ρ*h. For each incremental calculation, we then check whether T_i+1,j < T_f and Q_i,j > 0. If this is the case, we first decrease Q_i,j by ρ*c*(T_f -T_i+i,j), and if Q_i+1,j is still positive, we set T_i+1,j = T_f; otherwise, T_i+1,j = T_f + Q_i+1,j / (ρ*c).

Using the Code

The class ExplicitScheme.cs does the job. Each cell of the lattice is represented by the class Node.cs.

The method next() of the ExplicitScheme.cs does the traversal of the cells from the surface to the center of the earth to compute the new temperatures for each time increment. The class Node.cs does the calculation for the latent heat fusion of magma which is converted to heat during crystallization of magma.

The radiogenic heating at each time increment considering the radioactive decay of uranium, thorium and potassium isotopes is calculated as follows:

// calculates radiogenic heating at current time in W/kg
q = 3.95e-12 * Math.Exp(-0.158 * time / Ga) + 2.7e-12 * Math.Exp(-0.05 * time / Ga) + 
    5.7e-12 * Math.Exp(-0.577 * time / Ga);

The thermal diffusivity as a function of temperature and earth’s layers is calculated as follows:

// thermal diffusivity function of temperature and layers of the earth
if (TVector[i].Temperature < 363) 
{ 
   // thermal diffusivity of granite layer 
   D = Math.Max((19.64 - 0.0222 * TVector[i].Temperature) * 1e-7, 7.3e-7); 
}
else if (TVector[i].Temperature < 1050) 
{ 
   // thermal diffusivity of basement rocks 
   D = (7.85 - 0.00284 * TVector[i].Temperature) * 1e-7; 
}
else 
   D = 3e-7; //thermal diffusivity of liquid magma

Mathematical details for the algorithm and model parameters are described in the article "A Finite-Difference Model for the Thermal History of the Earth".

Reference

• A Finite-Difference Model for the Thermal History of the Earth