%A program to decompact sediments and then backstrip well A1NC198

clear
n = 10;
% initial parameters ( actual thicknesses from boreholes)
z = zeros(1,n);
z(1) = 252; % Thickness of the Continental Mesozoic, in m
z(2) = 215; % Thickness of the Post Tasselian in m
z(3) = 378; % Dalma
z(4) = 891; % Binem
z(5) = 44; % Tadrart
z(6) = 78; % Akakus
z(7) = 106; % Tannezuft
z(8) = 201; % Mamuniyat
z(9) = 1534; % PreGlacial Cambro-Ordovician
z(10) = 1; % Bottom meter of C-O so that the column starts at zero
% Surface porosity for each formation from Sclater & Christie
fi = zeros(1,n);
fi(1) = 0.50; % surface porosity of the Continental Mesozoics sandstones facies;
fi(2) = 0.51; % surface porosity of the Post Tasselian sandstones facies;
fi(3) = 0.51; % surface porosity of the Dalma formation;
fi(4) = 0.55; % surface porosity of the Binem Formation;
fi(5) = 0.515; % surface porosity of the Tadrart;
fi(6) = 0.55; % surface porosity of the Akakus;
fi(7) = 0.61; % surface porosity of the Tannezuft;
fi(8) = 0.50; % surface porosity of the Mamuniyat;
fi(9) = 0.49; % surface porosity of the PreGlacial Cambro Ordovician
fi(10)= 0.49; % Bottom metre of C-O
% fi-depth coefficent for each formation
c = zeros(1,n);
c(1) = 0.000280; % fi-depth coefficent in m^-1 for sandstones of the Continental Mesozoic;
c(2) = 0.000310; % fi-depth coefficent in m^-1 for sandstones of the Post Tasselian;
c(3) = 0.000310; % fi-depth coefficent in m^-1 for the Dalma Fm;
c(4) = 0.000405; % fi-depth coefficent in m^-1 for the Binem formation;
c(5) = 0.000295; % fi-depth coefficent in m^-1 for the Tadrart formation;

c(6) = 0.000420; % fi-depth coefficent in m^-1 for the Akakus formation;
c(7) = 0.000485; % fi-depth coefficent in m^-1 for the Tannezuft shaley sequence;
c(8) = 0.000280; % fi-depth coefficent in m^-1 for the Mamuniyat formation;
c(9) = 0.000270; % fi-depth coefficent in m^-1 for the PreGlacial Cambro Ordovician formation;
c(10)= 0.00027; % Bottom metre
%depth to the top of the formations
y = zeros(1,n);
y(1) = 0; % Continental Mesozoic Sandstones
y(2) = 252; % Post Tasselian Sandstones
y(3) = 467; % Dalma Fm
y(4) = 845; % Binem Fm
y(5) = 1736; % Tadrart Fm
y(6) = 1780; % Akakus Fm
y(7) = 1858; % Tannezuft shales
y(8) = 1964; % Mamuniyat Fm
y(9) = 2164; % Preglacial Cambro-Ordovician sandstones
y(10)= 3699; % Lowest Metre
yBas = 3700;
% density at the surface conditions
rho_w = 1030; % density of the sea water Kg/m-3;
rho_m = 3300; % density of the mantle
rho = zeros(1,n);
rho(1) = 2655; % grain density of the Continental Mesozoic Sandstones
rho(2) = 2660; % grain density of the Post tasselian Sandstones
rho(3) = 2660; % grain density of the Dalma Fm
rho(4) = 2690; % grain density of the Binem Fm
rho(5) = 2655; % grain density of the Tadrart Fm
rho(6) = 2680; % grain density of the Akakus Fm
rho(7) = 2710; % grain density of the Tannezuft Fm
rho(8) = 2655; % grain density of the Mamuniyat Fm
rho(9) = 2650; % grain density of the Cambro-Ordovician sandstones
rho(10) = 2650; % grain density of the bottom metre
%The paleosea-levels compared to the present day sea-level taken form the
%Exxon curves
sl(1) = 205; % The Cretaceous (Continental Mesozoic)
sl(2) = 20; % The Upper Jurrassic (Post Tasselian)
sl(3) = 20; % Mississippian (Dalma Fm)
sl(4) = 105; % Late Devonian (Binem Fm)
sl(5) = 175; % Early-Mid Devonian (Tadrart Fm)
sl(6) = 105; % Mid-Late Silurian (Akakus Fm)
sl(7) = 145; % Early Silurian (Tannezuft Fm)
sl(8) = 140; % Hirnantian (Mamuniyat Rm)
sl(9) = 165; % Cambro-Ordovician
sl(10) = 100; % Start of the Cambrian
sl = fliplr(sl);
%The Paleobathymetry or water depth.

Wd(1) = 0; % The Cretaceous (Continental Mesozoic)
Wd(2) = 0; % The Upper Jurrassic (Post Tasselian)
Wd(3) = 20; % Mississippian (Dalma Fm)
Wd(4) = 40; % Late Devonian (Binem Fm)
Wd(5) = 0; % Early-Mid Devonian (Tadrart Fm)
Wd(6) = 60; % Mid-Late Silurian (Akakus Fm)
Wd(7) = 100; % Early Silurian (Tannezuft Fm)
Wd(8) = 50; % Hirnantian (Mamuniyat Rm)
Wd(9) = 5; % Cambro-Ordovician
Wd(10) = 5; % Start of the Cambrian
Wd = fliplr(Wd);
x = zeros(n,n); %create an array for the depths to the top of the units over time
for it = 1:n
eps = 0.0001 ;
error = 2*eps ;
x_top = 0;
for in = n-(it-1):n
x_iter = 1;
eps = 0.0001 ;
error = 2*eps ;
while error > eps
x(in,it) = z(in) - fi(in)/c(in)*(exp(-c(in)*y(in))-exp(-c(in)*(y(in)+z(in)))) + fi(in)/c(in)*(exp(-c(in)*x_top)-exp(-c(in)*x_iter)) + x_top;
error = abs(x(in,it)-x_iter)/x_iter;
x_iter = x(in,it);
end
x_top = x(in,it); % reset for the next layer
end
end
xu = flipud(x); % This creates a matrix with time on the x-axis and each layer on the z-axis
t = [513 445.6 443.7 440.8 416 386 359.2 318 223.4 97]; %prepare for the plot
% carry out the backstripping
%calculating the porosity for each layer through time
fi_depth = zeros(n,n);
rho_i = zeros(1,n);
rho_b = zeros(1,n);
tecsub = zeros(1,n);
for it = 1:n
x_top = 0;
for in = n-(it-1):n
fi_depth(in,it) = (fi(in)/c(in)) * ((exp(-c(in)*x_top)-exp(-c(in)*x(in,it)))/(x(in,it) - x_top));
rho_i(in) = (fi_depth(in,it)*rho_w + (1-fi_depth(in,it))*rho(in))*(x(in,it)-x_top)/x(n,it);
x_top = x(in,it);
end
rho_i;
rho_b(it) = sum(rho_i);
tecsub(it) = x(n,it)*((rho_m - rho_b(it))/(rho_m - rho_w)) -sl(it)*(rho_w/(rho_m-rho_w)) + (Wd(it)-sl(it));
end
%Plot of the complete output of the program showing both the decompacted
%subsidence and the tectonic subsidence
figure(1,'position',[0,0,1024,800]), clf
title ('Subsidence of the Kufrah Basin Well A1-NC198')
subplot (2,1,1)
plot (-t,-xu(1,1:n),'bo-',-t,-xu(2,1:n),'bo-',-t,-xu(3,1:n),'bo-',-t,-xu(4,1:n),'kx-',-t,-xu(5,1:n),'gx-',-t,-xu(6,1:n),'rd-',-t,-xu(7,1:n),'rd-',-t,-xu(8,1:n),'c^-',-t,-xu(9,1:n),'y^-',-t,-xu(10,1:n),'y^-','LineWidth',2)
set (gca,'FontSize',14)
title ('The Decompacted Sedimetary thickness for the Kufrah basin','FontSize',18)
xlabel ('Time (Ma)','FontSize',16)
ylabel ('Depth of each uncompacted layer (m)','FontSize',16)
%legend ('Cambro-Ordivician','off','off','off','off','off','Black Shale','Silurian','hide','hide','hide','Devonian','hide','hide','hide','hide','Carboniferous','location','SouthWest')
subplot (2,1,2)
plot (-t,-tecsub,'o-','LineWidth',2)
set (gca,'FontSize',14)
title ('Water loaded tectonic subsidence','FontSize',18)
xlabel ('Time (Ma)','FontSize',16)
ylabel ('The Tectonic subsidence (m)','FontSize',16)
axis ([-550 0 -2000 0])