%% Pennink Series 6 - Density flow simulation of Pennnink's (1915) sand-box model tests
% -- SIXTH SERIES OF EXPERIMENTS (MARCH 1905)
%
%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 from a canal in the right and one in the left corner of the sand-box model
% to a canal in the center. The following text is cited from his book, p87:
%
% "C. Experiments like sub B, but central-canal as a means of concentration
% substituted by a well-pipe with gauze strainer, placed in the canal axis.
%  Sixth series of experiments (March 1905)
% The water is discharged by means of the well-pipe, by siphoning, on the
% left side of the apparatus. Moreover, glass gauging-pipes are place in
% the sand in order to determine the pressure of the water just above the
% milk, and beside the gause strainer."
%
% The exact construction of the "well-pipe" is not discribed by Pennink, neither  
% nor its depth and screen length. In the simulatons done here, a fully screened
% well is assumed to a depth decuced from the photos in the book.
%
% TO 090320 100520

clear variables; close all;

basename='PenninkSeries6'; fprintf('Model %s\n',basename);

%% Parameters
FRESHWATER=0;           % Relative conc of freshwater
SEAWATER  =1;           % Relative conc of saltwater (milk)
peff      =0.38;        % [-] calibrated effective porosity
k   =86500/(24*60);     % [cm/min] calibrated conductivity

% All coordinates relative to the LL corner of the sand mass in the model
MH      =96;    % [cm] Pennink specified this on p63
MW      =97;    % [cm] Pennink specified this on p63
D       =1.8;   % [cm] Thickness of model sand mass
hCanal0 =91;    % [cm] measured from photo p88
zIFace  =23;    % [cm]
zScrBot =zIFace+32; % [cm] 32 cm above original milk surface

% countor around sand mass
xSand =[
   -0.0
   97.0
   97.0
   48.5
   -0.0
];

zSand =[
   -0.0
   -0.0
   92.1
   90.7
   92.1
];

% saturated sand
xSat =[
   -0.0
   97.0
   97.0
   48.5
   -0.0
];

zSat =[
   -0.0
   -0.0
   88.0
   85.0
   88.0
];

% contour around left canal
xCanL =[
   -0.0
    0.9
    1.5
    2.0
    2.5
   -0.0
];


zCanL =[
   81.0
   81.0
   81.8
   88.0
   96.0
   96.0
];

xWatL =[
   -0.0
    0.9
    1.5
    2.0
   -0.0
];

zWatL =[
   81.0
   81.0
   81.8
   88.0
   88.0
];

xCanR=MW-xCanL;
zCanR =zCanL;

xWatR=MW-xWatL;
zWatR=zWatL;

% contour around center (mid) canal
xCanM=MW/2+[
   -3.5
   -3.0
   -2.2
   -0.9
   -0.0
    0.9
    2.2
    3.0
    3.5
];

zCanM =[
   96.0
   88.0
   80.0
   75.3  
   74.5
   75.3
   80.0
   88.0
   96.0
];

xWatM=MW/2+[
   -2.2
   -0.9
   -0.0
    0.9
    2.2
];

zWatM =[
   80.0
   75.3  
   74.5
   75.3
   80.0
];

xzMilk=[MW-1.5, 0.5];  % coordates of point where milk was supplied

%% Model grid
dx=1.0;       % [cm] grid horizontal cell size
dy=1.8;       % [cm] thickness of model as specified by Pennink
dz=1.0;       % [cm] grid vertical cell size

xGr=0:dx:MW;    % xGr-grid coordinates
yGr=[0 dy];     % one row of width of the model
zGr=MH:-dz:0;   % zGr line elevations

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

Z=zGr;

[XM,ZM]=meshgrid(xm,zm); % Full grid of cell centers

[ixM,iyM,izM]=xyzindex(xzMilk(1),D/2,xzMilk(2),xGr,yGr,zGr);

IBOUND=zeros(Nz,Nx);
%IBOUND(inpolygon(XM,ZM,xSand,zSand))=1;  % all sand active
IBOUND(inpolygon(XM,ZM,xSat ,zSat ))=1;  % only saturated sand active
IBOUND(inpolygon(XM,ZM,xCanL,zCanL) & inpolygon(XM,ZM,xSat,zSat))=2;
IBOUND(inpolygon(XM,ZM,xCanM,zCanM) & inpolygon(XM,ZM,xSat,zSat))=3;
IBOUND(inpolygon(XM,ZM,xCanR,zCanR) & inpolygon(XM,ZM,xSat,zSat))=4;
IBOUND(izM,ixM)=5; % milk supply point

HK=zeros(size(IBOUND));
HK(inpolygon(XM,ZM,xSand,zSand))=k/10;  % unsaturated all sand
HK(inpolygon(XM,ZM,xSat ,zSat ))=k;     % saturated   saturated sand
VK  = HK;

PEFF= ones(size(IBOUND))*peff;
STRTHD=ones(size(IBOUND))*hCanal0;

ICBUND=IBOUND;

STCONC=ones(size(IBOUND))*FRESHWATER; % Start with all cells no seawater
STCONC(ZM<zIFace)        =SEAWATER;   % Make cells below interface salt

% Turn 2D matrices into 3D arrays
IBOUND=permute(IBOUND,[3,2,1]);
ICBUND=permute(ICBUND,[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]);

%%  NOTICE
% We only model the situaton after Pennink pulled up his well pipe to the
% final position of 32 cm above the original interface.
%
% erroneous value mau cause the milk to fast in or outflow the model.
% Measuring directly from the photo was difficult and computing nog accurate
% enough because there are vertical gradients in the sand above the interface.
% Hence it was computed first and then fine-tuned in several trial runs.
% The applied head values are in the workbook on sheet PER.

%% Find the cells of the center well using figures given by Pennink
zW=interp1([zScrBot,MH],[zScrBot,MH],zm);
zW(isnan(zW))=[];
xW=ones(size(zW))*MW/2;
yW=ones(size(zW))*D/2;
[ixW,iyW,izW]=xyzindex(xW,yW,zW,xGr,yGr,zGr(:));

xWel=xGr([ixW(1) ixW(1) ixW(1)+1 ixW(1)+1]);  
zWel=zGr([izW(end) izW(1) izW(1) izW(end)]); zWel=zWel(:);

GHBCOND=1e8;

%% Setting boundaries for the stress periods
[PerNams,PerVals]=getperiods(basename);
hL=PerVals(:,strmatchi('hL',PerNams));
hM=PerVals(:,strmatchi('hM',PerNams)); % well
hR=PerVals(:,strmatchi('hR',PerNams));
hS=PerVals(:,strmatchi('hS',PerNams));

IL=cellIndices(find(IBOUND==2),size(IBOUND),'LRC'); uL=ones(size(IL(:,1)));
IM=cellIndices(find(IBOUND==3),size(IBOUND),'LRC'); uM=ones(size(IM(:,1)));
IR=cellIndices(find(IBOUND==4),size(IBOUND),'LRC'); uR=ones(size(IR(:,1)));
IW=[izW,iyW,ixW];                                   uW=ones(size(IW(:,1)));
IS=[izM,iyM,ixM];                                   uS=ones(size(IS(:,1)));

GHB=[]; PNTSRC=[]; ITYPE=5;
for iPer=1:size(PerVals(:,1));
    if iPer>1 && ...
         all([hL(iPer) hM(iPer) hR(iPer) hS(iPer)]==...
             [hL(iPer-1) hM(iPer-1) hR(iPer-1) hS(iPer-1)])
       GHB=[GHB; -iPer NaN(1,5)];
    else
        if ~isnan(hL(iPer)), GHB=[GHB; iPer*uL IL uL*[hL(iPer) GHBCOND]]; end
        if ~isnan(hR(iPer)), GHB=[GHB; iPer*uR IR uR*[hR(iPer) GHBCOND]]; end
      % if ~isnan(hM(iPer)), GHB=[GHB; iPer*uM IM uM*[hM(iPer) GHBCOND]]; end
        if ~isnan(hM(iPer)), GHB=[GHB; iPer*uW IW uW*[hM(iPer) GHBCOND]]; end
        if ~isnan(hS(iPer)), GHB=[GHB; iPer*uS IS uS*[hS(iPer) GHBCOND]]; end
    end
    if ~isnan(hS(iPer)), PNTSRC=[PNTSRC; iPer*uS IS SEAWATER ITYPE]; end
end

save underneath xSand zSand xSat zSat xCanL zCanL xCanR zCanR xCanM zCanM xWatL zWatL xWatR zWatR xWatM zWatM xWel zWel IL IR IM IW IS
