%% Calibration of Pennink's (1905) model as run in First Series  (Sept 1904)

% This m-script uses fdm2.m and fdm2path.m to show the streamlines
% of Pennink's the sand-box model. This script determines the
% conductivity of the sand used. It turned out to be very highg, 865 m/d.
% Also, using travel times, the porosity was estimated at 38%,
% this porosity is in good agreement with the the model and the
% displacmement of the tracers (indian ink) applied by Pennink and measured
% from his photos.
%
% There is no need to do a complex parameter opbimizaion here, because the
% only unknown is the hydraulic conductivity. Hence we can just try some
% values using the given recharge and optimized the conductivity using the
% measured head difference between the east and west size of the model,
% given by Pennink as 1:87, which comes down to 65 cm/87=0.75 cm=7.5 mm.
% In his series 3 he refers to 1:87 and in his series 1 to the 6 cm
% (instead of 7.5), these differnces 1-1.5 mm can be regarded reading
% inaccuracies and therefore are ignored.
%
% TO 090314 100727

clear variables; close all;

MW=65;         % [cm] model width
MH=65;         % [cm] model height
MD=1.8;        % [cm] model thickness
M =37;         % [cm] x of center of recharge
W =45;         % [cm] width of applied rechargee
k=86500/24/60; % [cm/min] initical conductivity
por=0.38;      % effective porosity;
hCanL=45.2;    % [cm] foto p32
rch= 4200/60/MD/45; % [cm/min] 4200 cm3/60min/modlwidth/rechargewidth(45cm) = cm3/cm2/min=cm/min
fringe=2;      % [cm] capillary fringe thickness

% Contour of sand mass obtained form photo by ginput
xSand =[
   -0.0
   65.0
   65.0
   55.7
   21.8
   20.4
   13.8
    7.6
   -0.0
];

zSand =[
   -0.0
   -0.0
   59.0
   59.0
   57.6
   56.4
   50.8
   44.7
   44.2
   ];

xCanL =[
   -0.0
    5.5
    6.8
    8.4
    9.2
    9.2
   -0.0
   ];

zCanL =[
   37.0
   37.1
   43.1
   45.5
   59.4
   64.5
   64.5
   ];

xWatL =[
   -0.0
    5.5
    6.8
   -0.0
   ];

zWatL =[
   37.0
   37.1
   45.5
   45.5
   ];


xCanR=[
   65.0
   60.8
   61.6
   62.2
   62.7
   63.4
   65.0
];

zCanR =[
   65.0
   62.0
   51.2
   46.2
   43.3
   41.5
   41.3
];

xGr=0:1:65;
zGr=65:-1:0;

[xGr,zGr,xm,zm,Dx,Dz,Nx,Nz]=modelsize(xGr,zGr); [Xm,Zm]=meshgrid(xm,zm);

IBOUND=zeros(Nz,Nx);
IBOUND(inpolygon(Xm,Zm,zSand,zSand))=1;
IBOUND(inpolygon(Xm,Zm,xCanL,zCanL))=-2;
IBOUND(inpolygon(Xm,Zm,xCanR,zCanR))=-3;

HK=zeros(Nz,Nx).*k;
HK(inpolygon(Xm,Zm,xSand,zSand) & Zm<hCanL+fringe)=k;
HK(inpolygon(Xm,Zm,xCanL,zCanL) & Zm>hCanL+fringe)=0;

VK=HK;
VK(HK==0 & (Xm>M-0.5*W & Xm<M+0.5*W)) =k;

FH=NaN(Nz,Nx); FH(IBOUND==-2)=hCanL;

md=1; %md=thickness of numeical model cross section = 1 cm (everything is in cm and min)
FQ=zeros(Nz,Nx); I=find(xm>M-0.5*W & xm<M+0.5*W); FQ(1,I)=rch.*Dx(I)*md;

[Phi,Q,Psi,Qx,Qz]=fdm2(xGr,zGr,HK,VK,FH,FQ);    % 2D matlab finite difference model

figure; hold on
patch([0 MW MW 0],[0 0 hCanL hCanL],'b');
patch(xSand,zSand,'y');
patch(xCanL,zCanL,'w');
patch(xWatL,zWatL,'b');
contour(xm        ,zm, Phi,20,'r');                  % 20 head contour lines
contour(xGr(2:end-1),zGr, Psi,20,'k');               % 20 stream lines

%% Compute flow for calibration

ObsPnt=[64.5, hCanL-0.5];
[ix,iz]=xyzindex(ObsPnt(1),ObsPnt(2),xGr,zGr);

Qout=sum(Q(IBOUND<0))*MD;
dh=Phi(iz,ix)-hCanL;

fprintf('We expect   Qout=                      -4.2 L/h given dh=65/87=0.75 cm\n');
fprintf('We obtained Qout=%8.1f cm3/min =%8.1f L/h; and  dh=    %6.2f cm.\n',Qout,Qout*60/1000,dh)
fprintf('This confirms the calibration of HK\n');
% 

xlabel('x cm'); ylabel('z cm');
title(sprintf('Calibration run Pennink Series 1, Q=%.1f L/h and dh=%.1f cm, k=%.0f m/d\n',Qout*60/1000,dh,k/100*60*24));


%% travel times
% This section compute travel times using particle tracking. It applies the
% matlab functin fdm2path.m in the mfiles directory. This is a function
% which was written entirely in Matlab, but it uses the same methods as
% Modpath. It is used here only to obtain an impression of the travel times
% invoved in this test.
% We first define the two points where Pennink released his indian ink.

startpoint(1).x=53.3; % [cm]
startpoint(1).z=56.5; % [cm]

startpoint(2).x=38.3; % [cm]
startpoint(2).z=56.5; % [cm]

% The for convenience of plotting, digitized a bit crudely the plumes form
% Penninks photo using the function MeasureOnPhoto in the Penninks/Photo
% directory: The times are minuts since the beginning of the ink release as
% deduced from the times the photos were taken. Times in minutes.

startpoint(1).plume(1).x=[54.3 59.3 59.3 57.6 53.7 49.8 43.1 38.1 31.9 27.1 22.8 22.6 31.9 37.7 43.3 47.6 49.6 50.5 50.9 51.3];
startpoint(1).plume(1).z=[58.4 40.6 34.6 28.9 23.2 21.1 18.5 17.4 17.2 16.9 19.0 20.3 21.6 25.5 30.7 38.6 44.8 51.3 56.3 59.2];
startpoint(1).plume(1).t=57; % 9:22-10:19

startpoint(1).plume(2).x=[55.5 57.2 58.3 58.7 58.5 56.6 52.9 47.4 39.3 30.5 21.1 14.7 6.5  4.2  6.3  12.5 23.4 29.6 35.4 41.4 46.5 49.1 50.4 51.4];
startpoint(1).plume(2).z=[59.6 52.7 45.7 38.8 31.9 26.8 19.9 16.7 13.5 12.2 14.8 18.3 24.4 28.4 26.6 22.8 20.7 20.7 23.6 29.5 36.9 45.5 53.2 59.6];
startpoint(1).plume(2).t=77; % 9:22-10:39

startpoint(1).plume(3).x=[56.2 59.5 59.9 58.6 55.5 51.4 42.9 33.7 23.7 13.4  8.2  4.9  1.7  1.4  1.0  3.0  9.1 16.5 25.2 33.7 40.7 45.3 48.8 51.0 51.2];
startpoint(1).plume(3).z=[58.6 44.0 35.3 28.7 23.7 18.4 12.6  9.7  9.7 12.6 17.9 22.1 29.8 32.4 36.1 31.4 24.5 20.8 20.5 23.7 28.5 35.6 43.5 50.4 58.8];
startpoint(1).plume(3).t=107; % 9:22-11:09

startpoint(1).plume(4).x=[55.9 59.0 59.5 58.0 51.3 41.9 30.5 19.6 11.7  6.0  2.2  0.3 -0.0  5.0 12.9 23.6 34.1 42.9 47.3 49.8 51.1 51.3];
startpoint(1).plume(4).z=[58.4 44.2 30.9 26.9 17.9 10.7  7.7  7.7 10.9 16.2 21.2 26.2 36.9 27.2 21.2 19.9 22.7 28.4 37.7 46.4 54.7 59.2];
startpoint(1).plume(4).t=137; % 9:22-11:39

% plume from second start point
startpoint(2).plume(1).x=[40.4 40.4 39.8 35.1 28.9 21.3 14.8 10.4  6.6  3.9  3.5  5.2  9.6 13.8 19.6 24.7 29.7 31.8 32.6 33.9];
startpoint(2).plume(1).z=[58.4 49.2 42.4 37.7 34.4 32.2 31.9 32.9 35.7 38.4 39.4 38.4 36.4 36.2 36.2 37.4 42.2 47.4 52.4 58.4];
startpoint(2).plume(1).t=27; % 11:12-11:39


%% Times to obtain particla position (these match the plume times. The
% times of StartPoint(2) have also intermedite times, otherwise the partile
% may have already been discharge at the time of the first photo after release.
% This is valid for the second (more left) release point.

startpoint(1).times=[0, 27, 57, 77, 107, 137]; % minutes
startpoint(2).times=[0:5:27 27]; % minutes
markers='o*^sp';

%% Plot patches to show the plumes (one on top of the other)
for iPnt=1:length(startpoint)
    for i=length(startpoint(iPnt).plume):-1:1
        patch(startpoint(iPnt).plume(i).x,startpoint(iPnt).plume(i).z,'r','facealpha',0.2*i/2);
    end
end

%% show particle path with markers at specified times
for iPnt=1:length(startpoint)
    [XP ZP TP]=fdm2path(xGr,zGr,1,Q,Qx,Qz,por,startpoint(iPnt).times,markers,startpoint(iPnt).x,startpoint(iPnt).z);
end

% REMARKS
% 1) The stream lines sharply bend at the top of the saturated zone
% 2) The downward vertical streamlines represent recharge
% 3) The upward vertical streamlines at the extraction canal are OK, a bit
% ugly but they are the flow out of the model. The lines are the so-cut
% branchcut of the stream lines due to sink or sources inside the medium,
% i.e. at the locations where water is removed from the model at the
% circumference of the canal. These streamlines can be exlcuded from
% printing by properly masking the Psi function with NaNs outside the model
% space. This is left to the user as an exercise. (Note that I did so in
% the examples Pennink/Series6.
