Periodic variations in seafloor temperature

We will explore how periodic variations in seafloor temperature affect the sub-seafloor temperature structure. For this we will learn about flux boundary conditions

import numpy as np
from scipy.sparse.linalg import factorized
from scipy.sparse import diags, csc_matrix
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML
import matplotlib as mpl

# Physical parameters
L       = 30           # Length of modeled domain [m]
Ttop    = 4            # Temperature at seafloor [C]
kappa   = 1e-6         # Thermal diffusivity of rock [m2/s]
day     = 60*60*24     # seconds per day
hf      = 0.06          # basement heat flow in W/m2
k       = 1.5           # thermal conductivity of sediment
ampl    = 1             # surface temperature variation
zeval   = -5
r_time  = 10*365*day

# Numerical parameters
ny      = 101           # Number of gridpoints in y-direction
steps   = 301;
Time,dt    = np.linspace(0, r_time, steps,  retstep=True)
Yvec,dy = np.linspace(-L, 0, ny,  retstep=True)    # Y coordinate vector, constant spacing

# Boundary and initial conditions
Tsurf   = Ttop + ampl*np.sin(2*np.pi/(365/2*day)*Time)
gradT   = -hf/k
T_init  = Ttop + gradT*Yvec
time    = 0
beta    = kappa*dt/dy**2

# hide: the code in this cell is hidden by the author

/Users/lruepke/miniconda3/envs/py312_fem_class/lib/python3.12/site-packages/scipy/sparse/_index.py:210: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil and dok are more efficient.
  self._set_arrayXarray(i, j, x)

# define diagonals
#lower = -beta * np.ones(ny - 1)
#diag  = (1 + 2*beta) * np.ones(ny)
#upper = -beta * np.ones(ny - 1)

# puzzle into sparse matrix
#A = diags([lower, diag, upper], offsets=[-1, 0, 1], format='csc')

# Boundary conditions
# constant temperature at seafloor
#A[-1, :] = ?; A[-1, -1] = ?

# constant gradient at the bottom
# ???

# factorize matrix A for fast solves
#solve_A = factorized(A)

print(A[0:6,0:6].toarray())

[[ 24.36 -23.36   0.     0.     0.     0.  ]
 [-11.68  24.36 -11.68   0.     0.     0.  ]
 [  0.   -11.68  24.36 -11.68   0.     0.  ]
 [  0.     0.   -11.68  24.36 -11.68   0.  ]
 [  0.     0.     0.   -11.68  24.36 -11.68]
 [  0.     0.     0.     0.   -11.68  24.36]]

def fdm_solve(Rhs: np.ndarray) -> np.ndarray:
    # fast per-timestep solve
    return solve_A(Rhs)

# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=(10,5))
ax = plt.axes(xlim=(Ttop-ampl-1, Ttop-L*gradT+1), ylim=(-L, 0))

# how many lines do you want to display
num_solutions       = 20
solution_list,lines = [],[]
for i in range(0,num_solutions):
    line, = ax.plot([], [], lw=1)
    lines.append(line)
timeLabel=ax.text(0.02,0.98,'Time: ',transform=ax.transAxes,va='top')
ax.set_xlabel('Temperature (°C)')
ax.set_ylabel('Depth (m)')

# Initialize Tnew
Tnew=T_init
# Animation function which updates figure data.  This is called sequentially
def animate(i):
    timeLabel._text='Time: %.1f day'%(i*dt/day)

    # use global keyword to store the latest solution and update it using fdm_solve function
    global Tnew
    Rhs = Tnew.copy()
    Rhs[0]= Rhs[0] + 2*beta*dy*hf/k;
    Rhs[-1]    = Tsurf[i];
    Tnew=fdm_solve(Rhs)
    # Update solution list and always keep latest solutions
    if(len(solution_list)<num_solutions):
        solution_list.append(Tnew)
    elif(len(solution_list)==num_solutions):
        solution_list.pop(0) #remove the oldest resolution
        solution_list.append(Tnew)
    # update lines
    for i in range(0,len(solution_list)):
        lines[i].set_data(solution_list[i], Yvec)
        lines[i].set_color('k')
        lines[i].set_alpha((num_solutions-i)/num_solutions)
# Call the animator.
anim = animation.FuncAnimation(fig, animate,frames=steps,interval=50)
plt.close(anim._fig)

# Call function to display the animation
#HTML(anim.to_html5_video())  # lower resolution
HTML(anim.to_jshtml())  # higher resolution

Once Loop Reflect