% MT3DMS Benchmark Problem see Zheng & Wang (1999), p140ff
% Two-dimensional Transport in a Radial Flow Field

%  TO 091114 091203

basename='2D-Radial'; 
problem='2D radial flow field (Zheng p140)';

% The test problem considered in this section concerns the two-dimensional
% transport of solute injected from a fully penetrating well. The problem
% is intended to test the accuracy of MT3DMS as applied to a radial flow
% system. The assumptions for this problem are:
% -- The injection rate of the well is constant
% -- The ambient groundwater velocity is negligible relative to the
%    velocity created by the injection
% -- The aquifer is homogeneous, isotropic, and infinite in areal extent.
% -- The flow fiel is steady state.

delx=10;      % m
dely=10;      % m
delz= 1;      % m
peff=0.3;     % [] effective porosity
aL=10;        % set in worsheet LAY
              % aT/aL = 1.0; set in worksheet LAY
QIn =100;     % m3/d Injection rate
cIn =1;       % ppm  Injection water concentration
              % simulation time 27d (set in worksheet PER)
%% Mesh

NLAY=1; NROW=31; NCOL=31;
DELX=ones(1,NCOL)* delx; xGr=            [0 cumsum(DELX)];
DELY=ones(NROW,1)* dely; yGr=sum(DELY) - [0; cumsum(DELY(:))];
DELZ=ones(NLAY,1)* delz; z  =sum(DELZ) - [0; cumsum(DELZ(:))];

% housekeeping the mesh. Notice that the mesh is centered around its center
% see input of modelsize
[xGr,yGr,z,xm,ym,zm,Dx,Dy,Dz,Nx,Ny,Nz]=modelsize3(xGr-mean(xGr),yGr-mean(yGr),z);

Z=NaN(NROW,NCOL,NLAY+1); for i=1:size(Z,3), Z(:,:,i)=z(i); end

%% Generate all other matrices
IBOUND=ones(NROW,NCOL,NLAY);
   IBOUND(:,[1 end],:)= -1;  % Fix west and east   boundary head
   IBOUND([1 end],:,:)= -1;  % Fix north and south boundary head
ICBUND=ones(NROW,NCOL,NLAY); % All one, don't need fixed conc boundary

%% Other inputs

kh= 10; kv=0;   % values are irrelevant in this 1d problem

HK       =ones(NROW,NCOL,NLAY)*kh;
VK       =ones(NROW,NCOL,NLAY)*kv;
PEFF     =ones(NROW,NCOL,NLAY)*peff;
STRTHD   =zeros(NROW,NCOL,NLAY);
STCONC{1}=ones(NROW,NCOL,NLAY)*0;

save cIn.mat cIn  % shore cIn because we need it in mf_analyze later on

%% Wells
% we only need the well
iPer=1;                 % stress period (only one)
welpos=[1 find(xm>=0,1,'first') find(ym<=0,1,'first')];       % position of the well  Layer Row Col
WEL=[iPer welpos QIn];  % well input

%% PNTSCR, to give wells their input concentration

%       iPer L R C CSS ITYPE CSSM(1:NCOMP)
PNTSRC=[iPer welpos cIn 2 cIn]; % well at 1 16 16, itype 2
