% Example Boussinesq (flow on an sloped base)
% TO 100414 100507

% Boussinesq studied the flow of unconfined groundwater on a sloped aquifer base.
% This can be modeled with modflow by letting the cells step down. However
% this way of modeling is not correct as follows from the math involved.
% The steerper the slope the more the two deviate. The difference is
% systematic and cannot be compensated for by reducing the size of the
% Modflow cells. See userguid for the background.
%
% This file implements a simple one-dimensional unconfined model (along the
% x-axis which steps down along a given slope (note that the waviness of
% the slope can be adjasted with the parameter A, A=0 giving a straight slope.
% The mfLab user guide compares the results obtained by MODFLOW with the
% analytical solution given by Steward in WRR (2007). Full analysis is given
% in the mfLab user guide.
%
% Note that the drying and rewetting mechanism faile easily in this type of
% problems, causing lack of convergence. This is a long standing
% frustrating issue with respect to MODFLOW. Others (Doherty (2000) and
% Painter (2008) have developed better, more robust mechanisms, which are
% also included in the MODFLOW that comes with GMS (Doherty's version) and
% within MODSURFACT. It is really a pity that this issue has not yet be
% decently solved.
%
% In the file below, convergence does occur, but only for as long as the
% parameter D, which defines the start heads at the boundaries is large
% enough (here 100 m). If you choose it to be lower than say 75 m,
% convergence fails. I found no obvious way to overcome this with the
% mechanisms provided by MODFLOW (see description in the mfLab user guide).
%
% To show that a solution is possible even under difficult situaitons, such
% as with a smal to medium stong precipitation on the slope, I implemented
% the steady state 3D finite difference model in Matlab, called fdm3.m in the
% separate m-script in this directory called
% DohertyBoussinesqRecharge.m
% which implements Doherty's method to mitigate non-linaerities due to the
% mentioned problems and which prevents cells from running dry altogether.
% The results of this are also described in the mentioned user guide.
%
% The current model defined below and in the accompanying workbook
% Boussinesq.xls, is transient, as this allows a gradual movement of some
% initial head configuration to the steady-state solution, thus preventing
% strange oscillations that might occur when solving the final steady-state
% solutions directly.
% 
% The model uses one stress period and one layer. The WETDRY values are
% specified below instead of in the workbook, thus allowing cell-by-cell
% values in the input and not just one value per layer.


clear variables;

%% The model name. Every model has a "basename" and all associated files
basename='Boussinesq';

clr='brgkmc';
k=2;   % conductivity (uniform)
Sy=0.1;
Ss=0.0005;
dx=1;  % grid cell width (uniform)
z0=0;  % z at xm(1)
D=100;  % dummy aqufier thickness (because unconfined)
slope =-1/2; % dB/dz, inclination with B the bottom of the aquifer
Q0 = 5;% m2/d discharge
%kstar=k;  % k*cos(theta)^2;
%% Specify grid line coordinates
xGr=0:dx:1000; % [m]
yGr=[-0.5 0.5];  % [m]
[xGr,yGr,xm,ym,DX,DY,Nx,Ny]=modelsize(xGr,yGr);

%% Aquifer bottom, top and center elevation

A=40; % set A=0 for constant slope

zB=z0+slope*(xm-xm(1))+A*sin(2*pi*(xm-xm(1))/(xm(end)-xm(1))*2);  % Aquifer bottom
%SLOPE=slope+A*cos(2*pi*(xm-xm(1))/(xm(end)-xm(1))*2)*2*pi/(xm(end)-xm(1)*2);

zT=ones(size(xm))*D;
Nz=1;

Z=NaN(Ny,Nx,Nz+1);

Z(:,:,1)    =zT;
Z(:,:,2)    =zB;
Z(:,:,end)  =zB;

Dz=-diff(Z,1,3);

ZB=Z(:,:,2:end);                    % make sure the second aquifer is always filled but has small k

zm=0.5*(Z(:,:,1:end-1)+Z(:,:,2:end));

%% Arrays
IBOUND=ones(Ny,Nx,Nz);  IBOUND(:,[1 end],1)=-1; % only the first layer

HK    =ones(Ny,Nx,Nz)*k; % hor conductivity      [m/d]
VK    =HK;                   % vertical conductivity [m/d]

SS    =ones(Ny,Nx,Nz)*Ss;
SY    =ones(Ny,Nx,Nz)*Sy;

%% We set the starting heads equal to the normal depth
STRTHD=ones(Ny,Nx,Nz)*D; STRTHD(IBOUND==-1)=ZB(IBOUND==-1)+D;

WETDRY=cell(Nz,1);
WETDRY{1}=0.05;

%% we try to put the fixed flow at x=0 or at x=L in both cases set h equal
iPer=1; iy=1; ix=1; iz=1;
WEL=[iPer iy ix iz Q0];
