%Exercise 3_3, Basin analysis
%Calculation of Palaeotemperatures in a Basin undergoing Thermal Subsidence.
%This program calculates the temperature of a chosen horizon as a function of time since the end of rifting.
% In this exercise the depth of the horizon and the time since since the end of rifting are obtained from...
...a decompaction and backstripping analysis of a borehole or surface stratigraphic section.
% Sonia Scarselli, ETH-Zurich, sonia.scarselli@erdw.ethz.ch

clear;

% Define some initial parameters
y_c     = 35000;                        % Initial crustal thickness in m              [m]
y_l     = 125000;                       % Initial lithospheric thickness in m         [m]
rho_m0  = 3330;                         % Density of the mantle at 0 degrees celsius  [kg/m^3]
rho_c0  = 2800;                         % Density of the crust at 0 degrees celsius   [kg/m^3]
rho_s   = 1000;                         % Density of the infilling material        	  [kg/m^3]
alpha_v = 3.28e-5;                      % Thermal expansion coefficient               [1/K]
Tm      = 1333;                         % Temperature of the mantle                   [C]
kappa   = 1e-6;                         % Thermal diffusivity                         [m^2/s]
Ks      = 1.25;                         % Thermal conductivity of sediments           [W/m/K]
Kb      = 3.0;                          % Thermal conductivity of basement            [W/m/K]
time_my = [0, 35, 45, 100];             % Time since end of rifting                   [my]
time_s  = time_my*365*24*3600*1e6;      % Time since end of rifing in seconds
depth   = [0.851, 2.098, 2.524, 3.976]; % Depth of in km
depth   = depth*1000;                   % Subsidence in m
beta1    =1.5 ;                         % Beta factor
beta2 = 2;
beta3 =4;

tau     = y_l^2/(pi^2*kappa);

for i=1:length(depth)
    y    = depth(i);
    t    = time_s(i);
    q1(i) = Kb*Tm/y_l*(1+2*beta1/pi*sin(pi/beta1)*exp(-t/tau))
    T1(i) = Kb/Ks*Tm*y/y_l*(1+2*beta1/pi*sin(pi/beta1)*exp(-t/tau))
    q2(i) = Kb*Tm/y_l*(1+2*beta2/pi*sin(pi/beta2)*exp(-t/tau))
    T2(i) = Kb/Ks*Tm*y/y_l*(1+2*beta2/pi*sin(pi/beta2)*exp(-t/tau))
    q3(i) = Kb*Tm/y_l*(1+2*beta3/pi*sin(pi/beta3)*exp(-t/tau))
    T3(i) = Kb/Ks*Tm*y/y_l*(1+2*beta3/pi*sin(pi/beta3)*exp(-t/tau))
end

plot(time_my,T1,'b+-');
hold on
plot(time_my,T2,'y+-');
hold on
plot(time_my,T3,'g+-');
hold off
xlabel('Time since end of rifting    [Ma]');
ylabel('Temperature at the base of the lower Cretaceaous [C]');
grid on