%% Circular polder

basename='Polder';

% The problem here is a circular polder with a radius of 1500 m situated
% above a semi confined aquifer of 200 m thickness.
% The aquifer is initially filled with brackish water
% having a concentration of 6 kg/m3. The polder has a radius RP=2000 m. The head at the top is
% maintained at -4 m within this radius and at 0 m outside this radius. The
% resistance against vertical flow at the op of the aqufier is 100 days in
% and 100 days outside the polder. Horizontal conductivity is 30 m/d and
% vertical conductivity 6 m/d. The effective porosity is 0.35. The density
% is computed from the concentration using a liner relationship yielding
% 1000 kg/m3 for fresh water and 1025 kg/m3 voor ocean water at 35 kg/m3.
% The question to solve here is how the salt water distribution evolves
% during the first 100 years after the polder was established.
% In this synthetic situation, the only boundary conditions are at the top
% of the system as indicated. The sides of the model and the bottom are
% no-flow boundaries.
% The model is 10000 by 10000 m. We use cells of 250x250x10 m uniformly.
% The situation is computed using a single stress period of 100 years,
% steady state, in time steps of 2 years, which are also saved for showing
% the dynamics of the development by means of a movie.
% One of the questions is how the saltwater and the freshwater discharge
% varies over time
% TO 091226  100313

%% specifics

peff=0.35; % effective  porosity
kh  =30;   % horizontal conductivity
kv  =kh/5; % vertical   conductivity
hIn =-4.0; % head inside polder
hOut= 0.0; % head outside polder
cOut=100;  % resistance cover layer outside polder
cIn = 50;  % resistance cover layer inside polder
RP =1500;  % polder radius
CONC_S=6;  % concentration of saltwater
CONC_F=0;  % concentration of freshwater

%% The grid
xGr=[-5125:250:-1750 -1750:125:1750 1750:250:5125];  % boudaries exactly on -5000 5000 and one cell in center
yGr=[-5125:250:-1750 -1750:125:1750 1750:250:5125];
zGr=permute(0:-10:-200,[1,3,2]);

%% Grid housekeeping
[xGr,yGr,zGr,xm,ym,zm,DELX,DELY,DELZ,NCOL,NROW,NLAY]=modelsize3(...
    xGr,yGr,zGr);

[XM,YM,ZM]=meshgrid(  xm,  ym,  zm);  % we use full matrices in GHB definition
[DX,DY,DZ]=meshgrid(DELX,DELY,DELZ); % same

Z=zGr(ones(NROW,1),ones(1,NCOL),(1:length(zGr(:)))); % full 3D Z-grid always necessary

%% Aquifer
IBOUND=ones(NROW,NCOL,NLAY);  % no fixed heads necessary because we use GHB
ICBUND=ones(size(IBOUND));   % no fixed conc necessary because inflow is fresh (C=0)
STRTHD =zeros(size(IBOUND)); % immaterial as we use steady state and GHB

kh_    =ones(1,1,NLAY)*kh;  % turn kh into a 3D vertical vector
kv_    =ones(1,1,NLAY)*kv;  % same for kv
HK     =kh_(ones(NROW,1),ones(1,NCOL),1:NLAY);  % generate 3D kh array
VK     =kv_(ones(NROW,1),ones(1,NCOL),1:NLAY);  % same for kv

PEFF   =ones(size(IBOUND))*peff;
STCONC=ones(size(IBOUND)) * CONC_F; STCONC(ZM<=-0)=CONC_S;

%% GHB --  the reistance at the top is modeled using the GHB package
R = sqrt(XM.^2+YM.^2); % distance from polder center, 3D matrix

Iout=find(R> RP & ZM>zGr(2));  % 3D indices
Iin =find(R<=RP & ZM>zGr(2));  % 3D indices

% Layer Row and Column indices of these areas
LRCin =cellIndices(Iin ,size(IBOUND),'LRC'); uin =ones(size(LRCin(:,1)));
LRCout=cellIndices(Iout,size(IBOUND),'LRC'); uout=ones(size(LRCout(:,1)));

vcontIn =DX(Iin ).*DY(Iin )/cIn;
vcontOut=DX(Iout).*DY(Iout)/cOut;

[PERnams,PERvals]=getPeriods(basename);
NPER=size(PERvals,1);

GHB=[ uin  LRCin  uin *hIn  vcontIn; ...
      uout LRCout uout*hOut vcontOut];
  
% Subsequent boundary conditions as previous
GHB=[GHB; -(2:NPER)' * ones(size(GHB(1,:)))]; % as previouw

%% Saving under the table some extra values needed to present the results
save pol.mat R RP Iout Iin hIn hOut