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
import numpy.matlib
from scipy.linalg import solve
from scipy.sparse import spdiags
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
# setup coefficient matrix with flux boundary condition at the bottom
data  = (np.ones((ny,1))*np.array([-beta, (1+2*beta), -beta ])).T
diags = np.array([-1, 0, 1])
A     = spdiags(data, diags, ny, ny).toarray()

# Boundary conditions
# constant temperature at seafloor
A[ny-1,ny-1] = 1
A[ny-1,ny-2] = 0

# constant gradient at the bottom
A[0,(0,1)]  = [(1+2*beta), -2*beta]

# setup coefficient matrix with flux boundary condition at the bottom
# data  = (np.ones((ny,1))*np.array([-beta, (1+2*beta), -beta ])).T
# diags = np.array([-1, 0, 1])
# A     = spdiags(data, diags, ny, ny).toarray()

# Set boundary conditions!!
# constant temperature at seafloor
# A[???] = ???
# A[???] = ???

# constant gradient at the bottom
# A[???]  = [???]

A[0:6,0:6]

array([[ 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]])

# Implicit FD solve with flux boundary conditions
def fdm_solve(Rhs):
    Tnew=solve(A,Rhs)
    return Tnew

# 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 ($^{\circ}$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