% MF_ADAPT -- Seawat Benchmark example
% See Guo and Langevin (2002) p72ff
%  TO 090101 091207

basename='Elder';

% The Elder problem was originally designed for heat flow (Elder, 1967),
% but Voss and Souza (1987) recast the problem as variable density
% groundwater problem in which fluid density is a funciton of salt
% concentration. The Elder problem is commonly used to verify variable
% density groundwater codes.
% The problem concerns a cross section, 600 m  wide and 150 m heigh. A
% constant concentration boundary is specified over 300 m at the center top
% of the model. Molecular diffusion is the sole mechanism for hydronamic
% dispersion during the simulation, which runs for 20 years. Salt diffusing
% into the model creates complex vortices that redistribute the salt mass.
% A concentant concentration boundary with value zero reighns the obottom
% of the model. Two outlet cells with constant head of zero are specified
% at the to corners, providing an outlef for fluid and salt.

%% Parameters

L     =300;     % m half width of model
H     =150;     % m height and top of model
peff  = 0.2;    % [ ] effective porosity
k     = 0.411;  % m/d hydraulic condictivity
Cbrine=285.7;   % kg/m3 brine concentration
Cfresh=0;       % kg/m3 freshwater concentraton
rhof  =1000;    % kg/m3 freshwater density
rhob  =1200;    % kg/m3 brine density
aL    =0;       % m longitudinal dispersivity
aT    =0;       % m transverse dispersivity
D     =0.308;   % m2/d molecular diffusion coefficient
hCorner=150;    % m (point water heads)


%% Mesh

xGr=-L:2*L/41:L;
zGr=H:-6:0;
yGr=[0 1];

[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:NZ+1);

%% Generate all other matrices
IBOUND=ones(NY,NX,NZ);
ICBUND=ones(NY,NX,NZ);
STRTHD=ones(NY,NX,NZ) * hCorner;
STCONC=ones(NY,NX,NZ) * Cfresh;
HK    =ones(NY,NX,NZ)   *k;
VK    =HK;
PEFF  =ones(NY,NX,NZ)   *0.1;

%% Adaptation to boundaries
IBOUND(:,xm<-L/2 | xm>L/2,1)= 0;
IBOUND(:,     [1 end]    ,2)=-1;   % only corner in to of model

ICBUND(:,xm>-L/2 & xm<L/2,1)=-1;
ICBUND(:,     :,        end)=-1;

STRTHD(:,[1 end],2)= 150;

STCONC(:,xm>-L/2 & xm<L/2,1)= Cbrine;
