% Exercise 9_1, Basin Analysis
% It is possible to reconstruct the decompacted layers, the changes of the initial parametrs
% representing the thicknesses derived by boreholes data permit to evaluate the decompacted
% layers in every natuaral example
%Sonia Scarselli, ETH-Zentrum, sonia.scarselli@erdw.ethz.ch

%--------------------------------------------------------------------------
clear
% initial parameters ( actual thicknesses from boreholes)
z_PreC= 1.052; % actual thickness of the PreCreatceous shaley sandstones , in Km
z_LowC= 0.459; % actual thickness of the Lower Cretaceous shales , in Km
z_UpC = 0.968; % actual thickness of the Upper Cretaceous chalk, in Km
z_Pal = 0.605; % actual thickness of the Palaeocene sandstones , in Km
z_EoPle= 1.944; % actual thickness of the Eocene Pleistocene shales , in Km
% Surface porosity for each formation
fi_PreC0 = 0.56; % surface porosity of the shaley sandstones facies;
fi_LowC0 = 0.63;% surface porosity of the shales facies;
fi_UpC0 =0.70;% surface porosity of the chalk facies;
fi_Pal0 = 0.49;% surface porosity of the sandstones facies;
fi_EoPl0 =0.63; % surface porosity of the shales facies;
% fu-depth coefficent for each formation
c1= 0.39; % fi-depth coefficent in km-1 for shaley sandstones;
c2= 0.51;% fi-depth coefficent in km-1 for shales;
c3 = 0.71;% fi-depth coefficent in km-1 for chalk;
c4 = 0.27;% fi-depth coefficent in km-1 for sandstones;
c5 =c2;
y6 = 0;% actual top of the layer5, surface
y5 = 1.944; % actual top of the layer4 and bottom of the layer 5, in Km
y4 = 2.549; % actual top of the layer3 and bottom of the layer 4, in Km
y3 = 3.517; % actual top of the layer2 and bottom of the layer3, in Km
y2 = 3.976 ; % actual top of the layer1 and bottom of the layer 2, in Km
y1 = 5.028; % actual bottom of the layer
y = 0; % surface conditions 
% density at the actual condition
rho_w = 1030; % density of the sea water; 
rho_UpC = 2710 ; % density for the layer 4, in kg/m3;
rho_Pal= 2650 ; % density for the 5th layer;
rho_EoPl= 2720; % density of the sixth layers;
rho_m = 3300; % density of the mantle

% CALCULATION OF THE DECOMPACTION FOR EACH FORMATION

% decompaction of the LAYER 1, 210-140MA (T1); x1_1 is the thickness of the layer 1 at time T1
eps = 0.0001 ;
error = 2*eps ;
 x1_1 =1;
while error > eps
    xPreC1=  z_PreC-fi_PreC0/c1*(exp(-c1*y2)-exp(-c1*y1))+ fi_PreC0/c1*(exp(-c1*y)-exp(-c1* x1_1));
    error = abs(xPreC1-x1_1 )/ x1_1
     x1_1 = xPreC1
end


% decompaction of the LAYER2, 140-100 MA (T2); x2_2 is the thickness of the layer 2 at time T2
eps = 0.0001 ;
error = 2*eps ;
x2_2=1;
while error > eps
    xLowC2=  z_LowC-fi_LowC0/c2*(exp(-c2*y3)-exp(-c2*y2))+ fi_LowC0/c2*(exp(-c2*y)-exp(-c2*x2_2));
    error = abs(xLowC2-x2_2)/x2_2
    x2_2 = xLowC2    
end


 % decompaction of the LAYER 1, 140-100MA (T2), x1_2is the  thickness of the layer 1 at time T2
eps = 0.0001 ;
error = 2*eps ;
x1_2 =1;
while error > eps
    xPreC2=  (z_PreC-fi_PreC0/c1*(exp(-c1*y2)-exp(-c1*y1))+ fi_PreC0/c1*(exp(-c1*xLowC2)-exp(-c1*x1_2))) + xLowC2; 
    error = abs(xPreC2-x1_2)/x1_2
    x1_2 = xPreC2
end


% decompaction of the LAYER3, 100-65 MA (T3), x3_3 is the thickness of the layer 3 at time T3
eps = 0.0001 ;
error = 2*eps ;
x3_3 =1;
while error > eps
    xUpC3=  (z_UpC-fi_UpC0/c3*(exp(-c3*y4)-exp(-c3*y3))+ fi_UpC0/c3*(exp(-c3*y)-exp(-c3*x3_3))) 
    error = abs(xUpC3-x3_3)/x3_3
    x3_3 = xUpC3
end


% decompaction of the LAYER2 (Lower Cretaceous, 100-65MA (T3), x2_3 is the thickness of the layer 2 at time T3
eps = 0.0001 ;
error = 2*eps ;
x2_3= 1;
while error > eps
    xLowC3=  z_LowC-fi_LowC0/c2*(exp(-c2*y3)-exp(-c2*y2))+ fi_LowC0/c2*(exp(-c2*xUpC3)-exp(-c2*x2_3)) + xUpC3 ;
    error = abs(xLowC3-x2_3)/x2_3
    x2_3= xLowC3    
end

 % decompaction of the LAYER 1, 100-65MA (T3), x1_3 is the thickness of the layer 1 at time T3
eps = 0.0001 ;
error = 2*eps ;
x1_3 =1;
while error > eps
    xPreC3=  (z_PreC-fi_PreC0/c1*(exp(-c1*y2)-exp(-c1*y1))+ fi_PreC0/c1*(exp(-c1*xLowC3)-exp(-c1*x1_3))) + xLowC3; 
    error = abs(xPreC3-x1_3)/x1_3
    x1_3 = xPreC3
end

% decompaction of the layer 4, 65-55 MA(T4), x4_4 is the thickness of the layer 4 at time T4
eps = 0.0001 ;
error = 2*eps ;
x4_4 =1;
while error > eps
    xPal4 = z_Pal-fi_Pal0/c4*(exp(-c4*y5)-exp(-c4*y4))+ fi_Pal0/c4*(exp(-c4*y)-exp(-c4*x4_4));
    error = abs( xPal4 -x4_4 )/x4_4 
    x4_4  =  xPal4 
end

% decompaction of the LAYER3, 65-55 MA (T3), x3_4 is the thickness of the layer 3 at time T4
eps = 0.0001 ;
error = 2*eps ;
x3_4 =1;
while error > eps
    xUpC4=  (z_UpC-fi_UpC0/c3*(exp(-c3*y4)-exp(-c3*y3))+ fi_UpC0/c3*(exp(-c3*xPal4 )-exp(-c3*x3_4))) + xPal4 ;
    error = abs(xUpC4-x3_4)/x3_4
    x3_4 = xUpC4
end

% decompaction of the LAYER2 (Lower Cretaceous, 65-55MA (T4), x2_4 is the thickness of the layer 2 at time T4
eps = 0.0001 ;
error = 2*eps ;
x2_4= 1;
while error > eps
    xLowC4=  z_LowC-fi_LowC0/c2*(exp(-c2*y3)-exp(-c2*y2))+ fi_LowC0/c2*(exp(-c2*xUpC4)-exp(-c2*x2_4)) + xUpC4 ;
    error = abs(xLowC4-x2_4)/x2_4
    x2_4= xLowC4    
end

% decompaction of the LAYER 1, 65-55MA (T4), x1_4 is the thickness of the layer 1 at time T4
eps = 0.0001 ;
error = 2*eps ;
x1_4 =1;
while error > eps
    xPreC4=  (z_PreC-fi_PreC0/c1*(exp(-c1*y2)-exp(-c1*y1))+ fi_PreC0/c1*(exp(-c1*xLowC4)-exp(-c1*x1_4))) + xLowC4; 
    error = abs(xPreC4-x1_4)/x1_4
    x1_4 = xPreC4
end

% decompaction of the layer 5 , 55-0 Ma (T5), x5_5 is the thickness of the layer 5 at time T5
eps = 0.0001 ;
error = 2*eps ;
x5_5 =1;
while error > eps
    xEoPl5=  z_EoPle- fi_EoPl0/c5*(exp(-c5*y6)-exp(-c5*y5))+ fi_EoPl0/c5*(exp(-c5*y)-exp(-c5*x5_5)) ; 
    error = abs(xEoPl5-x5_5)/x5_5
   x5_5 = xEoPl5
end

% decompaction of the layer 4, 55-0MA(T4), x4_4 is the thickness of the layer 4 at time T5
eps = 0.0001 ;
error = 2*eps ;
x4_5=1;
while error > eps
    xPal5 = ( z_Pal-fi_Pal0/c4*(exp(-c4*y5)-exp(-c4*y4))+ fi_Pal0/c4*(exp(-c4*xEoPl5)-exp(-c4*x4_5))) + xEoPl5;
    error = abs( xPal5 -x4_5 )/x4_5
    x4_5 =  xPal5
end

% decompaction of the LAYER3,55-0 MA (T5), x3_5 is the thickness of the layer 3 at time T5
eps = 0.0001 ;
error = 2*eps ;
x3_5 =1;
while error > eps
    xUpC5=  (z_UpC-fi_UpC0/c3*(exp(-c3*y4)-exp(-c3*y3))+ fi_UpC0/c3*(exp(-c3*xPal5 )-exp(-c3*x3_5))) + xPal5 ;
    error = abs(xUpC5-x3_5)/x3_5
    x3_5 = xUpC5
end

% decompaction of the LAYER2 (Lower Cretaceous, 55-0MA (T5), x2_5 is the thickness of the layer 2 at time T5
eps = 0.0001 ;
error = 2*eps ;
x2_5= 1;
while error > eps
    xLowC5=  z_LowC-fi_LowC0/c2*(exp(-c2*y3)-exp(-c2*y2))+ fi_LowC0/c2*(exp(-c2*xUpC5)-exp(-c2*x2_5)) + xUpC5 ;
    error = abs(xLowC5-x2_5)/x2_5
    x2_5= xLowC5    
end

% decompaction of the LAYER 1, 55-0MA (T5), x1_5 is the thickness of the layer 1 at time T5
eps = 0.0001 ;
error = 2*eps ;
x1_5 =1;
while error > eps
    xPreC5=  (z_PreC-fi_PreC0/c1*(exp(-c1*y2)-exp(-c1*y1))+ fi_PreC0/c1*(exp(-c1*xLowC5)-exp(-c1*x1_5))) + xLowC5; 
    error = abs(xPreC5-x1_5)/x1_5
    x1_5= xPreC5
end

%--------------------------------------------------------------------------
% CALCULATION OF BACKSTRIPPED POROSITY FOR EACH FORMATION
fi_UpC3 =  fi_UpC0/c3 *(((exp(-c3*y)) - exp(-c3*x3_3))/(x3_3-y)) % porosity of the Layer 3 at T3, 
fi_Pal4 =  fi_Pal0/c4 *(((exp(-c4*y)) - exp(-c4*x4_4))/(x4_4-y)) % porosity of the Layer 4 at T4, 
fi_UpC4 =  fi_UpC0/c3 *(((exp(-c3*x4_4)) - exp(-c3*x3_4))/(x3_4-x4_4)) % porosity of the Layer 3 at T4, 
fi_EoPl5 =  fi_EoPl0/c5 *(((exp(-c5*y)) - exp(-c5*x5_5))/(x5_5-y)) % porosity of the Layer 5 at T5, 
fi_Pal5=  fi_Pal0/c4 *(((exp(-c4*x5_5)) - exp(-c4*x4_5))/(x4_5-x5_5)) % porosity of the Layer 4 at T5,
fi_UpC5=  fi_UpC0/c3 *(((exp(-c3*x4_5)) - exp(-c3*x3_5))/(x3_5-x4_5)) % porosity of the Layer 3 at T3, 

%--------------------------------------------------------------------------
% INCREMENT OF THE DENSITY REFFERED TO THE FORMATIONS DURING THE TIME
% increment of the bulk density from 100Ma to 0 Ma
rho_t4= ((fi_UpC3*rho_w  + (1-fi_UpC3)*rho_UpC)/x3_3)* x3_3% increment due to the fourth layer;
% increment due to the fifth layer at time 5, 55MA
rho_t5 = (fi_UpC4*rho_w  + (1-fi_UpC4)*rho_UpC)*((x3_4 - x4_4)/x3_4) + (fi_Pal4*rho_w  + (1-fi_Pal4)*rho_Pal)*(x4_4/x3_4)
% increment due to the sixth layer at time 6, 0MA
rho_t6 = (fi_UpC5*rho_w  + (1-fi_UpC5)*rho_UpC )*((x3_5 - x4_5) /x3_5) + (fi_Pal5*rho_w  + (1-fi_Pal5)*rho_Pal)*((x4_5 - x5_5)/x3_5) + (fi_EoPl5*rho_w  + (1-fi_EoPl5)*rho_EoPl)*(x5_5/x3_5)

%BACKSTRIPPED TECTONIC SUBSIDENCE USING AIRY MODEL
Y0 = 0 ; % the subsidence started in the time 100 MA, so the thickness at that time is zero
Y1 = x3_3 * ((rho_m - rho_t4)/(rho_m-rho_w))
Y2 =x3_4 * ((rho_m - rho_t5)/(rho_m-rho_w))
Y3 = x3_5 * ((rho_m - rho_t6)/(rho_m-rho_w))

% plot of the data
subplot(1,2,1)
tPreCR = [210;140;  100; 65; 55;  0]; % plot of the decompacted sequences
YPreCR = [ Y0; x1_1;  x1_2;  x1_3; x1_4;  x1_5];
plot(tPreCR ,YPreCR,'bo-')
hold on;
tLowCR = [140;  100; 65; 55;  0];
YLowCR = [Y0; x2_2;  x2_3; x2_4; x2_5];
plot(tLowCR,YLowCR,'bo-')
hold on;
tUpCR= [  100; 65; 55;  0];
YUpCR = [ Y0;  x3_3; x3_4; x3_5];
plot(tUpCR ,YUpCR,'bo-')
hold on;
tPa = [ 65; 55;  0];
YPa = [Y0;x4_4; x4_5];
plot(tPa ,YPa,'bo-')
hold on;
tEoPl = [ 55;  0];
YEoPl = [ Y0;x5_5 ];
plot(tEoPl ,YEoPl,'bo-')
hold on;
title('Decompacted Depth versus Time');
xlabel('Time,My');
 ylabel ('Depth, Km');
set(gca,'XGrid','on','YGrid','on');       
 set(gca,'XDir','reverse','YDir','reverse');
hold off;
tTOT = [ 100; 65; 55;  0]; % plot of the tectonic subsidence
YTOT = [ Y0; Y1; Y2; Y3];
subplot(1,2,2)
 plot ( tTOT, YTOT,'ro-')
hold on;
title('Tectonic Subsidence versus Time');
xlabel('Time,My');
 ylabel ('Depth, Km');
set(gca,'XGrid','on','YGrid','on');       
 set(gca,'XDir','reverse','YDir','reverse');
 hold off