% MF_ADAPT -- Hydrocoin problem, example comes with Seawat
% Guo and Langevin (2002), p72ff
%  TO 090101

basename='Hydrocoin';

% The pupose of the Hydrologic Code Intercoparision (HYDROCOIN) project was
% to evaluate the accuracy of selected ground-water modeling codes.
% One of the porblems used for testing is the HYDROCOIN prolem presented
% here. It concerns a vertical cross section of 900 m long and 300 m
% height, has a sloping pressure boundary at the  top and no flow at the
% bottom and the sides. A constant concentration condition is applied on
% the middel 300 m of the bottom, to represent the top of a salt dome. Salt
% dispersis from the dome and is advected with the ambient groundwater
% flow.
% The applied mesh is 76 layer and 45 column by 1 row with cells of 20 m
% in xGr-direction, 4 in the vertical direction and 1 m along the y-direction.
% The constant head boundary along the top varies form 10.080 m at the
% center of the first column to 0.113 m at the center of the last.
% Fluid entering the model from the upper boundary has zero concentration
% and density 1000 kg/m3. The first and last 15 columns along of the bottom
% row are inactive. Between these inactive cells a constant concentration
% of 280 kg/m3 is maintainted.
% The uniform freshwater hydraulic conductivity is 0.8476 m/d. Only in the
% bottom layes was this conductivity set to 0.8475e-4 m/d to limit the
% advection from this layer.
% The problem has been published by Konikow a.a. (1977)
% The density of the Cbrine is 1200 kg/m3 at a concentration of 285 kg/m3
% salinity. Hence the drhodc=200/285 [dimensionless].
% We could just as well use relative concentration directly in SEAWAT or
% compute relative concentrations form the output by devision of the
% resulting concentrations by 285 kg/m3.
% The simulation time is 100 years in order to reach equilbrium. This will
% be computed in 10 stress periods state with 1 "time" step of increasing length
% (see worksheet for this problem) to get output at the end of each time
% step, whle a large number NPRS makes sure that also MT3D ouputs only at
% the end of each time step.
% TO 091208


%% Generate a Henry benchmark problem see Guo and Langevin different reports

peff=0.2;     % [  ] effective porosiy
k= 0.8476;    % m/d freshwater hydraulic conductivity
Cfresh =0;    % ;kg/m3] salinity of fresh water
Cbrine =1;  % [kg/m3] relative concentration
rhob  =1200;  % [kg/m3] brine density
rhof  =1000;  % [kg/m3] freshwater densit
aL    = 20;   % [m] longitudinal dispersivity
aT    =  2;   % [m] transverse dispersivity
Dm    =  0;   % [m2/d] molecular diffusion coefficient
strthd =0;    % [m] initial head
Ztop =   0;   % [m] to of model
Zbot =Ztop-300;% [m] bottom of model
xLeft=   0;   % [m] left side of model
xRight=xLeft+900; % [m] right size of model
HL    =10.08; % [m] prescribed head at left-most top cell
HR    = 0.113;% [m] prescribed head at right-most top cell
Ss    = 0;    % dummy specifice storage coefficient
dRhodC=(rhob-rhof)/(Cbrine-Cfresh)

save Cbrine.mat Cbrine
%% Mesh

xGr=xLeft:20:xRight;
yGr=[0 1];
zGr=Ztop:-4:Zbot;
    
[xGr,yGr,zGr,xm,ym,zm,DX,DY,DZ,NX,NY,NZ]=modelsize3(xGr,yGr,zGr);

Z=zGr(ones(NY,1),ones(1,NX),1:size(zGr,3));

%% Head along the top 

HeadTop=HL+((1:NX)-1)/(NX-1)*(HR-HL);

%% Generate all other matrices
M = ones(NY,NX,NZ);  % array of ones of the size of the model
IBOUND=M;
ICBUND=M;
STRTHD=M * strthd;
STCONC=M * Cfresh;
HK    =M * k;
PEFF  =M * peff;
SS    =M * Ss;

% Top of model, prescribed head (no need to prescibe concentration)
IBOUND(:,:,1)=-1;
STRTHD(:,:,1)=HeadTop;

% Bottom of model where C=Cbrine is prescribed
ICBUND(:,xm>300 & xm<600,end)=-1;
STCONC(:,xm>300 & xm<600,end)=Cbrine;
HK    (:,xm>300 & xm<600,end)= k/1e4;

% Bottom of model left and right of salt dome
IBOUND(:,xm<300 | xm>600,end)= 0;
ICBUND(:,xm<300 | xm>600,end)= 0;

% Vertical conductivity everywhere the same as horizontal concutivity
VK                           =HK;

