% Exercise 3_1, Basin analysis
% McKenzie model
% Sonia Scarselli ETH-Zurich, sonia.scarselli@erdw.ethz.ch
%This program calculates the syn-rift subsidence, the thermal (postrift) subsidence and the total subsidence, ...
...and plots the total subsidence as a function of time.
% The initial parameters can be changed to evaluate different synrift subsidence and thermal subsidence
% curve

clear;

% Define 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 celcius   [kg/m^3]
rho_c0  = 2800;     % Density of the crust at 0 degrees celcius    [kg/m^3]
rho_s   = 2066;     % Density of sediments                         [kg/m^3]
alpha_v = 3.28e-5;  % volumetric coefficient of thermal expansion  [1/K]
Tm      = 1333;     % Temperature of the mantle                    [C]
kappa   = 1e-6;     % Thermal diffusivity                          [m^2/s]
time_my = 0:150;    % Time                                         [my]
time_s  = time_my*365*24*3600*1e6;        % Time in seconds
beta    = 3 ;       % stretch factor
    
% Step 1: Calculate synrift subsidence
ys      = y_l*((rho_m0-rho_c0)*y_c/y_l*(1-alpha_v*Tm*y_c/y_l)-alpha_v*Tm*rho_m0/2)*(1-1/beta)/(rho_m0*(1-alpha_v*Tm)-rho_s)

% Step 2: Calculate thermal subsidence with time
E0      = 4*y_l*rho_m0*alpha_v*Tm/(pi^2*(rho_m0-rho_s));
tau     = y_l^2/(pi^2*kappa);

S       = E0*beta/pi*sin(pi/beta)*(1-exp(-time_s/tau));

% Calculate total subsidence
S_total   = ys + S;          %Total   subsidence in meters
S_thermal = S;               %Thermal subsidence in meters

S_total   = S_total  /1e3;     %Total   subsidence in km's
S_thermal = S_thermal/1e3;     %Thermal subsidence in km's

% Plotting
plot(time_my,S_total,'r-')      
title('Synrift and thermal subsidence, part a')
xlabel('Time since end of rifting [My]')
ylabel('Thermal subsidence [km]');
legend(num2str(beta))
grid on
axis ij