# Groundwater

The MACM 416 project page for Galen Marisco.

## Background

Since the late 1800's, mathematical models of groundwater flow have been used in the planning of well construction, and more recently the analysis of contaminant dispersion in aquifers. Before the advent of high powered computers, engineers and mathematicians were forced to make very restrictive assumptions in their models so that the problems would be tractable. Through the use of numerical methods however, modeling groundwater systems has come a long way such that consistent and powerful predictions of groundwater flow and contaminant dispersal can be made and used in an industrial setting.

This topic has special importance to British Columbians as B.C. has one of the largest concentrations of natural-gas in Canada as well as having a sizable portion of the worlds fresh water stored in its lakes and aquifers. With the demand for cheap energy and water as high as it is in current times it seems inevitable that this gas and water will be extracted at an increasing rate. The process of natural-gas extraction is well known in certain cases to cause contamination of groundwater systems(see fracking) and thus it is essential to properly model our groundwater systems if we wish to preserve them from contamination. Similarly, over extraction of water from our aquifers can have irreversible consequences on our water supply. Through proper modeling which relies on powerful numeric approximations of how contaminants disperse through groundwater systems and how aquifers react to welling, it may be possible to mitigate the pollution of natural-gas extraction and preserve B.C.'s abundance of fresh water.

## Project Outline

The focus of this project is to develop several finite difference models of groundwater flow about a pumping well. Initially, a simple model will be computed while subsequent models will each have an added level of complexity. All of the problems will deal with a simple rectangular domain centered around the well located at $$(x,y)=(x_0,y_0)$$. Initially we will perform some preliminary computations on the simple problem of well draw-down on a confined aquifer which is homogeneous and isotropic, after which we will proceed with the more interesting problems dealing with unconfined aquifers. The first (preliminary) model will be of an confined aquifer flow about a well. The hydraulic conductivity will be assumed to be homogeneous in this model for simplicity. The problem will first be solved in steady-state \begin{align} \frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial y^2} = \frac{R}{T} \tag{1} \end{align} \begin{align*} T = K b, \ R = \cases{ \frac{Q}{\Delta x \Delta y} & \text{if } (x,y) = (x_0,y_0)\cr 0 & \text{if } (x,y) \ne (x_0,y_0) }. \end{align*} After this we will move on to the transient problem \begin{align} \frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial y^2} = \frac{S}{T}\frac{\partial h}{\partial t} - \frac{R}{T} \tag{2} \end{align} which will be solved with an implicit Euler method. Once we have finished examining the confined aquifer we will move on to the unconfined aquifer in transient state.Our model is given by \begin{align} \frac{ K }{ 2 } (\frac{ \partial^2 v }{ \partial x^2 } + \frac{ \partial^2 v }{ \partial y^2 }) = \frac{S}{2\sqrt v}\frac{\partial v}{\partial t} -R(x, y)\tag{3} \end{align} \begin{align} v = h^2 \end{align} which we will solve using the Crank-Nicholson method. \begin{align} &h &=& \text{hydraulic head},\\ &K &=& \text{hydraulic conductivity},\\ &S &=& \text{storage coefficient},\\ &R &=& \text{recharge},\\ &Q &=& \text{flow rate of well},\\&T &=& \text{transmissivity},\\&b &=& \text{saturated thickness of aquifer}.\\ \end{align} Further, model calibration and verification will be discussed.

Using *Darcy's Law we can derive our model as such:

\begin{align} &q = -K \frac{dh}{dl}\\ \end{align} And in two dimensions: \begin{align*} &q_x = -K \frac{\partial h}{\partial x},&q_y = -K \frac{\partial h}{\partial y} \end{align*} Or in Vector Notation: \begin{align} {\bf q} = -K{\bf grad}h \end{align}

Now, if we assume water is incompressible and that conservation of mass holds, we get:

\begin{align} \frac{\partial q_x}{\partial x} + \frac{\partial q_y}{\partial y} = 0 \end{align} We can substitute Darcy's law into the above equation, piece by piece: \begin{align*} \frac{\partial }{\partial x}(-{\bf K}\frac{\partial h}{\partial x}) + \frac{\partial }{\partial y}(-{\bf K}\frac{\partial h}{\partial y}) = 0 \end{align*} Where: \begin{align*} {\bf K} = K(x,y) \\ \end{align*} If $${\bf K}$$ is independent of $$x, y$$, ie. if the medium is assumed homogeneous & isotropic, then the equation reduces to Laplace's equation: \begin{align*} \frac{\partial^2 h }{\partial x^2} + \frac{\partial^2 h }{\partial y^2} = 0 \end{align*}

## Preliminary Results

We start with a simple steady-state model of well draw-down on a confined aquifer which is both homogeneous and isotropic. In this model, the pumping well is located in the center of a $$4000 \cdot m^2$$ square field and is being pumped at a rate of $$2000 m^2 \cdot day^{-1}$$. It is assumed that the well is fully penetrating the aquifer which has a uniform depth $$b=100m$$. We can then assume that flow to the well will be horizontal and we can treat the well as a point sink. There will be no recharge from precipitation so $$R$$ will be zero everywhere outside the infinitesimal volume about the well. Observation wells along the outer edge of the region have been hypothetically measured to have head values of $$h = 10 m$$ and it will be assumed that the draw-down has no effect on the aquifer outside of this region. For a more accurate description, the boundary of the draw-down would be circular but for now we will approximate with a square domain.Our governing equation is given by $$(1)$$ and we can summarize our boundary conditions and parameters as follows: \begin{align} h = 10, \mbox{ on } x = 0,\ x = 4000\\ h = 10, \mbox{ on } y = 0,\ y = 4000\\ \end{align} \begin{align*} R = \cases{ \frac{Q}{\Delta x \Delta y} & \text{if } (x,y) = (2000,2000)\cr 0 & \text{if } (x,y) \ne (2000,2000) } \end{align*} We take sane parameters from : \begin{align} K &= 30m \cdot day^{-1},\\ b &= 10m,\\ T &= Kb = 300m^2 \cdot day^{-1},\\ Q &= 2000m^3 \cdot day^{-1},\\ \Delta x &= \Delta y = 200 m.\\ \end{align}

To approximate the solution we use a centered finite-difference scheme: \begin{align} \frac{H_{i+1,j} + H_{i-1,j} + H_{i,j-1} + H_{i,j+1} - 4H_{i,j}}{h^2} = \frac{R_{i,j}}{T} \end{align} Using this scheme we can build the solution in matlab with the $$\verb* delsq*$$ command.
First, we set up our grid:

   M =   21;
L =   4000;
h = L/(M-1);
[xg yg] = meshgrid(0:h:L,0:h:L);
G = numgrid('S',M);
center = (M-1)/2 +1;
N = sum(G(:)>0);
rhs = ones(N,1);


Now we define our discrete laplacian

   lapl = sparse(-delsq(G)/h^2);


and our parameters $$Q$$, $$T$$ and $$R_{well}=R(2000,2000)$$

   Q = 2000;
T =  300;
R_well = Q/(h^2);


To specify our recharge $$R$$ as was discussed above, we proceed as follows:

   i = 2;
while i<=M-1
rhs(G(2:M-1,i)) = 0;
i = i+1;
end
rhs(G(center,center)) = R_well/T;


It can be seen that the recharge is specified everywhere in the domain to be zero except at the location of the well. We can now set up the boundary conditions:

   HB = 10*ones(M,1);
HL = 10*ones(M,1);
HR = 10*ones(M,1);
HT = 10*ones(M,1);

rhs(G(  2,2:M-1)) = rhs(G(  2,2:M-1))-HB(2:M-1)/h^2;%bottom
rhs(G(M-1,2:M-1)) = rhs(G(M-1,2:M-1))-HT(2:M-1)/h^2;%top
rhs(G(2:M-1,  2)) = rhs(G(2:M-1,  2))-HL(2:M-1)/h^2;%left
rhs(G(2:M-1,M-1)) = rhs(G(2:M-1,M-1))-HR(2:M-1)/h^2;%right


We want to solve our problem now so we use the built in matlab system solver:

   h_ = lapl\rhs;


Plotting the results, we see that the hydraulic head at the well is given by $$H(2000,2000) = 5.8m$$. A visualization of the data shows the solution behaving as expected.
Error creating thumbnail: File with dimensions greater than 12.5 MP

Though this model has issues, it gives us a general idea of the solution.

## Transient Well Draw-down on a Confined Aquifer

### Analysis

We continue developing our model of a confined aquifer of uniform depth $$b = 100m$$ by considering transient conditions. Our model is given by $$(2)$$.

We use the same values of $$Q, K, b$$ as we did in the steady-state problem but now we introduce the storage term $$S = 0.002$$. To discretize our problem, we use implicit-Euler for time($$k = t^{n+1}-t^{n}$$) and a centered difference in space: \begin{align} \frac{H_{i+1,j}^{n+1} + H_{i-1,j}^{n+1}- 4H_{i,j}^{n+1} + H_{i,j-1}^{n+1} + H_{i,j+1}^{n+1} }{h^2} &= \frac{S}{T}\frac{H_{i,j}^{n+1} - H_{i,j}^{n}}{k} -\frac{R_{i,j}}{T}\\ \implies\frac{H_{i+1,j}^{n+1} + H_{i-1,j}^{n+1}- (4+\frac{h^2S}{Tk})H_{i,j}^{n+1} + H_{i,j-1}^{n+1} + H_{i,j+1}^{n+1} }{h^2} &= \frac{-S}{Tk} H_{i,j}^{n} -\frac{R_{i,j}}{T} \end{align} Now we will simplify the equation a bit \begin{align} \lambda = 4+\frac{h^2S}{Tk},\ \alpha = \frac{R}{T},\ \beta = \frac{S}{Tk} \end{align} to get the following system: \begin{align} \frac{H_{i+1,j}^{n+1} + H_{i-1,j}^{n+1}- \lambda H_{i,j}^{n+1} + H_{i,j-1}^{n+1} + H_{i,j+1}^{n+1} }{h^2} = -\beta H_{i,j}^{n} -\alpha. \end{align} If we discretize our domain as a $$5\times 5$$ grid \begin{align} \matrix{ & \bullet_{1} & \bullet_{6} & \bullet_{11} & \bullet_{16} & \bullet_{21} & \cr & \bullet_{2} & \bullet_{7} & \bullet_{12} & \bullet_{17} & \bullet_{22} & \cr & \bullet_{3} & \bullet_{8} & \bullet_{13} & \bullet_{18} & \bullet_{23} & \cr & \bullet_{4} & \bullet_{9} & \bullet_{14} & \bullet_{19} & \bullet_{24} & \cr & \bullet_{5} & \bullet_{10} & \bullet_{15} & \bullet_{20} & \bullet_{25} &} \end{align} it can be seen that the points directly adjacent to the boundary will be affected by our Dirichelet boundary conditions of $$H = 10$$. Before setting up our differentiation matrix, let's examine how the boundary conditions affect points on the corner and on the side. Points $$H_7$$ and $$H_8$$ are good candidates:
@ $$H_7$$:
\begin{align} \frac{H_8^{n+1} + H_{12}^{n+1} - \lambda H_7^{n+1} + 10 + 10 }{h^2} &= -\beta H_7^{n} -\alpha\\ \implies\frac{H_8^{n+1} + H_{12}^{n+1} - \lambda H_7^{n+1}}{h^2} &= -\beta H_7^{n} -\alpha - \frac{20}{h^2} \end{align}

@ $$H_8$$:
\begin{align} \frac{H_7^{n+1} + H_9^{n+1} - \lambda H_8^{n+1} + H_{13}^{n+1} + 10}{h^2} &= -\beta H_8^{n} - \alpha\\ \implies\frac{H_7^{n+1} + H_9^{n+1} - \lambda H_8^{n+1} + H_{13}^{n+1}}{h^2} &= -\beta H_8^{n} -\alpha- \frac{10}{h^2} \end{align} This leads us to the system we must solve at each time-step: \begin{align} \left[ \matrix{-\lambda&1&0&1&0&0&0&0&0 \cr 1&-\lambda&1&0&1&0&0&0&0 \cr 0&1&-\lambda&0&0&1&0&0&0 \cr 1&0&0&-\lambda&1&0&1&0&0 \cr 0&1&0&1&-\lambda&1&0&1&0 \cr 0&0&1&0&1&-\lambda&0&0&1 \cr 0&0&0&1&0&0&-\lambda&1&0 \cr 0&0&0&0&1&0&1&-\lambda&1 \cr 0&0&0&0&0&1&0&1&-\lambda } \right] \left[ \matrix{H_7^{n+1} \cr H_8^{n+1} \cr H_9^{n+1} \cr H_{12}^{n+1} \cr H_{13}^{n+1} \cr H_{14}^{n+1} \cr H_{17}^{n+1} \cr H_{18}^{n+1} \cr H_{19}^{n+1}} \right] &= \left[ \matrix{ -\beta H_7^{n}-\alpha-\frac{20}{h^2} \cr -\beta H_8^{n}-\alpha - \frac{10}{h^2} \cr -\beta H_9^{n}-\alpha - \frac{20}{h^2} \cr -\beta H_{12}^{n}-\alpha - \frac{10}{h^2} \cr -\beta H_{13}^{n} \cr -\beta H_{14}^{n}-\alpha - \frac{10}{h^2} \cr -\beta H_{17}^{n}-\alpha - \frac{20}{h^2} \cr -\beta H_{18}^{n}-\alpha - \frac{10}{h^2} \cr -\beta H_{19}^{n}-\alpha - \frac{20}{h^2}} \right] \end{align} We can now move on to solving our problem in MATLAB

### Computation

First, we set up our space and time grids as well as our grid of indices. Note that we choose our mesh to be $$41\times41$$ and our period to be $$60$$ days.

 % Grids
% Indices
M = 41;
G = numgrid('S',M);
N = sum(G(:)>0);
center = ceil(M/2);
% Space
L =   4000;  h = L/(M-1);
[xg yg] = meshgrid(0:h:L,0:h:L);
% Time
Period = 60;
n = 1*Period;
k = Period/n;
tg = 0:k:Period;


Now we specify our parameters for the simulation:

 % Parameters
% Physical
edge = 10;           % Water table height at boundaries
S = .002;            % Specific yield
T = 300;             % Transmissivity
Q_well = -2000;      % Pumping rate
R = 0;               % Recharge away from well
R_well = Q_well/h^2; % Recharge at the well
% Discretization
lambda = 4 + ((h^2*S)/(T*k)); % Main diagonal entry of diff. matrix
beta = S/(T*k);
alpha = R/T;
alpha_well = R_well/T;


Instead of using the $$\verb* delsq*$$ as we did in the stead-state problem, we will build our own:

 e = ones(N,1);
jump = sqrt(N);
i = jump;
sub_sup = ones(N,1);
while i <=N
sub_sup(i) = 0; % sub/sup diagonal
i = i+jump;
end
Diff = 0*ones(N);
Diff = Diff + diag(-lambda*e);
Diff = Diff + diag(e(1:N-jump),-jump);
Diff = Diff + diag(e(1:N-jump),jump);
Diff = Diff + diag(sub_sup(1:N-1),1);
Diff = Diff + diag(sub_sup(1:N-1),-1);
Diff = Diff/h^2;


Once we have our differentiation matrix, we can proceed by specifying our initial head value before the pumping begins. We take it to be what we measured along the boundaries:

 % Build RHS/Initial Solution
rhs = edge*ones(N,1);


We are now ready to commence the simulation so we will start looping through the time-steps:

 % Solve at each time-step
for t = tg
rhs = -beta*rhs - alpha;
rhs(G(center,center)) = rhs(G(center,center)) - alpha_well; % Well
%
% Boundaries
rhs(G(  2,2:M-1)) = rhs(G(  2,2:M-1))-edge/h^2;%bottom
rhs(G(M-1,2:M-1)) = rhs(G(M-1,2:M-1))-edge/h^2;%top
rhs(G(2:M-1,  2)) = rhs(G(2:M-1,  2))-edge/h^2;%left
rhs(G(2:M-1,M-1)) = rhs(G(2:M-1,M-1))-edge/h^2;%right
%
% Solve
u = Diff\rhs;
U = G;
U(G>0) = full(u(G(G>0)));
U(1,1:M) = edge;  U(2:M-1,1) = edge;    % fill in boundary conditions
U(M,1:M) = edge;  U(2:M-1,M) = edge;
%
% Update
rhs = u;
end


Here is our solution which evolves over two months, with time-steps taken at each day. We reach a steady-state solution with the head at the well being $$h = 4.92m$$: File:Confined Aquifer 60Days.gif

### Verification

We do not have an analytic solution to compare to so verification is somewhat tricky. We will take our solutions with the finest grid-spacing as our 'analytic' solution to examine the convergence. In doing so, we run into a problem almost immediately. Since we will be comparing meshes of different grid spacing, we cannot perform a check of the norm of the difference of the matrices without interpolating out numerical solution to the same number of points as the analytic solution. Thus we have decided to use the difference of the norms of each matrix to approximate the error. Though this is not precise, it will give us an idea of whether or not we are converging. Time permitting, we will investigate interpolating our numerical solutions to the correct mesh size.

#### Space

We first fix our Period and time-step. The Period will be set as being $$5$$ days with a time step $$k=0.07$$. We take $$M=80$$ to be the spacing of our 'analytic' solution. We then use an increasing dimension size of $$[4, 8, 16, 32, 64]$$ and make the comparison of norms. After performing the computation, we get the following results about convergence.
700px
This is not what was expected! Though the method appears to be converging, it is converging at a rate much faster than our second order space difference should be. It can only be assumed that this was caused by our calculation of the error. Had we been able to interpolate our solution, we would have hoped to have seen a line of slope $$-2$$.

#### Time

We continue similarly to the previous error computation exept this time we fix our mesh spacing and decrease the size of our time steps. We set our period to be at the steady-state solution$$t=60$$ and we then compare increasing values of time points, $$[4,8,16,32,64]$$ with our 'analytic' solution having $$100$$ time points. We expect linear convergence as we are using simple Euler time-stepping. But we know that there will be some issues as we still have the problem of conflicting mesh spacing.
700px
Here, we can see that there is quadratic convergence. It is not what we expected but still an indication that we are converging.

## Transient Well Draw-down on an Unconfined Aquifer

Moving forward we now wish to solve draw-down problem on an unconfined aquifer governed by $$(3)$$. The big difference between this model and the previous models is the non-linear term on the right hand side of the equation. There are also some differences in parameters. Further, in this case we will be using the Crank-Nicholson method instead of the implicit Euler method and we will not use matrix techniques to solve the system exactly. Instead, we will use Gauss-Seidel iterations to solve the system.

### Analysis

First, we re-write $$(3)$$ as follows: \begin{align} \frac{ \partial^2 v }{ \partial x^2 } + \frac{ \partial^2 v }{ \partial y^2 } = \frac{S}{K\sqrt v}\frac{\partial v}{\partial t} -\frac{2R(x, y)}{K}. \end{align} We now discretize our problem: \begin{align} \frac{1}{2}\frac{V_{i+1,j}^{n+1} + V_{i-1,j}^{n+1}- 4V_{i,j}^{n+1} + V_{i,j-1}^{n+1} + V_{i,j+1}^{n+1} }{h^2} + \frac{1}{2}\frac{V_{i+1,j}^{n} + V_{i-1,j}^{n}- 4V_{i,j}^{n} + V_{i,j-1}^{n} + V_{i,j+1}^{n} }{h^2}&= \frac{S}{K\sqrt{V_{i,j}^{n}}}\frac{V_{i,j}^{n+1} - V_{i,j}^{n}}{k} -\frac{2R_{i,j}}{K}\\\\ \implies V_{i+1,j}^{n+1} + V_{i-1,j}^{n+1}- 4V_{i,j}^{n+1} + V_{i,j-1}^{n+1} + V_{i,j+1}^{n+1} + V_{i+1,j}^{n} + V_{i-1,j}^{n}- 4V_{i,j}^{n} + V_{i,j-1}^{n} + V_{i,j+1}^{n} &= \frac{2h^2 S}{Kk\sqrt{V_{i,j}^{n}}}(V_{i,j}^{n+1} - V_{i,j}^{n}) -\frac{4h^2R_{i,j}}{K}\\ \end{align} Now we let \begin{align} \alpha = \frac{2h^2S}{kK}, \beta_{i,j} = \frac{4h^2R_{i,j}}{K} \end{align} and we solve for $$V_{i,j}^{n+1}$$: \begin{align} V_{i+1,j}^{n+1} + V_{i-1,j}^{n+1}- 4V_{i,j}^{n+1} + V_{i,j-1}^{n+1} + V_{i,j+1}^{n+1} + V_{i+1,j}^{n} + V_{i-1,j}^{n}- 4V_{i,j}^{n} + V_{i,j-1}^{n} + V_{i,j+1}^{n} = \frac{\alpha}{\sqrt{V_{i,j}^{n}}}(V_{i,j}^{n+1} - V_{i,j}^{n}) -\beta_{i,j}\\ \end{align} \begin{align} \implies 4V_{i,j}^{n+1}+ \frac{\alpha}{\sqrt{V_{i,j}^{n}}}V_{i,j}^{n+1} &= \frac{\alpha V_{i,j}^{n}}{\sqrt{V_{i,j}^{n}}}+ V_{i+1,j}^{n+1} + V_{i-1,j}^{n+1}+ V_{i,j-1}^{n+1} + V_{i,j+1}^{n+1} + V_{i+1,j}^{n} + V_{i-1,j}^{n}- 4V_{i,j}^{n} + V_{i,j-1}^{n} + V_{i,j+1}^{n} +\beta_{i,j}\\ \implies V_{i,j}^{n+1} &= \frac{1}{4+\frac{\alpha}{\sqrt{V_{i,j}^{n}}}}(\frac{\alpha V_{i,j}^{n}}{\sqrt{V_{i,j}^{n}}}+ V_{i+1,j}^{n+1} + V_{i-1,j}^{n+1}+ V_{i,j-1}^{n+1} + V_{i,j+1}^{n+1} + V_{i+1,j}^{n} + V_{i-1,j}^{n}- 4V_{i,j}^{n} + V_{i,j-1}^{n} + V_{i,j+1}^{n} +\beta_{i,j}).\\ \end{align} We now want to solve our system so we proceed with Gauss-Seidel iteration, where $$m$$ is the iteration step: \begin{align} V_{i,j}^{n+1,m+1} = \frac{1}{4+\frac{\alpha}{\sqrt{V_{i,j}^{n,m}}}}(\frac{\alpha V_{i,j}^{n,m}}{\sqrt{V_{i,j}^{n,m}}}+ V_{i+1,j}^{n+1,m} + V_{i-1,j}^{n+1,m+1}+ V_{i,j-1}^{n+1,m+1} + V_{i,j+1}^{n+1,m} + V_{i+1,j}^{n,m} + V_{i-1,j}^{n,m}- 4V_{i,j}^{n,m} + V_{i,j-1}^{n,m} + V_{i,j+1}^{n,m} +\beta_{i,j}).\\ \end{align}

### Computation

As before, we start by defining our grid:

  % Grids
% Space
M       = 41;
center  = ceil(M/2);
L       = 4000;
h       = L/(M-1);
[xg yg] = meshgrid(0:h:L,0:h:L);
% Time
Period = 60;
n      = 1*Period;
k      = Period/n;
tg     = 0:k:Period;


We now specify our parameters. We will be using approximately the same values as we did in the confined aquifer but there will be a few changes. We will increase the storage coefficient by a factor of $$100$$ to account for dewatering and saturation of the aquifer. Note also that for unconfined aquifers, it is known as the specific yield. We also specify the recharge on our domain as being zero everywhere except at the pumping location. This simulates there being no rainfall.

  % Parameters
% Physical
edge             = 10;         % Water table height at boundaries
S                = .1;         % Specific yield
K                = 30;         % Transmissivity
Q_well           = -2000;      % Pumping rate
R                = 0*ones(M);  % Recharge away from well
R_well           = Q_well/h^2; % Recharge at the well
R(center,center) = R_well;
% Discretization
alpha = (2*S*h^2)/(K*k);
beta   = (4*h^2*R)/K;
% Iteration
TOL = .0000001;


For this simulation, again we will take the initial conditions to be the state of the aquifer before pumping has begun.

  % Set up initial conditions
V_ot     = edge^2*ones(M); % old timestep
V_nt     = 0*ones(M);      % new timestep
V_nt_oit = edge^2*ones(M); % old iteration


Now that we have set up our parameters we will go about solving the problem. This technique is very different from the matrix method we used on the confined aquifer. We have a series of nested loops, the outter most being our time-step, nested within this loop is our Gauss-Seidel iteration loop which inside we have to traverse our solution matrix. to get our desired accuracy, the Gauss-Seidel method takes $$14$$ iterations to converge.

  % Time-step loop
for t = tg
flag = true;
% Gauss-Seidel Iteration
while flag
V_nt(1,1:M)   =  edge^2; % boundary conditions
V_nt(M,1:M)   =  edge^2;
V_nt(2:M-1,1) =  edge^2;
V_nt(2:M-1,M) =  edge^2;
i = 2;
while i<=M-1
j = 2;
while j<= M-1
V_nt(i,j) = ((alpha*V_ot(i,j))/sqrt(V_ot(i,j))+ ...
V_ot(i-1,j) + V_ot(i+1,j) - 4*V_ot(i,j)    + ...
V_ot(i,j-1) + V_ot(i,j+1)                  + ...
V_nt(i-1,j) + V_nt_oit(i+1,j)              + ...
V_nt(i,j-1) + V_nt_oit(i,j+1) + beta(i,j)) / ...
(4 + alpha/sqrt(V_ot(i,j)));
j = j+1;
end
i = i+1;
end
% Compare iterations
if norm(abs(V_nt-V_nt_oit))<=TOL
flag = false;
end
% Update Gauss-Seidel iteration
V_nt_oit = V_nt;
end
% Update time-step
V_ot = V_nt;
end


Let's run the simulation for 60 days,:
File:Unconfined Aquifer 60Days.gif
After the 60 days have passed, the height of the water table at the well is given by $$h=5.24m$$.

### Verification

We proceed exactly as we did for the confined aquifer. Our results are almost identical.

#### Space

Fixing our period to be $$6$$ with $$180$$ time steps and using the same tolerance scheme as we did in the confined example, we produce the following error plot:
700px
Just as before, we can see the convergence though it is a crude estimate.

#### Time

We use $$M=100$$ for our grid spacing and we set our period at $$t=60$$. We steadily increase the number of time-steps as we did for the confined aquifer and we can see similar convergence to the previous example though the slope is not quite at $$m=-2$$.
700px

## Discussion

As it stands, these models give a rudimentary look into the field of groundwater modeling. The models are a long way from being something useful in an engineering process: field testing and calibration must be done to ensure that the models are accurately describing the systems in which they claim to. To calibrate the model, computed head values must be compared with field tested values, generally parameters are tweaked by trial and error until the simulation matches experiment. Further, these models don't take into account and heterogeneity in the aquifer nor do they have any interesting boundary conditions. Later models would readily incorporate varying boundary conditions but adding heterogeneity to the problems would involve a major overhaul of the code.
One observation made during the coding process was that directly solving the systems through iteration of the domain array proved to be astoundingly faster than using the matrix methods. Using matrix methods limited us to rather coarse mesh spacing whereas the iterative approach did not.

## Appendix

### Confined Aquifer scripts

Computation:
Main
Differentiation Matrix
Verification:
Space
Time

Computation:
Main
Verification:
Space
Time