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

%  TO 091114 091203

basename='2D-Diagonal'; 

% The problem considered in this section is similar to the 2D-Uniform case
% except that the flow direction is oriented at a 45-deg angle to both
% the columns and the rows of the model grid. The grid Peclet number for this
% problem is 5 in the longitudinal direction and 50 in the transverse direction.
% The problem represents a challenging test for transport solution techniques
% because of the sharpness of the concentration front, compounded by the grid
% orientation effect.

% The model has 100 rows and columns and one layer. The parameters use
% are listed below. Note that mflab gets aL and aT from the worksheet LAY
% not from this m-file.

delx=10;      % m
dely=10;      % m
delz= 1;      % m
peff=0.14;    % [] effective porosity
v   =1;       % m/d seepage velocity 45degrees NE
aL  =2;       % set in worsheet LAY
              % aT/aL = 0.1; set in worksheet LAY
QIn =0.01;    % m3/d Injection rate
cIn =1000;    % ppm  Injection water concentration
              % simulation time 1000d (set in worksheet PER)
%% Mesh

NLAY=1; NROW=100; NCOL=100;
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(:))];

[xGr,yGr,z,xm,ym,zm,Dx,Dy,Dz,Nx,Ny,Nz]=modelsize3(xGr,yGr,z);

[Xm,Ym]=meshgrid(xm,ym);  % centers of all cells
xUR=xGr(end); yUR=yGr(1);

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

%% Uniform flow field under 45degrees in NE direction

% The uniform ambient flow is at 45 degrees to with the grid, but Zheng
% specifies no boundaries. Hence we have to estimate them ourselves. The
% problem is independent of the conductivity (it is assumed uniform) and
% iwth the given groundwater seepage of 1 m/d we have a gradient i=q/k
% Assuming Zheng uses the word seepage for Darcy flow, we have with k=10
% i=1/10=0.1 m/m. So the gradient along the 45degree boundary is sqrt(1/k).
% so starting with a head 0 at the upper corner of the model, we have the
% following boundary head everywhere
%STRTHD=0+(xGr-xUR)*(qx/k)+(yGr-yUR)*(qy/k), with ix=iy=-sqrt(q/k)

angle=pi*45/180; % degrees
q =v*peff;
qx=q*cos(angle);  % from seepage velocity to specific yield in x-dir
qy=q*sin(angle);  % from seepage velocity to specific yield in y-dir

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

STRTHD=0-(Xm-xUR)*qx/kh-(Ym-yUR)*qy/kh;

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

%% Wells
% we only need the well
iPer=1;                 % stress period (only one)
welpos=[1 90 10];       % 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
