% Document Type: Latex
% Master File: ground.tex
\section{Introduction}
Groundwater is a valuable and vital resource. Water covers over 70\% of the earth's surface and although less than 5\% of this water is defined to be groundwater, it is the largest source of freshwater which is vital to all life. The natural superior quality of groundwater as well as the reliability of its supply makes it a resource to be valued. Groundwater supplies must be studied and analyzed to efficiently use the resource without exploiting the source.

When evaluating a groundwater site for the construction of a possible well, the negative effects of the pumping must be weighed against the benefits of the yield of water from the well.  Drawdown or the decline in hydraulic head is an important consequence in the consideration of a well site as it can indicate a lowered water level that outweighs the benefit of the well; an unacceptable change to the hydrologic cycle; or an adverse effect on other wells tapping the same aquifer.

This paper first develops the theory behind the calculation of the drawdown cone or cone of depression in the potentiometric surface around a potential well.  The empirical law, Darcy's Law is discussed first to set the foundation for the analysis of groundwater flow.  Since groundwater flow is essentially three dimensional through the void spaces in the formations under the surface of the earth, Section~\ref{sec-three} derives the 3-dimensional form of Darcy's Law that forms the foundation for the drawdown equations.

A simple model of the hydraulic head drawdown as a mathematical boundary-value problem will be developed in Section~\ref{sec-flow}. In conclusion, Section~\ref{sec-theis} presents an analytical solution in the form of the Theis Equation in the specific example of the drawdown in a confined aquifer of a non-steady flow.

The appendices contain detailed discussions of topics mentioned throughout the paper as well as a glossary.  Code for a program written in Matlab that plots the drawdown versus the radial distance of a cone of depression is also included.  


\section{Darcy's Law}
In 1856 Henri Darcy published in an appendix to his book {\em Italiass Les Fontaines Publiques de la Ville de Dijon}, his results of an experiment he conducted for the city of Dijon, France. From the generalized results of his experiments we have the empiricism known as Darcy's Law which describes the rate of flow through filter sands and can be expressed as
\[q= -K*\frac{(h_2 -h_1)}{L},\]
where
\begin{tabbing}
$q$ \hspace{.4in}\= = volume of water crossing a unit cross sectional area per unit time  [LT$^{-1}$],\\ 
$K$ \>= factor of proportionality or hydraulic conductivity (Appendix~\ref{append:hydr})
 [LT$^{-1}$],\\
$h_2,h_1$ \>=\= water heights above a reference level measured by  manometers\\ 
\>\>which terminate above and below the sand column [L], and\\
$L$ \>= length of the flow path [L].\\ 
\end{tabbing}

Using an apparatus similar to the one in Figure~\ref{apparatus},
Darcy filled a circular tube with sand with an input tube at the top and a matching output tube at the bottom. 
\begin{figure}[hbt]
\centerline{\psfig{figure=appC.eps,height=3in}}
\caption{Experimental apparatus to illustrate Darcy's Law}
\label{apparatus}
\end{figure}
Manometers, tubes used to determine the elevation of water level or the measurement of the hydraulic head, are fitted into the cylinder at elevations $z_1$ and $z_2$. The elevations of the fluid inside the manometers are $h_1$ and $h_2$ with their difference in elevation equal to $\Delta h$ and the length between them $\Delta l$.  Water is directed into the input tube until all the pores are filled with water and consequently the inflow rate $Q$ is equal to the outflow rate. Let $A$ be a cross section of the cylinder with dimensions $L^2$.  Since the flow rate $Q$ has dimensions $\frac{L^3}{T}$, then $v$, the specific discharge throughout the cylinder defined as $v = \frac{Q}{A}$ has dimensions $\frac{L}{T}$. Note that we are calculating the specific discharge on a macroscopic scale which is easily measured as opposed to the microscopic view of the almost impossible measurement of the individual velocities between the grains of the porous medium. This type of ``sand-tank'' model is still being used today since Darcy's Law is still based on empirical evidence although numerous attempts have been made to derive his law from known physical laws.

Darcy's experiments showed conclusive evidence that the flow of water is directly related to  $\frac{(h_2 -h_1)}{L}$, the slope of the hydraulic gradient. $K$, although a constant of proportionality actually depends on the permeability of the soil as well as the fluid flowing through the soil. Darcy was able to show that $v$ is directly proportional to $h_1 - h_2$ when $\Delta l $ is held constant and inversely proportional to $\Delta l$ when $h_1 - h_2$ is held constant.~\cite{hand}\\ 
If we define  $\Delta h = h_2 - h_1$,  \hspace{.25in}  and with $v \propto -\Delta h,$  \hspace{.25in} and $v\propto \frac{1}{\Delta l}$,\\ 
we now have 
\[v = -K \frac{\Delta h}{\Delta l}.\] 
The negative in the equation implies that the groundwater flows in the direction of head loss.  Although Darcy's original experiments were conducted with the tube in a vertical position such that the flow was perpendicular to the variable $z = 0$, the law holds for groundwater flow in any direction even if the flow is forced upward from the datum $ z = 0 $ against gravitational forces.~\cite{cherry}

Although there are some theoretical and practical limitations to the use of Darcy's Law, it ``provides an accurate description of the flow of groundwater in almost all hydrogeological environments. In general, Darcy's Law holds for
\begin{enumerate}
\item saturated flow and for unsaturated flow,
\item steady-state flow and for transient flow,
\item flow in aquifers and for flow in aquitards,
\item flow in homogeneous systems in anisotropic media, and
\item flow in both rocks and granular media.~\cite{cherry}
\end{enumerate}
\subsection{Reynolds Number}
In 1883, Osborne Reynolds expanded the results of the experiments of Darcy and defined the ``critical velocity''.  Darcy's experiments were based on laminar flow, the movement of fluid essentially in parallel lines when viewed on a macroscopic scale. Critical velocity defines the transition between laminar flow and turbulent flow where fluid moves in irregular and rotational paths.
The number that describes the transition between the velocities is a dimensionless ratio named the Reynolds number and defines a range for the validity of Darcy's Law in turbulent flow. (Appendix~\ref{append:Rey})\\
Table~\ref{reytable} gives the ranges.
\begin{table}[hb]
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
Reynolds Number & &  \\ \hline \hline
& & \\ \hline
$Re \leq 10$ & Darcy's Law is valid & \\ \hline \hline 
& & \\ \hline
$10 < Re \leq 600-700$ & Partially turbulent flow exists & $\Rightarrow$ Darcy's Law does \underline{not} hold \\ \hline
\hline
& & \\ \hline
$Re \geq 600-700$ & Turbulent flow & $\Rightarrow$ Darcy's Law does \underline{not} hold \\ \hline
\end{tabular}
\end{center}
\caption{\label{reytable} \em Reynold's Number }
\end{table}

\subsection{Darcy's Law in Three Dimensions}
\label{sec-three}
The one dimensional form of Darcy's Law can be generalized to three dimensions such that the head potential, $\frac{\Delta h}{\Delta l}$ is expanded to be a function of the three space coordinates, $x, y,$ and $z$. Consequently the velocity $v$ is a vector with components $v_x, v_y$, and $v_z$ such that the potential head $h = h(x,y,z)$ is dependent on the position.

The three dimensional generalization of Darcy's Law, $v = -K \frac{\Delta h}{\Delta l}$,  can now be written as:
\[v_x = - K_x \frac{\partial h}{\partial x},\] 
\[v_y = - K_y \frac{\partial h}{\partial y},\] 
\[v_z = - K_z \frac{\partial h}{\partial z}.\]
Each component takes the partial derivative of $h$ with respect to the corresponding direction, since $h$ is dependent on $x, y$ and $z$.
\hspace{.5in}
\begin{figure}[hbt]
\centerline{\psfig{figure=cubeC.eps,height=4in}}
\caption{Elementary Cube}
\label{cube}
\end{figure}

Groundwater flow is a function of time and is essentially 3-dimensional in space.  The velocity vector at any point has components along three mutually perpendicular axes $x, y,$ and $z$ of a Cartesian coordinate system, thus the velocity is a function of $x, y, z$ and $t$.
%\cite
Consider the flow through a small elementary cube whose sides are of length $\delta x, \delta y$, and $\delta z$ (Figure~\ref{cube}). The law of conservation of matter holds for fluids. Specifically,
\begin{quote}
The net excess of mass flux, per unit of time, into or out of any infinitesimal volume element in the fluid system is exactly equal to the change per unit of time of the fluid density in the element multiplied by the free volume of the element.~\cite{hand}
\end{quote}
Therefore, the sum of the fluid entering the three faces of the cube is equal to the sum of the fluid leaving the cube from the opposite faces plus storage change.  The mass inflow in each of the $x, y, z$ directions is:
\[\mbox{mass inflow = fluid density $\times$ velocity  $\times$ cross-sectional area},\]
or for each of the specific directions:\\
mass inflow =
\[\rho v_x \delta y \delta z \hspace{1in} \mbox{$x$ direction},\]
\[\rho v_y \delta x \delta z \hspace{1in} \mbox{$y$ direction},\]
\[\rho v_z \delta x \delta y \hspace{1in} \mbox{$z$ direction},\]
where
\begin{tabbing}
$v_x, v_y, v_z$ \hspace{.6in} \= = fluid velocity (Darcian) in the $x, y, z$ directions respectively [LT$^{-1}$], and\\ 
$\rho$ \> = fluid density.
\end{tabbing}
Since the mass outflow is equal to the mass inflow plus storage change, we have:
\[\left(\rho v_x + \frac{\partial (\rho v_x)}{\partial x} \delta x \right) \delta y \delta z  \hspace{1in} \mbox{$x$ direction},\] \[\left(\rho v_y + \frac{\partial (\rho v_y)}{\partial y} \delta y \right) \delta x \delta z \hspace{1in} \mbox{$y$ direction}, \]
\[\left(\rho v_z + \frac{\partial (\rho v_z)}{\partial z} \delta z \right) \delta x \delta y \hspace{1in} \mbox{$z$ direction},\]
giving the total net inward flux as the sum of differences between the inflow and outflow in each of the three directions.
\[\mbox{Total net inward flux} = - \left( \frac{\partial (\rho v_x)}{\partial x} +\frac{\partial (\rho v_y)}{\partial y} +\frac{\partial (\rho v_z)}{\partial z} \right) \delta x \delta y \delta z. \]
Obeying the law of conservation of matter implies this total must equal the change in mass with respect to time:
\[- \left( \frac{\partial (\rho v_x)}{\partial x} +\frac{\partial (\rho v_y)}{\partial y} +\frac{\partial (\rho v_z)}{\partial z} \right) \delta x \delta y \delta z  = \frac{\partial(M)}{\partial t},\]
where
\begin{tabbing}
$M$ \hspace{.6in}\= = mass in the elementary cube,\\
\>=$\rho \theta \partial x \partial y \partial z $, and \\ 
$\theta$ \>= effective porosity, 
\end{tabbing}
giving us the {\em Equation of Continuity}:
\begin{equation}
\frac{(\partial \rho v_x)}{\partial x} + \frac{(\partial \rho v_y)}{\partial y} + \frac{(\partial \rho v_z)}{\partial z} = - \frac{\partial (\rho\theta)}{\partial t},
\label{eq:Cont}
\end{equation}
where
\begin{tabbing}
$\rho$ \hspace{.6in}\= = fluid density [ML$^{-3}$],\\ 
$v_x, v_y, v_z$ \>= fluid velocity (Darcian) in the $x, y, z$ directions respectively [LT$^{-1}$], and \\ 
$\theta$ \>= effective porosity. 
\end{tabbing}
Assuming density ($\rho$) constant and substituting the Darcian velocities(Appendix~\ref{append:Darcian})~\cite{hand}
\[v_s = - K \frac{\partial h}{\partial s},\]
with the specific directions $x, y, z$ in Equation~\ref{eq:Cont} we obtain\begin{equation}
K_x\frac{\partial ^2h}{\partial x^2} +K_y\frac{\partial ^2h}{\partial y^2}+K_z\frac{\partial ^2h}{\partial z^2} =S_s\frac{\partial h}{\partial t}. \label{eq:Kgov} 
\end{equation}

\section{Radial Flow to a Well}
\label{sec-flow}
Groundwater resources in a confined aquifer with a nonsteady-state flow can be evaluated for the consideration of the construction of wells. A confined aquifer is a primary source for well tapping and is defined to be an aquifer bounded by an upper and lower bed of material that allows only small amounts of groundwater to penetrate through. A nonsteady\-state flow is a flow in which the velocity changes direction or magnitude at some point in time.  Equation~\ref{eq:Kgov} is the {\em diffusion} equation. The solution $h(x, y, z, t)$ gives the value of the hydraulic head at any point in a flow field at any time. In this form, the equation describes a flow through a saturated anistrophic porous medium.  Incorportating the two basic assumptions of essentially horizontal flow obeying Darcy's Law which reduces the equation to a two dimensional form and that the flow is in an aquifer that is homogeneous, saturated and isotropic (i.e. $K = K_x = K_y = K_z$), Equation~\ref{eq:Kgov} reduces to:
\begin{equation}
\frac{\partial ^2h}{\partial x^2} +\frac{\partial ^2h}{\partial y^2}= \frac{S_s}{K} \frac{\partial h}{\partial t} = \frac{S}{T} \frac{\partial h}{\partial t}.
\label{eq:gov}
\end{equation}
It is assumed that the flow towards the well, where water is being removed, is radial. Hence, we convert Equation~\ref{eq:gov} into radial coordinates (Appendix~\ref{append:polar}) to describe the hydraulic head in terms of the drawdown around the well.
\begin{equation}
\frac{\partial^2s}{\partial r^2} + \frac{1}{r}\frac{\partial s}{\partial r} =\frac{1}{\kappa} \frac{\partial s}{\partial t},
\label{eq:radial}
\end{equation}
where
\begin{tabbing}
$s$ \hspace{.4in} \= = drawdown($h_0-h$) [L],\\ 
$r$ \>= radial distance from the well [L],\\
$t$ \>= time since pumping began [L],\\ 
$\kappa$ \>= hydraulic diffusivity T/S [L$^2$T$^{-1}$],\\ 
$T$ \>= transmissivity [L$^2$T$^{-1}$],\\
\>= function of Hydraulic Conductivity and the thickness of the saturated aquifer, and\\
$S$ \>= storativity\\ 
\>= function of the thickness of the saturated aquifer and the specific storativity. 
\end{tabbing}
\newpage
\begin{figure}[ht]
\centerline{\psfig{figure=conef.eps,height=4.5in}}
\caption{Cone of Depression in a Confined Aquifer}
\label{cone}
\end{figure}

Equation~\ref{eq:radial} provides solutions in terms of the drawdown based on the physics of the movement of the flow towards a well during pumping. The drawdown is computed in terms of the radial distance from the cone of depression which surrounds the pumping well. The cone of depression is the potentiometric or gravity free pressure surface and extends radially outward from the well and is a function of storativity, transmissivity, the pumping rate of the well and time. 
\begin{figure}[htb]
\centerline{\psfig{figure=radialC.eps,height=2.3in}}
\caption{Radial Flow from a Well}
\label{radial}
\end{figure}
Figure~\ref{radial} shows the radial flow from  $r = 0$ at the well to $r = \infty$ as the cone of depression radiates out from the well.

For the theoretical analysis of the hydraulic head drawdown in a proposed well, we will make the following assumptions:
\begin{itemize}
  \item A single pumping well in the aquifer,
  \item Pumping rate is constant with time,
  \item Well diameter is sufficiently small so that the volume of water removed from the well bore during pumping does not affect the total amount of water in the aquifer,
  \item  Well penetrates the entire thickness of the aquifer and receives water by a horizontal flow, and
  \item Hydraulic head or initial water level is horizontal throughout the aquifer prior to pumping~\cite{hand}. 
\end{itemize}
This simple model becomes a boundary-value problem as follows:\\
{\bf Governing Equation} (Equation~\ref{eq:radial})
\[\frac{\partial^2s}{\partial r^2} + \frac{1}{r}\frac{\partial s}{\partial r} =\frac{1}{\kappa} \frac{\partial s}{\partial t}.\]
{\bf Initial Condition}
\[s(r,0) = 0,\]
\hspace{1.5in}  (drawdown at any radial distance at time = 0).\\
\vspace{.3in}
{\bf Boundary Conditions} for $s(r,t)$
\[s(\infty,t) = 0,\]
\hspace{1.5in} (drawdown at an infinite distance for any time).
\[r_w\frac{\partial(s(r_w,t))}{\partial r} = - \frac{Q}{2\pi T},\] \hspace{2in} (discharge condition at the well).

\subsection{Theis Equation}
\label{sec-theis}
C.V. Theis first published in 1935 ``The Relation Between the Lowering of the Piezometric surface and the Rate and Duration of Discharge of a Well Using Groundwater Storage''.  He developed an analytic solution for the drawdown  for a non-steady flow in a confined aquifer.  Theis found the non-steady flow of groundwater to be analagous to the unsteady flow of heat in a homogeneous solid.  The Theis equation has become the most widely used equation in transient groundwater hydraulics and the solution in terms of drawdown is
\[s(r,t) = \frac{Q}{4 \pi T} W(u),\]
where
\begin{tabbing}
$s(r,t)$ \hspace{.4in}\= = drawdown at distance($r$) at time ($t$) after the start of pumping [L],\\
$Q$ \> = discharge rate [L$^3$T$^{-1}$], and\\
$W(u)$ \>= well function of Theis [Dimensionless]. 
\end{tabbing}

Specifically
\[W(u) = \int_{u}^{\infty} \frac{e^{-y}}{y} \partial y 
\hspace{3mm} = - \gamma - log_e u + u - \frac{u^2}{2 \cdot 2!} + \frac{u^3}{3 \cdot 3!} - \cdot \cdot \cdot \cdot+ (-1)^{n+1} \frac{u^n}{n \cdot n!} + \cdot \cdot \cdot \cdot, \] 
\[u= \frac{r^2S}{4Tt},\]
where\\
$\gamma$ \hspace{.4in} = Euler's constant = .577215664901532860606512.\\
Conventional Units
\[s(r,t) = \frac{114.6Q}{T} W(u),\]
\[u = \frac{1.87r^2S}{Tt},\]
where
\begin{tabbing}
$t$ \hspace{.4in}\= = time since pumping started(days),\\
$s(r,t)$ \> = drawdown [ft] at distance $r$ and time $t$ after pumping began,\\
$T$ \> = aquifer transmissivity\ [gpd/ft],\\
$S$ \> = aquifer storativity, and \\
$Q$ \> = discharge or well-pumping rate [gpm].
\end{tabbing}

The exponential integral is easily calculated with $u$ defined as above and is known as the {\em Well Function, W(u)}.  Storativity is the addition or release of water to the storage space due to the increase or decrease of hydraulic head, while transmissivity is a function of the hydraulic conductivity and the thickness of the aquifer and describes how easily the aquifer moves groundwater through its pore spaces. If $T, S,$ and the pumping rate $Q$ are known for the aquifer, the drawdown can be easily calculated.  A plot of the cone of depression can be calculated using values of $h_0 - h$ or the drawdown at various values of $r$ for a given time $t$. (See Appendix~\ref{append:plot})
``For a given aquifer the cone of depression increases in depth and extent with increasing time. Drawdown at any point at a given time is directly proportional to the pumping rate and inversely proportional to aquifer transmissivity and aquifer storativity.''\cite{cherry} Figure~\ref{examples} illustrates these relationships between high and low transmissivity and storativity. Low transmissivity produces a tight ``v'' shape, while high transmissivity pulls the cone out into a wider, more shallow shape.  High and low storativity both create a wide cone with the low storativity having a deeper v-shape.
\newpage
\begin{figure}[hb]
\centerline{\psfig{figure=exC.eps,height=4.5in}}
\caption{Effects of Storativity and Transmissivity on the Cone of Depression}
\label{examples}
\end{figure}
The calculation of the cone of depression for different pumping rates for various periods of time is just one small part of a many faceted and detailed evaluation of a potential well site. To ensure the performance and efficiency of the well and the protection of the aquifer being pumped from requires many tests and a careful study of all the factors invloved. Assumptions underlying all of the derivations must be taken into consideration also. Due to the nature of the assumptions in the model for drawdown using the Theis equation, it is most commonly used for single well analysis.
\nocite{seepage}
\nocite{eng}
\nocite{rec}
\nocite{analy}
\nocite{intro}
\nocite{con}
\nocite{south}
\nocite{theis}
\nocite{rel}
