%% Pennink1 - density flow in Pennnink's (1905) sand box model
% Experiments series 2 (Septemner 1904)
%
%see http://www.citg.tudelft.nl/live/pagina.jsp?id=68e12562-a4d2-489a-b82e-deca5dd32c42&lang=en
%
% in this experiment, Pennink studies freshwater flow between a recharge
% canal on the right of the model to an extraction canal at the left.
% He then injects ink at different depths near the recharge canal and
% shows in a series of photos the movement of the ink.
%
% The flow is from the right to the left canal and the flow is
% 1.25 L/h (Pennink (1905) p485 with a (water table) gradient of 1:130
%
%
% This model can be run by mt3dms as well as by seawat.
% To easier match the actual times of the tests, we will use time in minutes
% instead of days.
%
% You may want to try VDF 1 instead of -1 in the NAM worksheet
% which cause the Ink to be of seawater density
% You may also swicth off VDF in that sheet and run mt3dms to do this
% switch off VDF adn GCG and switch on FTL and LMT in the NAM worksheet
%
% TO 090312

clear variables;
close all;

basename='Pennink_Series2';
save name basename

%% Parameters

por=0.38;      %   calibrated
k=86500/24/60; % [cm/min] Conductivity calbrated from data in Pennink's series 1 experiments 
FRESHWATER=0;  % Freshwater concentration
SEAWATER  =1;  % not used in this simulation

%% Times of Photos taken as printed in the book Pennink(1915)

Times=[1904,9,1, 9,0;  % start of experiment
    1904,9,1, 9,40;    % ink added to injection points
    1904,9,1, 9,55
    1904,9,1,10,20
    1904,9,1,10,55
    1904,9,1,11,25
    1904,9,1,13,20     % addition of ink stopped at this time
    1904,9,1,14,15
    1904,9,2, 7,50     % last photo
];

% %contour of sand mass (obtained with ginput form photo)
xSand =[
   -3.3152
   68.5366
   69.4442
   41.4598
   31.4762
   24.8205
   19.6774
   14.6856
    9.9963
    7.7273
   -3.3152
];

zSand =[
   -2.9318
   -2.4127
   61.2618
   59.5316
   59.1855
   58.3204
   56.7631
   55.3789
   52.6104
   50.3610
   49.3229
];

xCanL =[
   -2.9049
    4.3070
    6.4706
    7.9129
    9.1750
   -5.4291
];

zCanL =[
   39.6954
   39.5234
   43.9957
   51.2203
   62.9171
   65.4973
];

xCanR=[
   71.1975
   60.3797
   61.1009
   61.6418
   63.0841
   64.1659
   66.6901
];

zCanR =[
   67.5615
   66.1854
   53.9725
   47.6080
   43.1357
   42.7916
   43.3077
];

%% The grid, the box is 96x96 cm and the sand-box model is 1.8 cm thick

MW   =65.0;       % [cm] Width of sand-box model, see Pennink p6
MH  =65.0;        % [cm] Top of model
D   = 2.0;        % [cm] Thickness of model
zCap=51.0;        % [cm] Top of full capillary zone (see descripition)
hCanL=45.2;       % [cm] see photo p32
hCanR=45.7;       % [cm] Stage in right-hand canal, from given gradient and left canal stage

%% Point injection locations
xzInk=[
   50.6   50.0
   52.6   44.9
   53.8   33.7
   56.1   23.0
];

%% Grid everything in m
dx=1.0;      % [cm] model grid cell width
dy=  D;      % [cm] model grid cell in row direciton, taken model as 1 m deep for easier water budget computation
dz=1.0;      % [cm] model grid cell heigth

xGr=0:dx:MW;     % mesh grid lines
yGr=[0 D];       % mesh grid lines, however the real model was 1.8 cm thick
zGr=(MH:-dz:0)'; % mesh grid lines

[xGr,yGr,Z,xm,zm,zm,Dx,Dy,Dz,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr);  % mesh houskeeping

[XM,ZM]=meshgrid(xm,zm);

[ixInk,izInk]=xyzindex(xzInk(:,1),xzInk(:,2),xGr,zGr);

%% Model arrays

IBOUND=zeros(Nz,Nx);
IBOUND(inpolygon(XM,ZM,xSand,zSand))=1;
IBOUND(inpolygon(XM,ZM,xCanL,zCanL))=-2;
IBOUND(inpolygon(XM,ZM,xCanR,zCanR))=-3;
IBOUND(IBOUND==-2 & ZM>hCanL)=0;
IBOUND(IBOUND==-3 & ZM>hCanR)=0;

HK = ones(size(IBOUND))*k; HK(IBOUND==0)=0;
VK =HK;
HK(HK>0 & ZM>zCap)=k/10; % HK=0 above full capillary zone

PEFF=ones(size(IBOUND))*por;
PEFF(ZM>zCap)=por/3;           % unsaturated

STRTHD=ones(size(IBOUND))*0.5*(hCanL+hCanR);
STRTHD(IBOUND==-2)=hCanL;  % gegeven peil links fixed head
STRTHD(IBOUND==-3)=hCanR;  % gegeven peil rechts fixed head
STCONC=ones(size(IBOUND))*FRESHWATER;
ICBUND=ones(size(IBOUND));

%% Permute Arrays to 3D
IBOUND=permute(IBOUND,[3,2,1]);
HK    =permute(HK    ,[3,2,1]);
VK    =permute(VK    ,[3,2,1]);
PEFF  =permute(PEFF  ,[3,2,1]);
STRTHD=permute(STRTHD,[3,2,1]);
STCONC=permute(STCONC,[3,2,1]);
ICBUND=permute(ICBUND,[3,2,1]);

%% Generate the source-boundary conditions for the ink injection
%  taking the injection points as constant concentration cells
ITYPE=-1;  % fixed concentration source

[PerNams,PerVals]=getPeriods([basename,'.xls']);
C_Ink=PerVals(:,strmatchi('C_Ink',PerNams));

NPER=size(PerVals,1);
u    =ones(size(ixInk(:,1)));

PNTSRC=[];
for iPer=1:NPER
    PNTSRC=[PNTSRC; [u*iPer izInk u ixInk u*[C_Ink(iPer) ITYPE C_Ink(iPer)]]];
end

save underneath xSand zSand xCanL zCanL xCanR zCanR % needed in mf_analyze to set PSIMASK for plotting
        