Unconfined Flow
===============

One-dimensional, steady-state, constant-property, saturated flow in an unconfined aquifer
-----------------------------------------------------------------------------------------

Flow in an unconfined aquifer is inherently multi-dimensional, 
which generally precludes simple closed-form analytic solutions. 
In many practical applications however, vertical gradients and velocities are small relative
to their horizontal counterparts and can be neglected following Dupuit (1863) and Forchheimer (1930).
The Dupuit-Forchheimer theory of free-surface flow assumes that (Freeze and Cherry, 1979):
1) flow is horizontal and equipotential lines are vertical, and 
2) hydraulic gradient is equal to the slope of the free surface. With these assumptions and 
assuming the Cartesian :math:`x`- and :math:`y`-coordinates align with the 
principal axes of the hydraulic conductivity tensor,
the general expression for saturated flow in an unconfined aquifer becomes (e.g. de Marsily, 1986)

.. math:: 
   \frac{\partial}{\partial x} \left[\int_\sigma^h K_{xx}dz\frac{\partial h}{\partial x}\right]
   + \frac{\partial}{\partial y} \left[\int_\sigma^h K_{yy}dz\frac{\partial h}{\partial y}\right]
   = \omega_d \frac{\partial h}{\partial t} + Q
   :label: generalODE

where :math:`h` = hydraulic head (and elevation of the water table) [L],
:math:`t` = time [T],
:math:`K` = saturated hydraulic conductivity tensor [L/T], 
:math:`\sigma` = elevation of aquifer base [L], 
:math:`\omega_d` = specific yield or drainage porosity [-], and
:math:`Q` = volumetric flow rate per unit area withdrawn from the aquifer (sink) [L/T]. 
The hydraulic or total head (:math:`h`, [L]) is the sum of pressure head (:math:`P/\rho g`, [L]) 
and elevation (:math:`z`, [L])

.. math:: h = \frac{P}{\rho g}+z

where :math:`\rho` = density [M/L\ :sup:`3`\ ], and :math:`g` = gravitational acceleration [L/T\ :sup:`2`\ ]. 
For constant, isotropic, hydraulic conductivity, a horizontal aquifer base at 
:math:`\sigma = z = 0`, and one-dimensional flow in the :math:`x`-direction, Equation :eq:`generalODE` becomes simply

.. math:: \frac{d^2h^2}{dx^2} = \frac{2Q}{K}
   :label: ode

The ordinary differential equation :eq:`ode` is linear in :math:`h^2` and easily solved by 
direct integration as

.. math:: h^2 = C_1 x + C_2
   :label: generalSoln

where the integration constants :math:`C_1` and :math:`C_2` depend on the boundary conditions
and volumetric sink :math:`Q`.  
Following Aleman (2007, Section 4.2), two particular solutions are considered
for testing *Amanzi* implementation of prescribed hydraulic head, no-flow, and recharge boundary conditions, 
Darcy's law, and mass conservation on a problem involving variably saturated conditions and quasi 2D flow
(approximately horizontal in the saturated zone and approximately vertical in the unsaturated zone).


References
~~~~~~~~~~

Aleman, S. 2007. *PORFLOW Testing and Verification Document*. Savannah River National Laboratory technical report WSRC-STI-2007-00150 Rev 0. 193 p.

de Marsily, G. 1986. *Quantitative Hydrogeology: Groundwater Hydrology for Engineers*, Academic Press, Inc., Orlando Florida.

Dupuit, J. 1863. *Estudes Theoriques et Pratiques sur le mouvement des Eaux dans les canaux decouverts et travers les terrains permeables* (Second Edition ed.). Paris: Dunod.

Forchheimer, P. 1930. *Hydraulik*. 3d ed., Teubner, Leipzig, Berlin. 

Freeze, R.A. and J.A. Cherry. 1979. *Groundwater*. Prentice-Hall, Inc. Englewood Cliffs, NJ. 604 p.

Mualem, Y. 1976. *A new model predicting the hydraulic conductivity of unsaturated porous media*. Water Resour. Res. 12:513-522.

Phifer, M. A., M. R. Millings, and G. P. Flach. 2006. *Hydraulic Property Data Package for the E-Area and Z-Area Soils, 
Cementitious Materials, and Waste Zones*. Savannah River National Laboratory technical report WSRC-STI-2006-00198 Rev 0. 325 p.

Richards, L.A. 1931. *Capillary conduction of liquids through porous mediums*. Physics 1 (5): 318-333.

van Genuchten, M. Th. 1980. *A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils*. Soil Sci. Soc. Am. J. 44: 892-898.


Analytic solution: Prescribed head at boundaries
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

When hydraulic head is prescribed at both boundaries as

.. math:: 
   h(0) &= h_0\\
   h(L) &= h_L
   :label: BCs1

the analytic solution :eq:`generalSoln` for hydraulic head becomes

	.. math:: h^2 = h_0^2 + (h_L^2 - h_0^2) \frac{x}{L} + \frac{Q_{src}L^2}{K}\left( \frac{x}{L} \right) \left(1 - \frac{x}{L} \right)
		:label: specificSoln1

where :math:`L` = domain length [L], and :math:`Q_{src} \equiv -Q`. 


Amanzi verification test problem
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
To generate numerical results the following specifications are considered for 
analytic solution #1 (Aleman 2007, Figure 4.2.1):

.. image:: schematic/porflow_4.2.1.png
   :scale: 35 %
   :align: center

With the Dupuit approximation the analytic solution given by Equation :eq:`specificSoln1`
is one-dimensional in the horizontal coordinate and describes only the saturated zone. 
Because *Amanzi* does not directly solve a reduced governing equation set
equivalent to Equation :eq:`generalODE`, a two-dimensional :math:`(x,z)` simulation of the combined 
saturated and unsaturated zones using the Richards (1931) equation is required. 
Thus a vertical hydraulic conductivity and parameters defining moisture characteristic
curves for the unsaturated zone are required beyond the material properties implied by 
Equation :eq:`specificSoln1`. To minimize vertical gradients consistent with the Dupuit
assumption, the vertical hydraulic conductivity is set 10x higher than the
horizontal conductivity. To minimize non-vertical flow in the unsaturated zone
(and preserve the uniform distribution of recharge applied to the top of the model domain),
van Genuchten (1980) - Mualem (1976) parameters consistent with a gravel
are selected. Input parameters for the numerical simulation are summarized as:

* Domain (2D)

	* :math:`x_{min} = z_{min} = 0`
	* :math:`x_{max} = L = 100 ft,\ z_{max} = 60 ft`

* Boundary conditions

	* no-flow prescribed at the :math:`z_{min}, z_{max}` boundaries (:math:`Q_{src} = 0`)
	* prescribed hydraulic head at the x-coordinate boundaries: :math:`h(0) = 40 ft,\ h(L) = 20 ft`

* Material properties

	* :math:`\rho = 998.2 \: kg/m^3,\ \mu = 1.002 \cdot 10^{-3} \: Pa\cdot s,\ g = 9.807 \: m/s^2` 
	* :math:`K_{xx} = 10^{-3} ft/s`
	* :math:`K_{zz} = 10 K_{xx}` 
	* van Genuchten (1980) - Mualem (1976) parameters for a gravel based on Phifer et al. (2006):
		* :math:`\alpha = 0.143 cm^{-1} (1.46 \cdot 10^{-3} Pa^{-1})`
		* :math:`S_r = 0.052`
		* :math:`m = 0.314`

* Model discretization

	* :math:`\Delta x = 1 ft,\ \Delta z = 1 ft`

For these input specifications, *Amanzi* simulation output is expected to closely match

.. math:: h = \sqrt{1600 - 12 x}
   :label: expectedH1

from Equation :eq:`specificSoln1`. 


Amanzi verification test results
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The hydralic head computed using both Amanzi and the above analytic formulas
is shown on the figure and in the table below.

.. figure:: figures/hydraulic_head.png 
   :figclass: align-center
   :width: 600 px

    **Concentration transverse to the flow field from plume centerline at x=0.**


.. include:: table_values.txt



