Here is another way to determine the temporal change of heat content in the substratum.
rho*cp*dT/dt = d/dx (k*dT/dx)
Integrating with respect to x gives
d/dt (integral_{x=-20}^{x=0} (rho*cp*T) dx ) = k*dT/dx @h=0 - k*dT/dx @h=-20
Above, we tried to approximate the right-hand side. Here, we try to approximate the left-hand side.
close all
clc
load("time_in_days.mat")
load("Depth_Profiles.mat")
load("Temperature_Profiles.mat")
num_intervals = size(Temperature_Profiles, 2);
depth_index = find(Depth_Profiles <= -20, 1, 'last');
rhocp = 1.0; % Specify rho*cp of the substratum [J/(m^3*K)]
for i = 1:num_intervals
temperature_profile = Temperature_Profiles(:, i);
energy_content(i) = trapz(Depth_Profiles(depth_index:end),temperature_profile(depth_index:end))*rhocp; %[J/m^2]
end
figure(1)
plot(time_in_days,energy_content) %[J/m^2]
change_of_energy_content = gradient(energy_content)./gradient(time_in_days*24*3600);%[W/m^2]
figure(2)
plot(time_in_days,change_of_energy_content*1e3) %[mW/m^2]
