% Salt dilution test simulated with MT3D

% TO 091029 091204
%
% We simulate the flow in a pipe with a square cross section which is
% connected to a box (cube) and a subsequent pipe with the same XSection.
% The right end of the model has a fixed head.
% The left end has a fixed flow boundary.
% Other boundaries are all closed
% Initial concentrations are zero, the flow is steady state.
% To simulate flow in a pipe and box we set porosity to 1.
% The idea is to inject salt during a short period and measure the
% salinity curve as it passes by observation points in the system.
% We then use these curves to extimate the total flow through the system,
% which is the essence of a salt dilution test.
% This configuration resemples the flow in extraction galleries in the
% Amsterdam Water Supply Dunes in which we carried out salt dilution tests
% to measure the flow along the 300 m long gallery at its obsrvation
% points, which are 20 cm wide pipe connected to the mentioned bos.
% Clearly real flow in the gallery is turbulent, which we cannot simulate
% with a groundwater transport model. however we may see if putting the
% salinity sensors at different places in the box. In this simulation
% the flow situation is completely known and therefore suitable to test.
%
% We run three stress periods. First 5 seconds with zero concentration
% input and then the salt is dosed during 5 seconds, after which the input
% concentration is again zero.
% After 5 seconds salt is dosed during 5 seconds. Then the salt is washed
% out by the flow.
% Output is produced at every second.
% You may want to experiment with dispersion etc.
% You may also add observation points or show the concentration in certain
% points.
%
% TO 091204
%
%% recharge on a rectangular area with fixed head boundary at left and right
clear variables;
close all;

basename='SaltTest';

%% The model grid
% The pipes are 15 cells long as is the box:
NROW=15; NCOL=3*15; NLAY=15;
rib=0.04; % m                    rib is rublengte van model cell

xGr=(0:NCOL)*rib;
yGr=(0:NROW)*rib;
z=(NLAY:-1:0)*rib;   z=z-mean(z);

[xGr,yGr,z,xm,ym,zm,DELX,DELY,DELZ,Nx,Ny,Nz]=modelsize3(xGr-mean(xGr),yGr-mean(yGr),z-mean(z));

[X, Y, Z ]=meshgrid(xGr,yGr, z);
[XM,YM,ZM]=meshgrid(xm ,ym ,zm);
[DX,DY,DZ]=meshgrid(DELX,DELY,DELZ);

%% fixed head boundaries. Are interpreted as local point heads in fdm2dens
STRTHD=zeros(NROW,NCOL,NLAY);

%% HK,VK according to layer specification given above (sec)
HK =ones(NROW,NCOL,NLAY);  % doesn't matter here
VK =ones(NROW,NCOL,NLAY);

IBOUND=zeros(NROW,NCOL,NLAY);  % initialize all inactive cells
IBOUND(6:10,  :  ,6:10)= 1;    % activate the pipe cells
IBOUND( :  ,16:30, :  )= 1;    % activate the box  cells
IBOUND(6:10,end  ,6:10)=-1;    % set right boundary to fixed heads

%% Visualisation of cross section

figure;
contour(permute(XM(8,:,:),[3,2,1]),permute(ZM(8,:,:),[3,2,1]),permute(IBOUND(8,:,:),[3,2,1]),[0.5,-0.5]);  % show layers
hold on

xlabel('xGr [m]'); ylabel('z [m]'); title('Cross section');

ICBUND=ones( size(IBOUND));
PEFF  =ones( size(IBOUND));
STCONC=zeros(size(IBOUND));

%% Get NPER from workbook
[PERparnams,PERparvals]=getExcelData([basename,'.xls'],'PER','HOR');
NPER=size(PERparvals,1);
% To specify this injection concentration in a most general fashion,
% we add a new variable to the PER sheet, called cInj and specify the
% tracer concentrations for all stress periods (either zero or the value we
% desire. Having done so, we can retrieve the injected tracer concentration
% as follows:
CInj=PERparvals(:,index('CInj',PERparnams));

% The duration of the stress periods comes from
PERLEN=PERparvals(:,index('PERLEN',PERparnams));

% Ok but how much mass have we injected?
% To compute this, we use the water velocity, the cross section, which is
% given here:
v0 =0.15;       % m/s   water velocity

% To compute the cross section of the pipe, we make use of IBOUND, which
% has value 1 for active and 0 for ineactive cells. There are different
% ways to compute this cross section
% One way is to turn a cross seciton toward us (align z vertically and rows
% from left to right while taking a cross section (a given column, for
% instance the second one). This is done using permute, which yield the
% cross section plane wigh zeros for the inactive cells and ones for the
% active cells. To show it try:
spy(permute(IBOUND(:,3,:),[3,1,2])); xlabel('iy'); ylabel('iz');

% If we multiply this element by element (.*) with the cross section area
% of the individual cells in the plance, obtained by the matrix
% multiplicaton (*) of DELZ vertically and DELY horizontally, wer're done:
A=sum(sum(permute(IBOUND(:,3,:),[3,1,2]).*(DELZ(:)*DELY')));
% this is somewhat complicated. There is a somewhat easier way to compute
% the cross section area as is shown in the next section

% Ass all cInj in stress periods without tracer are zero, we can just
% compute the entire injected mass over all stress periods as follows:
M=v0*A*sum(CInj.*PERLEN);  % g if conc is in g/cm3 here 

%% Mass loading at single entry cell or injection with a given

% How to inject the tracer?
% We may do it by means of a uniform concentration over the entire cross
% seciton, i.e. in the most ideal way or we may do it by injection pure
% tracer at one or more of the entry cells. This can be done by increasing
% the concentrataion at these points (ITYPE=2 in SSM) or by injection mass
% at these points (ITYPE=15) in the SSM module of MT3DMS.
% To demonstrate the two ITYPE options we simulate both methods by means of
% the switch mass_loading. If on, then mass is injected at a single cell at
% the entry face of the model, if off we have a constant tracer
% concentration across the entry face. The tracer is only injected during
% the second stress period.
%
% We use the same parameter in both cases. In the case of mass_loading CInj
% is interperted as g/s instead of concentration
%

mass_loading=1;   % use 0 for "off" and 1 for "on"

%% concentration during the second stress period across the entire zone

% Which are the left hand side cell in the pipe??
I=find(IBOUND==1 & XM<xGr(2));

% Alternative way to compute the injection mass
A=sum(DY(I).*DZ(I));

mask=ones(size(I));  % to swtich cells at model entry on and off

if mass_loading    % we select a random cell from the model entry face
    
    sr=sortrows([rand(size(I)) I]);
    irand=sr(1,2); % random cell of model entry face

    mask(I~=irand)=0;  % switch off all other entry face cells

    % Compute the tracer entry mass in each stress period
    M=     sum(CInj.*PERLEN); % where now CInj is interpreted as mass inje 
else
    % If no mass_loading, treat CInj as uniform conc over the entry face
    M=v0*A*sum(CInj.*PERLEN);
end

% Total discharge of model given v0 and configuration
Q0=v0*DY(I).*DZ(I); % m3/s

save Q0.mat Q0 M      % for later use in mf_analyze

% Get a list of cells at the pipe at left hand size of the model
LRC=cellIndices(I,size(IBOUND),'LRC');  % Use cellIndices to get a list of L R C

if mass_loading
    ITYPE=15; % ITYPE=15 in SSM package means mass loading source
else
    ITYPE= 2; % ITYPE=2  in SSM package means "well" 
end

WEL=[];                % initialize list of wells
PNTSRC=[];             % initialize list of pointscources for SSM module 
u=ones(size(I));       % a column of ones as long as I
for i=1:NPER           % Then for each stress period
    WEL   =[WEL;   [i*u, LRC Q0]                       ]; % generate wells
    PNTSRC=[PNTSRC;[i*u, LRC (u.*mask)*[CInj(i) ITYPE CInj(i)]]]; % generate pointsources
    % NOTE that ITYPE signals whether or not mass_laoding is on
    % and CInj is mass injection if mass loading is on
    % mask switches all point sources of exectp the selected one
end
PNTSRC(PNTSRC(:,6)==0,:)=[];

% note that mask contains the information on mass loading, in fact that
% only a single point will get tracer injected (in case of mass loading)

% That's it, model defined !

% A convenient way to run the model including the analysis is this one
% close('all'); cleanup; mf_setup; mf_analyze