%% Pennink1 - density flow in Pennnink's (1915) sand box model
% Experiments series 3
%
%see http://www.citg.tudelft.nl/live/pagina.jsp?id=68e12562-a4d2-489a-b82e-deca5dd32c42&lang=en
%
% in this experiment, Pennink studies simultaneous freshwater and saltwater
% flow in a cross section fed by rain and discharging to a canal at the
% westsided of the section.
% He first establiches a saltwater zone by letting milk enter slowly at the
% bottom of the model, regulated by a fixed-head reservoir.
% After the milk volum as reached 13 cm height in the model, he starts the
% precipitation 4.2 L/h over a width of 45 cm.
% This causes the interface to reach equilibrium.
% During this phase, milk is injected continuously at a rate of 2.14
% ml/min, the rate that fills the model to 13 cm in 5 hours time, as was
% the time Pennink took to let the milk enter the model. This rate seems
% also ok in the subsequent simulation. Pennink mentions that the milk was
% added drop by drop. THerefore, no fixed saltwater head was applied, but a
% well was used as the boundary condition for the milk entry point instead.
% After a day, Pennink stops adding milk and then after one more day, he
% takes the last photo of this test series.
% As it takes about two hours to establish virtual equilibrium of the
% saltwater interface, I shortened the simulation to 2 hours for the time
% to reach equilbrium with the milk injection and then continu during a
% full day without the milk injecition to simulate the washing out of the
% milk from the model. This should allow for calibration of the
% dispersivity in Penninks original model. This calibration was not really
% attempted as the situation after a full day without milk addition was
% similar enough the that of the last photo in Pennink's series 3 test.
% A somewhat better result may be obtained by raising the dispersivity a
% bit.
% Computation time about 30 min. To reduce the compuation time you may use
% a coarser grid, say 2 cm instead of 1 cm.
%
% As alway, the models can be readily run with mt3dms as with seawat by
% switching the appropriate packages in the NAM worksheet. An alternative
% to running mt3dms is to run SEAWAT without density active. This is done
% by setting -1 for the on switch for the VDF package in the NAM worksheet
% instead of 1.
%
% TO 090312 100523 100721

clear variables; close all;

basename='PenninkSeries3';
save name basename

fprintf('Tests Pennink (1915), series 3\n');

xSand =[
   -2.2878
   67.6080
   67.7985
   43.8015
   30.0890
   23.2327
   17.7096
   11.4247
   10.4725
    8.5679
   -2.2878
   ];

zSand =[
   -2.3540
   -2.5445
   60.3133
   60.1229
   59.1705
   57.8371
   55.9323
   52.6942
   52.6942
   50.4085
   49.2656
   ];

xCanL =[
    -2.2878
    3.2353
    5.7112
    7.4252
    7.8061
   10.4725
   -2.2878
   ];

zCanL =[
   39.7417
   38.7893
   41.0750
   45.2656
   50.7894
   65.8372
   65.6467
   ];

xCanR=[
   68.3698
   66.2748
   61.1326
   61.3231
   61.8944
   63.0372
   65.5130
   69.1316
   68.7507
   ];

zCanR =[
   63.3610
   66.2182
   65.2658
   51.7418
   46.9799
   42.4084
   41.2655
   41.4560
   64.3134
];


%% The grid, the box is 96x96 cm and 1.8 cm thick
FRESHWATER=0;    % Relative minimum concentration
SEAWATER  =1;    % Relative maximum concentration
k=86500/(24*60); % cm/min calibrated
peff=0.38;       % calibrated

MW   =65;         % [cm] Width of model. Pennink p6
MH   =65;         % [cm]
D    = 2;         % [cm] thickness of model
Fringe=2.0;       % [m] Capillary fringe see page 34
zIface=13;        % [cm]

xzMilkSupply=[63.5,0.1]; % point of milksupply

%% Grid is always 3D
dx=2.0;              % [cm] grid cell width
dy=  D;              % [cm] grid cell length = depth of model
dz=2.0;              % [cm] grid cell size vertical

xGr=0:dx:MW;         % [cm] grid
yGr=[0 D];           % [cm] grid
zGr=(MH:-dz:0);      % [cm] grid

[xGr,yGr,zGr,xm,ym,zm,DX,DY,DZ,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr);

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

[ixMilk,izMilk]=xyzindex(xzMilkSupply(1),xzMilkSupply(2),xGr,zGr);

Z=zeros(Ny,Nx,Nz);
for i=1:length(zGr)
    Z(:,:,i)=zGr(i);
end

%% Model arrays
hCanL=45; hCanR=46; % approximate canal water elevation (see workbook sheet PER)

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;  % not a fixed head
IBOUND(IBOUND==2 & ZM>hCanL)=0;
IBOUND(IBOUND==3 & ZM>hCanR)=0;
IBOUND(izMilk,ixMilk)=5;         % milk injection point

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

PEFF=ones(size(IBOUND))*peff;
PEFF(ZM>hCanL+Fringe)=peff/3;           % unsaturated
PEFF(IBOUND==0)=0;             % inactive

% 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]);

ICBUND=ones(size(IBOUND));

STRTHD=ones(size(IBOUND))*0.5*(hCanL+hCanR);

STCONC=ones(Ny,Nx,Nz)*FRESHWATER;
STCONC(:,:,zm<zIface)=SEAWATER;  % zm < interface as read from photo in mm above base

[PerNams,PerVals]=getPeriods(basename);  % get stress periods from workbook
NPER = size(PerVals,1);                     % total number of stress periods
%% Get period data to set CHD boundaries for MODFLOW and PNTSRC for MTRDMS/SSM

[PerNams,PerVals]=getPeriods(basename);  % get stress periods from workbook

NPER=size(PerVals,1);
hL   =PerVals(:,strmatchi('hL',   PerNams));
QMilk=PerVals(:,strmatchi('QMilk',PerNams));

save underneath xSand zSand

IL=cellIndices(find(IBOUND==2),size(IBOUND),'LRC'); uL=ones(size(IL(:,1)));
IS=cellIndices(find(IBOUND==5),size(IBOUND),'LRC'); uS=ones(size(IS(:,1)));

ITYPE=2; % well

iPer=1;
CHD=[ iPer*uL IL uL*[hL(iPer) hL(iPer)]]; % Left   canal
WEL=[ iPer*uS IS uS*QMilk(iPer) ]; % Milk injection point

PNTSRC=[iPer*uS IS uS*[SEAWATER ITYPE SEAWATER]]; % milk injection point

for iPer=2:NPER
    if   hL(iPer) == hL(iPer-1), CHD=[CHD; -iPer NaN(1,5)];
    elseif hL(iPer)>0,           CHD=[CHD;  iPer*uL IL uL*[hL(iPer) hL(iPer)]];
    end
    
    if  QMilk(iPer)==QMilk(iPer-1), WEL=[WEL; -iPer NaN(1,4)];
    elseif abs(QMilk(iPer))>0,      WEL=[WEL;  iPer*uS IS uS*QMilk(iPer)];
    end
    
    PNTSRC=[PNTSRC; iPer*uS IS uS*[SEAWATER ITYPE SEAWATER]];
end

%% RCH Recharge see below N=4.2 L/h according to Pennink

W=45; % [cm] width of rain added to top of model
M=37; % [cm] center of rain added to top of model (mid between 2 screws see photo)
N=4.2*1e3/60/W/dy;  % recharge in cm/min: 4.2 L/h / 60 min/h / W / dy

RECH=zeros(Ny,Nx,NPER);
RECH(:,xm>=M-0.5*W & xm<=M+0.5*W,:)=N;

save underneath xSand zSand xCanL zCanL % needed to set PSIMASK
