Flow in porous media 195
% Calculation of new depth values of the lower layer;
for i=2:nx-1
Qwu=0;
if i==iwu
Qwu=qu;
end
qr=-K*(ho(i+1)+hu(i+1)-ho(i)-hu(i)-rdf*(ho(i+1)ho(i)))/Dx*(hu(i)+hu(i+1))/2;
ql=-K*(ho(i)+hu(i)-ho(i-1)-hu(i-1)-rdf*(ho(i)ho(i-1)))/Dx*(hu(i)+hu(i-1))/2;
hun(i)=hu(i)-Dt*(qr-ql-Qwu)/Dx/p;
if hun(i)<0
hun(i)=0;
end
end
% Boundary conditions;
hun(nx)=hun(nx-1);
% Updating depth values;
for i=1:nx
ho(i)=hon(i);
hu(i)=hun(i);
end
chk=0;
xtoe=0;
% Estimation of saline wedge toe location;
for i=1:nx
if hu(i)<0.1
xtoe=i;
break
end
end
toe(k)=xtoe;
index=k
end
i=1:nx;
GWt=ho(i);
SWe=hu(i);
ground = GWt+SWe;
plot(1:k,toe,'b','Linewidth',1.5)
xlabel('Time steps x 0.1 [days]')
ylabel('Position of the saline wedge toe x 200 [meters]')
figure
plotyy(1:nx,SWe,1:nx,ground)
xlabel('Distance from the coast x 200 [meters]')
ylabel('Interface [meters]')
text(3,15,'Saline wedge');
text(25,8.8,'Fresh water table');