%% Saltwater Intrusion Packge, example 2

% The second examples is a varition of the first example, but includes three
% zones and no constant flow. Consider two-dimensional, confined flow in a
% vertical cross-section. The origin of a Cartesian x,z, coordinate system
% is chosen at the uper left hand  corner of the section; the aquifer is 40
% m thick. The hydraulic conductivity is k=2 m/d and the effective porosity
% is 0.2. A 300 m long section of the aquifer is considered and all
% boundaries are impermeabile. The groundwater is divided into 3 zones:
% fresh, brackish, and salt water, with dimensionless desnities of 0,
% 0.0125 and 0.025, respectively; at first, flow is treated as stratified.
% Initially, at time t=0, both interfaces are straight and make a 45 degree
% angle with the horizontal and run from (x,z)=(150,0) to (x,z)=(190,-40),
% and from (x,z)=(110,0) to x,z)=(150,-40). The head is fixed to 0.05 at
% the top ceell 1. The brackish zonee will rotate to a horizontal position
% through time.
% The Domain is discretized intno 60 cells of 5 m wide, and the timestep is
% 2 d. There are three zones and thus the initial elevations of two
% surfaces area specified. The shape of the brackish zone is computed at
% t=2000 days, by taking 1000 time steps. The maximum slop of the toe and
% tip is chosen to be 0.4, and the deltazeta and deltazetamin parameters
% are computed according to equation 16. All sources and sinks are fresh
% water.

% For illustration puposes, the same porblem is solved with the variable
% density option. This requires tow changes in the SWI input file. First,
% the ISTRAT parameters must be set to 0 to indicated variable density flow
% (rather than stratified flow). Second, the NU variable now needs 4 input
% values: the value of nu along the top of zone 1, the tow surfaces and the
% bottom of zone 3, so that the fourth lien in the swiex2.swi files
% becomes: 0 0 0.025 0.025. Notice that this gives an average value
% NU=0.0125 in the brackish zone.

basename='swiex2';

%% Specifics of this example
k     =2;    % hydraulic conductivity
peff  =0.2;  % effective porosity
strthd=0.05; % specified head in cell at right
Qo    =2;    % source in cell 1

%% Specify mesh
xGr=0:5:300; yGr=[0 2]; zGr=[0 -40];

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

Z=zGr(ones(NROW,1),ones(1,NCOL),1:NLAY+1);

%% Matrices of computation grid
IBOUND=ones(NROW,NCOL,NLAY);
HK    =ones(size(IBOUND))*k;
VK    =HK;

%% Head boundary conditions
IBOUND(:,1,1)=-1;                % fixed head at left hand side
STRTHD=ones(size(IBOUND))*strthd;  % start using fixed head

%% SWI specific matrices
SSZ=ones(size(IBOUND))*peff;     % effective porosity called SSZ in SWI
ISOURCE=zeros(size(IBOUND));     % source zone incase of sinks 

%% ZETA planes
NSURF=2;  % We have only one surface here, so NSURF=1 in this case

% A plane with elevations is necessary for all interfaces
% We just specify this single interface as simply as possible using
% coordinates in space (x,z)
ZETA(1).x=[0 150 190 300];  % x points for interface
ZETA(1).z=[0  0  -40 -40];  % z points for interface
ZETA(2).x=[0 110 150 300];  % x points for interface
ZETA(2).z=[0  0  -40 -40];  % z points for interface

% ZETA must be specified for every interface in all layers, but we have
% only one layer in this example. We use the above specified surface and
% interpolate its elevations to the cell centers of all (1) layers
for iSurf=1:NSURF
    ZETA(iSurf).values=interp1(ZETA(iSurf).x,ZETA(iSurf).z,xm); % zeta all cells
end
% Actual distribution over the layers is done by writeSWI