% Conductivity Excercise

% Insert input parameteres (columns: Depth, porosity, b, conductivity; rows:
% roks types (sand, siltstone, clay, marl,limestone)

Tab=[1300,0.35,0.284,1.79;3500,0.30,0.379,1.58;4800,0.25,0.293,1.43;7200,0.15,0.298,1.68;10000,0.25,0.396,2.3];

% Define the water conductivity
Kw=0.6;

% Define the depth of the stratigraphy from 0 to 10 km
z=[0:100:10000];

% Define the depth of the lithology change
Depth=[0;Tab(:,1)];

% Initialize the matrix of the conductvity
MagK1=zeros(100,1);
MagK2=zeros(100,1);

% Initialize the matrix of the thermal resistance mean
MagThermResMean1=zeros(5,1);
MagThermResMean2=zeros(5,1);

% Initialize the matrix of the harmonic mean
MagHarmM1=zeros(5,1);
MagHarmM2=zeros(5,1);

% Create a figure
figure

for i=1:5

NewDepth=(z(find(z==Depth(i)))+100:100:z(find(z==Depth(i+1))))./1000; 

% Estimate the porosity variation (function of depth)

fi=Tab(i,2).*exp(-Tab(i,3).*NewDepth);

% Estimate the thermal conductivity of each rock K1 (geometric mean) and K2 (square
% root mean)

K1=(Tab(i,4).^(1-fi)).*Kw.^fi;
K2=((1-fi).*(Tab(i,4)).^0.5+fi.*Kw^0.5).^2;

% plot the rocks' conductivities

plot(K1,NewDepth.*-1,'Marker','diamond','LineWidth',1)
hold on
plot(K2,NewDepth.*-1,'Marker','v','LineWidth',1)
hold on

% Create xlabel
xlabel('Thermal Conductivity (W/mK)','FontSize',20,'FontWeight','bold');

% Create ylabel
ylabel('Depth (km)','FontSize',20,'FontWeight','bold');

% Estimate the aritmetic mean of the thermal conductivity for each
% lithology

K1Mean=mean(K1);
K2Mean=mean(K2);

% Estimate the fracion of each lithology and of the thermal resistance

Fraction=(Depth(i+1)-Depth(i))./Depth(6);
ThermResMean1=Fraction./K1Mean;
ThermResMean2=Fraction./K2Mean;

MagThermResMean1(i)=ThermResMean1;
MagThermResMean2(i)=ThermResMean2;

% Estimate the fracion of the harmonic mean

HarmM1=((Depth(i+1)-Depth(i))/K1Mean);
MagHarmM1(i)=HarmM1;

HarmM2=((Depth(i+1)-Depth(i))/K2Mean);
MagHarmM2(i)=HarmM2;

% Write in two vectors the conductivities of the rocks

Z=find(MagK1==0);
S1=size(K1);

MagK1(Z(1):Z(1)-1+S1(2))=K1;
MagK2(Z(1):Z(1)-1+S1(2))=K2;

end

% Estimate the thermal resistance mean

ThermResMean1=sum(MagThermResMean1);
ThermResMean2=sum(MagThermResMean2);

% Estimate the harmonic mean

HarmonMean1=(Tab(5,1)./sum(MagHarmM1));
HarmonMean2=(Tab(5,1)./sum(MagHarmM2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ecercise 2

% Construct the geotherm with a variable gradient (from 30C/km to 21C/km at a depth of 10 km)

% Initialize the matrix and associate the first two values

T0=0;
T1=30;

MagT1=zeros(11,1);
MagT1(1)=T0;
MagT1(2)=T1;

% Construct the geotherm every km 

for h=1:9
   
T1=T1+(30-h);
MagT1(h+2)=T1;

end

% Interpolate the values of temperatures every 100 m (z)

Depthkm=[0:1000:10000];
NewTemp=interp1(Depthkm,MagT1,z);

% Estimate the thermal conductivity of the water for the two conditions

if NewTemp<=137
    Kw1=0.5648+1.878*10^-3.*NewTemp-7.231.*10^-6*(NewTemp).^2;
else
    Kw1=0.6020+1.309*10^-3.*NewTemp-5.140.*10^-6*(NewTemp).^2;
end

% Estimate the conductivity of the matrix

figure

for j=1:5

% Estimate the thermal conductivity of the rocks'matrix as function of
% temperature

Km=1.8418+(Tab(j,4)-1.8418).*(1./(0.002732.*NewTemp+0.7463)-0.2485);

% Define the depth of each rock (in m and km) with a step of 100 m

DepthRock=[z(find(z==Depth(j)))+100:100:z(find(z==Depth(j+1)))];
NewDepthRock=DepthRock./1000; 

% Estimate the size of the vector of the depth of each rock

SizeDR=size(DepthRock);
SizeDR1=SizeDR(2);

% Initialize the matrices of the thermal conductivity (with the size of the
% vector of the depth of each rock) and associate the values of the thermal
% conductivity of the rock matrix and water

MagKm1=zeros(SizeDR1,1);
MagKw1=zeros(SizeDR1,1);


for l=1:101
    PosDepth=find(z(l)==DepthRock);
    if PosDepth ~=0
    Km1=Km(PosDepth);
    MagKm1(PosDepth)=Km1;
    Kw2=Kw1(PosDepth);
    MagKw1(PosDepth)=Kw2;
    else
    end

% Estimate the porosity variation (function of depth)

fi=Tab(j,2).*exp(-Tab(j,3).*NewDepthRock');

% Estimate the thermal conductivity of each rock K1 (geometric mean) and K2 (square
% root mean)

Kmw1=(MagKm1.^(1-fi)).*MagKw1.^fi;
Kmw2=((1-fi).*(MagKm1).^0.5+fi.*MagKw1.^0.5).^2;

end

% plot the thermal conductivity (geometric mean and square root mean)

plot(Kmw1,NewDepthRock.*-1,'Marker','<','LineWidth',1)
hold on
plot(Kmw2,NewDepthRock.*-1,'Marker','>','LineWidth',1)
hold on

% Create xlabel
xlabel('Thermal Conductivity (W/mK)','FontSize',20,'FontWeight','bold');

% Create ylabel
ylabel('Depth (km)','FontSize',20,'FontWeight','bold');

end