%% Pennink - density flow in Pennnink's (1915) sand box model
% -- FOURTH 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 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, p63:
%
% "B. Phenomena of movement of the liquids of different specific gravity,
% viz water and milk in case the discharge canal is placed in the centre
% and the water symmetrically supplied on both sides"
% Height of the apparatus   = 0.96 m
% Breath of the apparatus   = 0.97 m
% Thickness of the sandmass = 0.018 m
% The reservoir which can be raised and lowered is applied on the right
% side, it was filled with a heavier liquid, viz milk (speicif gravity =
% +/- 1.03), which can enter ton the right hand side from the bottom.
% Water is supplied in the upper corners. Quantity of water = 28 liter
% hourly."
%
% But the surface of the sand has a v-shape with the lowest point in the
% center. THis shape can be used to fix a head gradient at the top of the
% model.
% Initially the water level is above the sand surface and, therefore,
% horizontal, with no flow in the model. The initial level of the milk,
% used as a surrogate for saltwaer, is 25 cm above the bottom of the model.
% The milk is added via a tube near the lower right bottom of the model
% using a reservoir kept at a suitable % elevation to create the initial
% density distribution.
% Next, the water flow is initiated and maintained by adding water to the
% left and right edges on top of the sand surface and extraction it at the sand
% surface in the center of the model. This creates a head gradient equal to
% the inclination of the sand surface, 1:16 according to Pennink on P64.
% Therefore, the flow is driven by the head, which makes the used water
% flow to useless information. It was 28 L/h as mentioned by Pennink.
%
% Because there is no flexibility in this model, there is only one photo in
% Penninks book, showing the situation after some time, which can be
% regarded as the final equilibrium.

% TO 090312 100523 100719

clear variables; close all;

basename='PenninkSeries4';
save name basename

fprintf('Tests Pennink (1915), Test series 4, 17 March 1905\n');

%% 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;       % [-] effective porosity calibrated
MW   =97;        % [cm] Width of model. Pennink p6
MH   =96;        % [cm] Height of model
D   =1.8;        % [cm] thickness of sand in model

%% Contour of sand body, all coords relative to LL of model sand
xSand =[
   -3.1383
   97.0917
   97.3429
   48.1071
   -0.1238
   -3.1383
];

zSand =[
   -3.4174
   -4.8567
   93.8822
   90.4278
   94.7459
   95.6095
   ];

grad=1./16;  % gradient of sand top acc to Pennink

% top of sand surface, gradient 1/16 acc to Pennink
xTop=   [ 0   0.5 1   ]*MW;
zTop=91+[ 0.5 0   0.5 ]*MW/16;  % 91 cm is top of sand in center of model
    
zIface = 26; % initial elevation of interface

xzMilkSupply=[91.5, 1.5];  % coordates of point where milk was supplied
ITYPE=1;                   % well source used for inflow from CHD nodes

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

xGr=0:dx:MW;         % [m] grid
yGr=[0 dy];          % [m] grid
zGr=(MH:-dz:0)';     % [m] grid

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

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

[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

IBOUND=zeros(Nz,Nx);
IBOUND(inpolygon(XM,ZM,xSand,zSand))= 1;
IBOUND([IBOUND(1,:); diff(IBOUND,1,1)]==1)=-1;
IBOUND(izMilk,ixMilk)=5;  %milk injection point

HK = ones(size(IBOUND))*k;
HK(IBOUND==0)=0;             % this is redundant as cells are inactive
VK=HK;

PEFF=ones(size(IBOUND))*peff;   % effective porosity
PEFF(IBOUND==0)=0;              % inactive, redundant

STRTHD=ones(size(zm))*interp1(xTop,zTop,xm);

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

ICBUND=ones(size(IBOUND));

STCONC=ones(Ny,Nx,Nz)      * FRESHWATER;
STCONC(:,:,zm<zIface)=SEAWATER;  % zm < interface  --> SEAWATER

%% Get period data to set CHD boundaries for MODFLOW and PNTSRC for MTRDMS/SSM
% This is only desired here for the milk injction point which has a fixed
% saltwater head to suppy or extrat milk automatically from and to the
% milk reservoir

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

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

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

ITYPE=1; % constant head cell (milk injection point)

iPer=1;
CHD=   [iPer*uS IS uS*[hS(iPer) hS(iPer)]];    
PNTSRC=[iPer*uS IS uS*[SEAWATER ITYPE SEAWATER] % milk injection point
];

for iPer=2:NPER
    if   hS(iPer) == hS(iPer-1)
         CHD   =[CHD   ; -iPer NaN(1,5)];
    else
        CHD=[CHD;
            iPer*uS IS uS*[hS(iPer) hS(iPer)]
        ];
    end
    PNTSRC=[PNTSRC; iPer*uS IS uS*[SEAWATER ITYPE SEAWATER]];

end


save underneath xTop zTop
