Oceanic Geotherm
The 1-D temperature profile of an oceanic geotherm can be calculated by solving the conductive part of the 1-D temperature conservation equation using variable thermal parameters with a proper conserving finite difference scheme (so far only including a radiogenic heat source). For the sake of continuity, we use the 1-D solver for variable thermal parameters, even though, we assume a constant thermal conductivity in this example.
A proper conservative finite difference scheme in 1-D means that the temperature is defined on the centroids and the vertical heat flux $q_y$ and the thermal conductivity $k$ on the vertices.
The 1-D temperature equation is given by:
\[\begin{equation} \rho c_{p} \frac{\partial{T}}{\partial{t}} = \frac{\partial{q_y}}{\partial{y}} + \rho H, \end{equation}\]
where the heat flux is defined as:
\[\begin{equation} \left. q_{y,m} = -k_m \frac{\partial T}{\partial y}\right\vert_{m},\ \textrm{for}\ m = 1:nv, \end{equation}\]
with $\rho, c_{p}, T, t, k, H, y$ and $nv$ are the density [kg/m³], the specific heat capacity [J/kg/K], the temperature [K], the time [s], the thermal conductivity [W/m/K], the heat generation rate per mass [W/kg], the depth [m], and the number of vertices, respectively.
For more details on how to discretize the equation using an explicit, forward Euler finite difference scheme, please see the documentation.
First, one needs to load the required packages:
using Plots, SpecialFunctions
using GeoModBox.HeatEquation.OneDLet's start with the definition of the geometrical, numerical, and physical constants:
# Constants --------------------------------------------------------- #
H = 200e3 # Hight of the model [ m ]
nc = 200 # Number of central grid points
Δy = H/nc # Grid resolution
# Depth ---
yc = LinRange(-H + Δy/2.0,0.0 - Δy/2.0,nc)
yv = LinRange(-H,0,nc+1)
Py = (
ρm = 3000, # Density [ kg/m^3 ]
cpm = 1000, # Heat capacity [ J/kg/K ]
km = 3.0, # Conductivity [ W/m/K ]
HM = 0, # Heat generation rate [W/kg]; Q = ρ*H0
)
# ------------------------------------------------------------------- #In the following, one needs to define the initial and boundary condition:
- Temperature at the surface and bottom.
- Linear increasing temperature profile assuming a certain adiabatic gradient and potential mantle temperature.
# Initial Condition ------------------------------------------------- #
T = (
Tpot = 1315 + 273.15, # Potential temperautre [ K ]
ΔTadi = 0.5, # Adiabatic temperature gradient [ K/km ]
Ttop = 273.15, # Surface temperature [ K ]
T_ex = zeros(nc+2,1),
)
T1 = (
Tbot = T.Tpot + T.ΔTadi*abs(H/1e3), # Bottom temperature [ K ]
T = T.Tpot .+ abs.(yc./1e3).*T.ΔTadi, # Initial T-profile [ K ]
)
T = merge(T,T1)
Tini = zeros(nc,1)
Tini .= T.T
T.T_ex[2:end-1] .= T.T
# ------------------------------------------------------------------- #One can choose between Dirichlet or Neumann thermal boundary conditions at the surface and bottom.
# Boundary conditions ----------------------------------------------- #
BC = (
type = (N=:Dirichlet, S=:Dirichlet),
val = (N=T.Ttop[1],S=T.Tbot[1])
)
# If Neumann boundary conditions are choosen, the following values result in the given heatflux for the given thermal conductivity k.
# S = -0.03; # c = -k/q -> 90 mW/m^2
# N = -0.0033; # c = -k/q -> 10 mW/m^2
# ------------------------------------------------------------------- #Now, one needs to define the multiplication factor fac of the diffusion stability criterion for the explicit thermal solver. This factor controls the stability criterion. If it is larger then 1 the solver becomes instable.
# Time stability criterion ------------------------------------------ #
fac = 0.8 # Courant criterion
tmax = 60 # Lithosphere age [ Ma ]
tsca = 60*60*24*365.25 # Seconds per year
age = tmax*1e6*tsca # Age in seconds
# ------------------------------------------------------------------- #Let's check that the initial and boundary conditions are properly defined by plotting the temperature profile.
# Plot Initial condition -------------------------------------------- #
plotparam = 1
q = plot(Tini,yc./1e3,
label="",
xlabel="T [ K ]", ylabel="z [ km ]",
title="Initial Temperature",
xlim=(T.Ttop,T.Tbot),ylim=(-H/1e3,0))
display(q)
# ------------------------------------------------------------------- #
Figure 1. Initial temperature profile.
Since a thermal solver for variable thermal parameters is used, one needs to expand the scalar to a vector with the dimensions of the number of centroids nc. Additionally, the thermal diffusivity $\kappa$ and initialize the vertical heat flux q need to be defined.
# Setup Fields ------------------------------------------------------ #
Py1 = (
ρ = Py.ρm.*ones(nc,1),
cp = Py.cpm.*ones(nc,1),
k = Py.km.*ones(nc+1,1),
H = Py.HM.*ones(nc,1)
)
Py = merge(Py,Py1)
Py2 = (
# Thermal diffusivity [ m^2/s ]
κ = maximum(Py.k)/minimum(Py.ρ)/minimum(Py.cp),
)
Py = merge(Py,Py2)
T2 = (
q = zeros(nc+1,1),
)
T = merge(T,T2)
# ------------------------------------------------------------------- #Now, one can calculate the time stability criterion.
# Time stability criterion ------------------------------------------ #
Δtexp = Δy^2/2/Py.κ # Stability criterion for explicit
Δt = fac*Δtexp # total time step
nit = ceil(Int64,age/Δt) # Number of iterations
time = zeros(1,nit) # Time array
# ------------------------------------------------------------------- #Now all required parameter and constants are defined to solve the 1-D temperature equation for each time step in a for loop.
The temperature conservation equation is solved via the function ForwardEuler1D!(), which updates the temperature profile T.T for each time step using the extended temperautre field T.T_ex, which include the ghost nodes. The temperature profile is plotted for a certain time.
# Calculate 1-D temperature profile --------------------------------- #
count = 1
for i = 1:nit
if i > 1
time[i] = time[i-1] + Δt
elseif time[i] > age
Δt = age - time[i-1]
time[i] = time[i-1] + Δt
end
ForwardEuler1D!(T,Py,Δt,Δy,nc,BC)
if i == nit || abs(time[i]/1e6/tsca - count*5.0) < Δt/1e6/tsca
println(string("i = ",i,", time = ", time[i]/1e6/tsca))
q = plot!(T.T,yc./1e3,
label=string("t = ",ceil(time[i]/1e6/tsca),"[Ma]"),
xlabel="T [ K ]", ylabel="z [ km ]",
title="Oceanic Geotherm",
xlim=(T.Ttop,T.Tbot),ylim=(-H/1e3,0))
display(q)
count = count + 1
end
end
# ------------------------------------------------------------------- #
Figure 2. Evolution of the temperature profile with depth in 5 Ma steps.
For the final time step, a depth profile for the vertical heat flux is calculated. Therefore, one needs to update the temperature at the ghost nodes to calculate the heat flux at the boundary.
# Calculate heaf flow ----------------------------------------------- #
# South ---
T.T_ex[1] = (BC.type.S==:Dirichlet) * (2 * BC.val.S - T.T_ex[2]) +
(BC.type.S==:Neumann) * (T.T_ex[2] - BC.val.S*Δy)
# North ---
T.T_ex[end] = (BC.type.N==:Dirichlet) * (2 * BC.val.N - T.T_ex[nc+1]) +
(BC.type.N==:Neumann) * (T.T_ex[nc+1] + BC.val.N*Δy)
for j=1:nc+1
T.q[j] = -Py.k[j] *
(T.T_ex[j+1] - T.T_ex[j])/Δy
end
# ------------------------------------------------------------------- #Finally, one can caluclate the temperature profile for an oceanic geotherm using the analytical expression of an infinite half-space cooling model for a certain age. The analytical solution is plotted in the final figure, together with the final numerical temperature profile and the heat flux profile.
# Plot -------------------------------------------------------------- #
if BC.type.N==:Dirichlet && BC.type.S==:Dirichlet
Tana = zeros(nc,1)
Tana .= Tini .+ (T.Ttop - T.Tpot).*erfc.(-yc./(2*sqrt(age*Py.κ)))
Tana[end] = T.Ttop
end
p = plot(T.T,yc./1e3,
label=string("t = ",ceil(maximum(time)/1e6/tsca),"[Ma]"),
xlabel="T [ K ]", ylabel="z [ km ]",
title="Oceanic Geotherm",
xlim=(T.Ttop,T.Tbot),ylim=(-H/1e3,0),
layout=(1,2),subplot=1)
if BC.type.N==:Dirichlet && BC.type.S==:Dirichlet
plot!(p,Tana,yc./1e3,
label="T_HSCM",linestyle=:dash,
layout=(1,2),subplot=1)
end
p = plot!(T.q.*1e3,yv./1e3,
label="",
xlabel="q [ mW ]", ylabel="z [m]",
title="Heat Flux",
ylim=(-H/1e3,0),
subplot=2)
display(p)
savefig(p,"./examples/DiffusionEquation/1D/Results/OceanicGeotherm_1D.png")
savefig(q,"./examples/DiffusionEquation/1D/Results/OceanicGeotherm_1D_evolve.png")
# ======================================================================= #
Figure 3. Final temperature and heat flux profiles.