% mt3dms benchmark problem:
% One-dimensional Transport with Nonlinear or Nonequilibrium Sorptions
% See  Zheng & Wang (1999) p135ff
%  TO 091112

basename='1D-nonlin';

% Grove and Stollenwerk (1984) presented a computer code for modeling
% one-dimensional advective-dispersive transport with nonlinear equilibium
% controlled sorption and ion-exchange. MT3DMS is used here to solve the
% same problem involving Freundlich and Langmuir isotherms.
%
% The entry concentration changes at time zero from zero to a constant value
% and returns to zero after a fixed time.
%
% The list of parameters to sovle the problem is given here:

kh= 10; kv=10;   % values are irrelevant in this 1d problem
peff=0.37;       % []
dcol= 0.16;      % cm
vx=0.1; % cm/s,  % true velocity = vx/peff=
rhob=1.587;      % g/cm3                --- rhob in LAY worksheet
Kfreundlich = 0.3; % (ug/g)(L/mf)^alfa  --- SP1 in LAY worksheet
alfa = 0.7;        % []                 --- SP2 in LAY worksheet
KLangmuir = 100;   % L/mg               --- SP1 in LAY worksheet
Slangmuir = 0.003; % ug/g               --- SP2 in LAY worksheet
Co = 0.05;       % mg/L
t0  = 160;       % s -- pulse duration  --- PER in PER worksheet

% The values here are not directly used except for the conductivities. We
% use the values given here in the worksheet LAY where we specify the
% sorption properties. t0 is used in the worksheet PER to define the length
% of the first stress period. The concentration Co is used in the
% definition of the PNTSRC for the SSM module of MT3DMS. vx, the seepage
% velocity is used below to compute the inflow in the left most cell, after
% which this value, then a WEL, was used in the WEL worksheet for the
% WELL-flow.

% We cannot use different sorption processes in different layers because
% the ISOTHM property in the RCT module is set for the entire model for all
% species; only the spcific process variable values are set per layer or
% even per cell. Hence all species in the model are either subject to
% Freundlich or Langmuir sorption and, therefore te simulate both processe,
% the model has to be run at as many times. However, we may combine
% sorption with an irreversible reaction such as decay.

% Hence, contrary to the previous example, there is no value in have more
% than a single layer in the model.

%% Mesh
% Get the number of active layers from the excel sheet
[LAYparnams,LAYparvals]=getExcelData([basename,'.xls'],'LAY','Horizontal');
LAYparvals(LAYparvals(:,1)<1,:)=[]; % remove inactive layers (with layer number <1)

NLAY=size(LAYparvals,1);

% Specifying layers is very flexible in mflab.
%  -- only specify the layers that are different. If all layers are the
%  same than specifying one layer is enough (if layer number is 1).
% --- use layer number in your specification
%    ---- any layers with layer number <1 and > NLAY are ignored
%    ---- this allows easy switching between alternatives
%    ---- layers don't need to be sorted, but there must always be one
%    ---- with layer number 1 or else some layers will not be defined.
%    ---- Use omitted layers will be copied from the defined layer with the
%    ---- highest number < than the actual layer. So in a model of 100
%    ---- layers where a change in parameters is at layers 1 (always) and
%    ---- at layers 32 and 54, and only layers 1, 32 and 54 need to be
%    ---- define. Layers 2 - 31 will be the same as layer 1, layers 33 to 53 the
%    ---- same as layer 32 and the layers beyond 54 will be the same as layer
%    ---- 54.

% The worksheet LAY has 2 layers one numbered 0 and one numbered 1. One is
% for the Freundlich isotherm and one for the Langmuis isotherm. To switch
% between them change the layer numbers to make to correct one active and
% in worksheet MT3D chagne the value of ISOTHM to 2 for Freundlicht and to
% 3 for Langmuir.

%% The computational mesh

NROW=1; NCOL=70;            % MESH definition
DELX=ones(1,NCOL) * dcol;           % in cm, note DELC is column width
DELY=1;                             % in cm, note DELR is row width
DELZ=ones(NLAY,1) * 1;              % in cm

%% Grid xGr and yGr, centralizing the fist cell center at xGFr=0
xGr=-DELX(1)/2 + [0 cumsum(DELX)];   % center of first cell at xGr=0
yGr=[0 cumsum(DELY)];

z=sum(DELZ)+ [0; -cumsum(DELZ(:))]'; z=reshape(z,[1,1,length(z)]); 
Z=z(ones(NROW,1),ones(NCOL,1),:);

%% Generate all other matrices
IBOUND=ones( NROW,NCOL,NLAY);  IBOUND(:,end,:)= -1;  % fix right hand side
ICBUND=ones( NROW,NCOL,NLAY);
STRTHD=zeros(NROW,NCOL,NLAY);

STCONC{1}=ones(NROW,NCOL,NLAY)*0;  % start concentration species 1 for Freundlich

HK    =ones(NROW,NCOL,NLAY)   *kh;   % horizontal conductivity
VK    =ones(NROW,NCOL,NLAY)   *kv;   % vertical   conductivity
PEFF  =ones(NROW,NCOL,NLAY)   *peff; % effective porosity

%% Wells at tehe west boundaries
% Using the WEST and EAST list it is now straigthforward to generate the
% list of wells, i.e. cells with given flow input. These are the cells at
% the west end of the model, as we fix the head at the right end.

q=vx.*DELZ.*DELY*peff;  % compute the cell inflow at the west end from the given seepage velocity

% The result of this q was put in the worksheet.
% As will be clear, defining WEL and other boundary conditions in the
% worksheet, as it is non longer transparant how parameters are defined and
% it is contrary to the desire to script the modeling process for the sake
% of  reproducibility.
% However, we use the worksheet to demonstrate that it works.
%
% Notice that we need to define the WEL for two stress periods, the same
% number as in the PER sheet, where the stress periods are defined  

%% PNTSCR
% The SSM package of MT3DMS needs the pont sources that carry a
% concentraton upon infow. The only point sources to be defined are at the
% west end of the mdodel, i.e. the wells (ITYPE=2, see SSM documentation):
%
% In this example we also specify the point sources for MT3DMS in the
% workbook worksheet "PNTSRC"
%
% Notice that, because the inflow concentration in the second stress period
% is zero, we only need to define the PNTSRC for the first stress period as
% was done in the worksheet.