groundwater_flow.tex

% Groundwater Flow Analysis Template
% Topics: Darcy's law, aquifer hydraulics, pumping tests, well hydraulics, numerical methods
% Style: Technical report with computational aquifer analysis

\documentclass[a4paper, 11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{subcaption}
\usepackage[makestderr]{pythontex}
\usepackage{geometry}
\geometry{margin=1in}
\usepackage{hyperref}
\usepackage{float}

% Theorem environments
\newtheorem{definition}{Definition}[section]
\newtheorem{theorem}{Theorem}[section]
\newtheorem{example}{Example}[section]
\newtheorem{remark}{Remark}[section]

\title{Groundwater Flow Analysis: Darcy's Law and Well Hydraulics}
\author{Hydrogeology Research Group}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This technical report presents comprehensive computational analysis of groundwater flow in
porous media, focusing on aquifer hydraulics and well dynamics. We examine the fundamental
principles of Darcy's law, develop analytical solutions for radial flow to wells using the
Theis equation and Cooper-Jacob approximation, and implement numerical finite difference
methods for solving the 2D groundwater flow equation. Pumping test analysis demonstrates
determination of hydraulic parameters including transmissivity ($T = 1.2 \times 10^{-3}$ m²/s)
and storativity ($S = 2.5 \times 10^{-4}$), with visualization of hydraulic head distributions
and groundwater flow fields under various pumping scenarios.
\end{abstract}

\section{Introduction}

Groundwater flow through porous media is governed by Darcy's law, which relates hydraulic
gradient to discharge velocity through the hydraulic conductivity of the formation. Understanding
these flow dynamics is essential for water resource management, contaminant transport modeling,
and wellfield design. The mathematical framework combines conservation of mass with Darcy's
empirical relationship to yield partial differential equations describing hydraulic head
distributions in both steady-state and transient conditions.

\begin{definition}[Darcy's Law]
For one-dimensional flow through a porous medium, the specific discharge $q$ is proportional
to the hydraulic gradient:
\begin{equation}
q = -K \frac{dh}{dx}
\end{equation}
where $K$ is hydraulic conductivity (m/s) and $h$ is hydraulic head (m). The negative sign
indicates flow from high to low head.
\end{definition}

\section{Theoretical Framework}

\subsection{Governing Equation for Groundwater Flow}

Combining Darcy's law with the continuity equation for incompressible flow yields the governing
equation for transient groundwater flow in a confined aquifer. Consider a representative
elementary volume in an aquifer where water is released from storage as head declines.

\begin{theorem}[2D Groundwater Flow Equation]
For horizontal flow in a confined aquifer with constant transmissivity, the governing equation is:
\begin{equation}
\frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial y^2} = \frac{S}{T} \frac{\partial h}{\partial t}
\end{equation}
where $T = Kb$ is transmissivity (m²/s), $b$ is aquifer thickness (m), and $S$ is storativity
(dimensionless). For steady-state conditions, this reduces to Laplace's equation.
\end{theorem}

\subsection{Radial Flow to a Well}

When a well pumps water from a confined aquifer, flow converges radially toward the well
creating a cone of depression. The Theis solution provides the fundamental analytical framework
for analyzing this transient drawdown response.

\begin{definition}[Theis Equation]
The drawdown $s$ at radial distance $r$ from a pumping well at time $t$ is:
\begin{equation}
s(r, t) = \frac{Q}{4\pi T} W(u)
\end{equation}
where $Q$ is pumping rate (m³/s), $W(u)$ is the well function (exponential integral), and
$u = r^2 S / (4Tt)$ is the dimensionless time parameter \cite{theis1935relation}.
\end{definition}

The well function is defined as:
\begin{equation}
W(u) = \int_u^\infty \frac{e^{-y}}{y} dy = -0.5772 - \ln(u) + u - \frac{u^2}{2 \cdot 2!} + \frac{u^3}{3 \cdot 3!} - \cdots
\end{equation}

\begin{theorem}[Cooper-Jacob Approximation]
For small values of $u$ (large $t$ or small $r$), the series expansion truncates to:
\begin{equation}
s \approx \frac{2.30Q}{4\pi T} \log_{10}\left(\frac{2.25Tt}{r^2 S}\right)
\end{equation}
This provides a straight-line relationship on semi-log plots used for pumping test analysis
\cite{cooper1946drawdown}.
\end{theorem}

\section{Computational Implementation}

We now implement computational methods for analyzing groundwater flow, including analytical
solutions for well hydraulics and numerical finite difference approximations for complex
boundary conditions. The analysis begins with establishing fundamental hydraulic properties
and developing functions for the Theis well function calculation.

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import expi, exp1
from scipy.optimize import curve_fit
from scipy.integrate import odeint
import matplotlib.patches as patches

np.random.seed(42)

# Configure matplotlib for LaTeX rendering
plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.size'] = 10

# Define hydraulic parameters for confined aquifer
K = 1.5e-5  # Hydraulic conductivity (m/s)
b = 25.0    # Aquifer thickness (m)
T = K * b   # Transmissivity (m^2/s)
S = 2.5e-4  # Storativity (dimensionless)
Q = 0.015   # Pumping rate (m^3/s = 15 L/s)

def well_function(u):
    """
    Theis well function W(u) = exponential integral.
    For small u, use series expansion for numerical stability.
    """
    u = np.asarray(u)
    w = np.zeros_like(u, dtype=float)

    # For u < 0.01, use series expansion
    small_u = u < 0.01
    if np.any(small_u):
        u_small = u[small_u]
        w[small_u] = -0.5772 - np.log(u_small) + u_small - u_small**2/(2*2) + u_small**3/(3*6)

    # For larger u, use exponential integral
    large_u = ~small_u
    if np.any(large_u):
        w[large_u] = -expi(-u[large_u])

    return w

def theis_drawdown(r, t, Q, T, S):
    """
    Calculate drawdown using Theis equation.

    Parameters:
    r: radial distance from well (m)
    t: time since pumping started (s)
    Q: pumping rate (m^3/s)
    T: transmissivity (m^2/s)
    S: storativity (dimensionless)

    Returns:
    s: drawdown (m)
    """
    u = r**2 * S / (4 * T * t)
    w = well_function(u)
    s = Q / (4 * np.pi * T) * w
    return s

def cooper_jacob_drawdown(r, t, Q, T, S):
    """
    Cooper-Jacob approximation for large t, small r.
    Valid when u < 0.01.
    """
    s = 2.30 * Q / (4 * np.pi * T) * np.log10(2.25 * T * t / (r**2 * S))
    return s

# Radial distances for analysis
r_obs = np.array([10, 25, 50, 100, 200, 300])  # observation wells (m)
time_pumping = np.logspace(1, 5, 100)  # 10 s to ~1 day

# Calculate drawdown at multiple observation points
drawdown_curves = {}
for r in r_obs:
    drawdown_curves[r] = theis_drawdown(r, time_pumping, Q, T, S)
\end{pycode}

\subsection{Drawdown Analysis at Multiple Observation Wells}

Pumping test analysis relies on measuring drawdown at multiple observation wells located at
varying radial distances from the pumping well. These time-drawdown curves reveal the transient
response of the aquifer and enable determination of hydraulic properties through curve matching
or straight-line methods. The following analysis examines drawdown evolution over time periods
ranging from minutes to over 24 hours of continuous pumping.

\begin{pycode}
# Create drawdown vs time plots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

# Plot 1: Time-drawdown curves (semi-log)
colors = plt.cm.viridis(np.linspace(0, 0.9, len(r_obs)))
for i, r in enumerate(r_obs):
    ax1.semilogx(time_pumping / 3600, drawdown_curves[r],
                 linewidth=2.5, color=colors[i], label=f'$r = {r}$ m')
ax1.set_xlabel('Time (hours)', fontsize=11)
ax1.set_ylabel('Drawdown (m)', fontsize=11)
ax1.set_title('Time-Drawdown Response at Observation Wells', fontsize=12, fontweight='bold')
ax1.legend(loc='lower right', fontsize=9)
ax1.grid(True, alpha=0.3, which='both')
ax1.set_xlim(time_pumping[0]/3600, time_pumping[-1]/3600)

# Plot 2: Distance-drawdown at fixed times
time_snapshots = np.array([1, 6, 24]) * 3600  # 1, 6, 24 hours in seconds
r_fine = np.logspace(0.5, 2.8, 200)  # 3 m to 600 m

colors_time = ['red', 'blue', 'green']
for i, t_snap in enumerate(time_snapshots):
    s_snap = theis_drawdown(r_fine, t_snap, Q, T, S)
    ax2.semilogx(r_fine, s_snap, linewidth=2.5, color=colors_time[i],
                 label=f'$t = {int(t_snap/3600)}$ hr')
ax2.set_xlabel('Distance from well (m)', fontsize=11)
ax2.set_ylabel('Drawdown (m)', fontsize=11)
ax2.set_title('Radial Drawdown Distribution', fontsize=12, fontweight='bold')
ax2.legend(loc='upper right', fontsize=9)
ax2.grid(True, alpha=0.3, which='both')
ax2.invert_yaxis()

plt.tight_layout()
plt.savefig('groundwater_flow_drawdown.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{groundwater_flow_drawdown.pdf}
\caption{Drawdown analysis for pumping test in confined aquifer with $T = 3.75 \times 10^{-4}$ m²/s
and $S = 2.5 \times 10^{-4}$. Left panel shows time-drawdown curves at six observation wells
(10--300 m from pumping well), demonstrating the characteristic logarithmic increase in drawdown
with time. Right panel displays radial drawdown profiles at 1, 6, and 24 hours, illustrating
the expanding cone of depression. Maximum drawdown at the 10-m observation well reaches
approximately 2.8 m after 24 hours of pumping at 15 L/s. The semi-logarithmic plotting reveals
the straight-line behavior predicted by the Cooper-Jacob approximation at late times.}
\end{figure}

\section{Cooper-Jacob Analysis for Parameter Estimation}

The Cooper-Jacob method provides a practical approach for determining transmissivity and
storativity from pumping test data without requiring curve matching to type curves. By plotting
drawdown versus the logarithm of time, the late-time data forms a straight line whose slope
and intercept yield the aquifer parameters. This graphical method is widely used in field
hydrogeology for its simplicity and computational efficiency.

\begin{pycode}
# Generate synthetic pumping test data with measurement noise
t_measured = np.array([1, 2, 3, 5, 8, 10, 15, 20, 30, 45, 60, 90, 120, 180,
                       300, 480, 720, 1080, 1440, 2160, 4320, 8640, 14400]) * 60  # minutes to seconds
r_test_well = 50  # observation well at 50 m
s_measured = theis_drawdown(r_test_well, t_measured, Q, T, S)
# Add realistic measurement noise
noise_std = 0.02  # 2 cm standard deviation
s_measured_noisy = s_measured + np.random.normal(0, noise_std, len(s_measured))

# Cooper-Jacob analysis: plot s vs log10(t)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

# Semi-log plot for Cooper-Jacob
ax1.plot(t_measured / 3600, s_measured_noisy, 'bo', markersize=8,
         markeredgecolor='black', markeredgewidth=1, label='Field measurements')

# Fit straight line to late-time data (Cooper-Jacob approximation valid for u < 0.01)
u_values = r_test_well**2 * S / (4 * T * t_measured)
late_time_idx = u_values < 0.01
t_late = t_measured[late_time_idx]
s_late = s_measured_noisy[late_time_idx]

# Linear fit to log(t) vs s
log10_t_late = np.log10(t_late)
coeffs = np.polyfit(log10_t_late, s_late, 1)
slope_cj = coeffs[0]  # m per log cycle
intercept_cj = coeffs[1]

# Calculate T and S from slope and intercept
T_estimated = 2.30 * Q / (4 * np.pi * slope_cj)
# From intercept: t_0 is where s = 0
t_0 = 10**(-intercept_cj / slope_cj)
S_estimated = 2.25 * T_estimated * t_0 / r_test_well**2

# Plot fitted line
t_fit = np.logspace(np.log10(t_late[0]), np.log10(t_measured[-1]), 50)
s_fit = slope_cj * np.log10(t_fit) + intercept_cj
ax1.semilogx(t_fit / 3600, s_fit, 'r-', linewidth=2.5, label='Cooper-Jacob fit')
ax1.set_xlabel('Time (hours)', fontsize=11)
ax1.set_ylabel('Drawdown (m)', fontsize=11)
ax1.set_title('Cooper-Jacob Semi-Log Analysis', fontsize=12, fontweight='bold')
ax1.legend(loc='lower right', fontsize=9)
ax1.grid(True, alpha=0.3, which='both')
ax1.text(0.05, 0.95, f'$\\Delta s$ per log cycle = {slope_cj:.3f} m\\n$T$ = {T_estimated:.2e} m$^2$/s\\n$S$ = {S_estimated:.2e}',
         transform=ax1.transAxes, fontsize=9, verticalalignment='top',
         bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))

# Residual analysis
s_predicted = slope_cj * np.log10(t_measured) + intercept_cj
residuals = s_measured_noisy - s_predicted
ax2.plot(t_measured / 3600, residuals * 100, 'go', markersize=8,
         markeredgecolor='black', markeredgewidth=1)
ax2.axhline(y=0, color='red', linestyle='--', linewidth=2)
ax2.axhline(y=noise_std*100, color='gray', linestyle=':', linewidth=1.5, alpha=0.7)
ax2.axhline(y=-noise_std*100, color='gray', linestyle=':', linewidth=1.5, alpha=0.7)
ax2.set_xlabel('Time (hours)', fontsize=11)
ax2.set_ylabel('Residual (cm)', fontsize=11)
ax2.set_title('Cooper-Jacob Fit Residuals', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.set_xscale('log')

plt.tight_layout()
plt.savefig('groundwater_flow_cooper_jacob.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{groundwater_flow_cooper_jacob.pdf}
\caption{Cooper-Jacob straight-line analysis of pumping test data from observation well at 50 m
distance. Left panel shows drawdown measurements (blue circles) plotted versus logarithm of time,
with the fitted straight line (red) representing the Cooper-Jacob approximation. The slope of
\py{f'{slope_cj:.3f}'} m per log cycle yields transmissivity $T = \py{f'{T_estimated:.2e}'}$ m²/s,
comparing favorably with the true value of \py{f'{T:.2e}'} m²/s. Right panel displays fit residuals
(green circles) showing random scatter around zero within $\pm$2 cm (gray dashed lines), confirming
good agreement between measured and predicted drawdown. The early-time deviations reflect the
breakdown of the Cooper-Jacob approximation when $u > 0.01$.}
\end{figure}

\section{Two-Dimensional Hydraulic Head Distribution}

Numerical solution of the 2D groundwater flow equation enables visualization of hydraulic head
distributions and flow patterns under complex boundary conditions. Using finite difference
approximations on a discrete grid, we solve for steady-state head distributions in a confined
aquifer with a pumping well and specified constant-head boundaries. The resulting head contours
and flow vectors reveal the convergent radial flow pattern characteristic of well hydraulics.

\begin{pycode}
# Solve 2D steady-state groundwater flow equation using finite differences
# Domain: 1000 m x 1000 m confined aquifer with central pumping well

nx, ny = 101, 101  # Grid resolution
Lx, Ly = 1000, 1000  # Domain size (m)
dx = Lx / (nx - 1)
dy = Ly / (ny - 1)

x = np.linspace(0, Lx, nx)
y = np.linspace(0, Ly, ny)
X, Y = np.meshgrid(x, y)

# Initial head (m above datum)
h0 = 100.0
h = np.ones((ny, nx)) * h0

# Boundary conditions: constant head on all boundaries
h[0, :] = h0      # bottom
h[-1, :] = h0     # top
h[:, 0] = h0      # left
h[:, -1] = h0     # right

# Well location at center
i_well, j_well = ny // 2, nx // 2
r_well = 0.15  # well radius (m)

# Effective well cell area
A_well = np.pi * r_well**2
# Pumping rate creates sink term
well_flux = Q / (dx * dy)  # m/s

# Iterative solver: Gauss-Seidel with successive over-relaxation (SOR)
omega = 1.5  # relaxation parameter
max_iter = 5000
tolerance = 1e-6

for iteration in range(max_iter):
    h_old = h.copy()

    for i in range(1, ny - 1):
        for j in range(1, nx - 1):
            # Standard finite difference stencil
            h_new = 0.25 * (h[i+1, j] + h[i-1, j] + h[i, j+1] + h[i, j-1])

            # Apply well sink if at well location
            if i == i_well and j == j_well:
                # Modify equation for pumping well
                source_term = -well_flux * dx * dy / T
                h_new = 0.25 * (h[i+1, j] + h[i-1, j] + h[i, j+1] + h[i, j-1] - source_term)

            # Apply SOR
            h[i, j] = omega * h_new + (1 - omega) * h[i, j]

    # Check convergence
    max_change = np.max(np.abs(h - h_old))
    if max_change < tolerance:
        break

# Calculate drawdown
drawdown = h0 - h

# Calculate hydraulic gradients and specific discharge
dh_dx = np.zeros_like(h)
dh_dy = np.zeros_like(h)
dh_dx[:, 1:-1] = (h[:, 2:] - h[:, :-2]) / (2 * dx)
dh_dy[1:-1, :] = (h[2:, :] - h[:-2, :]) / (2 * dy)

# Darcy velocity (specific discharge)
qx = -K * dh_dx
qy = -K * dh_dy
q_magnitude = np.sqrt(qx**2 + qy**2)

# Create visualization
fig = plt.figure(figsize=(14, 12))

# Plot 1: Hydraulic head contours
ax1 = plt.subplot(2, 2, 1)
levels_head = np.linspace(h.min(), h0, 25)
cs1 = ax1.contourf(X, Y, h, levels=levels_head, cmap='Blues_r')
cs1_lines = ax1.contour(X, Y, h, levels=15, colors='black', linewidths=0.8, alpha=0.4)
ax1.clabel(cs1_lines, inline=True, fontsize=7, fmt='%0.1f')
cbar1 = plt.colorbar(cs1, ax=ax1, label='Hydraulic head (m)')
ax1.plot(x[j_well], y[i_well], 'r*', markersize=20, markeredgecolor='black', markeredgewidth=1.5,
         label='Pumping well')
ax1.set_xlabel('x (m)', fontsize=11)
ax1.set_ylabel('y (m)', fontsize=11)
ax1.set_title('Steady-State Hydraulic Head', fontsize=12, fontweight='bold')
ax1.legend(loc='upper right', fontsize=9)
ax1.set_aspect('equal')

# Plot 2: Drawdown contours
ax2 = plt.subplot(2, 2, 2)
levels_dd = np.linspace(0, drawdown.max(), 25)
cs2 = ax2.contourf(X, Y, drawdown, levels=levels_dd, cmap='Reds')
cs2_lines = ax2.contour(X, Y, drawdown, levels=12, colors='black', linewidths=0.8, alpha=0.4)
ax2.clabel(cs2_lines, inline=True, fontsize=7, fmt='%0.2f')
cbar2 = plt.colorbar(cs2, ax=ax2, label='Drawdown (m)')
ax2.plot(x[j_well], y[i_well], 'k*', markersize=20, markeredgecolor='white', markeredgewidth=1.5)
ax2.set_xlabel('x (m)', fontsize=11)
ax2.set_ylabel('y (m)', fontsize=11)
ax2.set_title('Drawdown Cone of Depression', fontsize=12, fontweight='bold')
ax2.set_aspect('equal')

# Plot 3: Flow field (quiver plot)
ax3 = plt.subplot(2, 2, 3)
skip = 5
ax3.contourf(X, Y, drawdown, levels=levels_dd, cmap='Reds', alpha=0.3)
Q_plot = ax3.quiver(X[::skip, ::skip], Y[::skip, ::skip],
                     qx[::skip, ::skip], qy[::skip, ::skip],
                     q_magnitude[::skip, ::skip], cmap='viridis', scale=5e-6)
cbar3 = plt.colorbar(Q_plot, ax=ax3, label='Specific discharge (m/s)')
ax3.plot(x[j_well], y[i_well], 'r*', markersize=20, markeredgecolor='black', markeredgewidth=1.5)
ax3.set_xlabel('x (m)', fontsize=11)
ax3.set_ylabel('y (m)', fontsize=11)
ax3.set_title('Groundwater Flow Field', fontsize=12, fontweight='bold')
ax3.set_aspect('equal')

# Plot 4: Radial profile comparison
ax4 = plt.subplot(2, 2, 4)
# Extract drawdown along centerline
centerline_dd = drawdown[i_well, j_well:]
r_centerline = x[j_well:] - x[j_well]
# Analytical steady-state solution (Thiem equation for confined aquifer)
r_analytical = np.linspace(r_well, 500, 200)
# Thiem: s = Q/(2*pi*T) * ln(R/r) where R is radius of influence
R_influence = 500  # m
s_analytical = Q / (2 * np.pi * T) * np.log(R_influence / r_analytical)

ax4.plot(r_centerline[1:], centerline_dd[1:], 'b-', linewidth=2.5, label='Numerical solution')
ax4.plot(r_analytical, s_analytical, 'r--', linewidth=2, label='Thiem equation')
ax4.set_xlabel('Radial distance from well (m)', fontsize=11)
ax4.set_ylabel('Drawdown (m)', fontsize=11)
ax4.set_title('Radial Drawdown Profile', fontsize=12, fontweight='bold')
ax4.legend(loc='upper right', fontsize=9)
ax4.grid(True, alpha=0.3)
ax4.set_xlim(0, 500)

plt.tight_layout()
plt.savefig('groundwater_flow_2d_field.pdf', dpi=150, bbox_inches='tight')
plt.close()

# Store key results
max_drawdown = drawdown.max()
drawdown_at_100m = drawdown[i_well, j_well + int(100/dx)]
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{groundwater_flow_2d_field.pdf}
\caption{Two-dimensional numerical solution of steady-state groundwater flow in 1000 m × 1000 m
confined aquifer with central pumping well (red star). Top-left panel shows hydraulic head
distribution with contour lines indicating head decline from 100 m at boundaries to
\py{f'{h[i_well, j_well]:.2f}'} m at the well. Top-right panel displays the cone of depression
with maximum drawdown of \py{f'{max_drawdown:.2f}'} m at the well location. Bottom-left panel
presents the convergent flow field as velocity vectors colored by specific discharge magnitude,
clearly showing radial flow toward the pumping well. Bottom-right panel compares the numerical
solution (blue line) with the analytical Thiem equation (red dashed), demonstrating excellent
agreement beyond the near-well region where numerical discretization effects become negligible.
At 100 m distance, drawdown is \py{f'{drawdown_at_100m:.3f}'} m.}
\end{figure}

\section{Pumping Test Parameter Sensitivity}

Understanding the sensitivity of drawdown to hydraulic parameters is crucial for interpreting
pumping test results and quantifying uncertainty in parameter estimates. We examine how variations
in transmissivity and storativity affect the time-drawdown response, revealing which parameters
dominate at different time scales. Early-time drawdown is primarily controlled by storativity,
while late-time response depends mainly on transmissivity.

\begin{pycode}
# Sensitivity analysis: vary T and S independently
r_sensitivity = 50  # observation well distance

# Transmissivity variation (S constant)
T_range = np.array([0.5, 0.75, 1.0, 1.5, 2.0]) * T
colors_T = plt.cm.Reds(np.linspace(0.4, 0.9, len(T_range)))

# Storativity variation (T constant)
S_range = np.array([0.5, 0.75, 1.0, 1.5, 2.0]) * S
colors_S = plt.cm.Blues(np.linspace(0.4, 0.9, len(S_range)))

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

# Panel 1: Transmissivity sensitivity
for i, T_var in enumerate(T_range):
    s_T = theis_drawdown(r_sensitivity, time_pumping, Q, T_var, S)
    label_T = f'$T = {T_var/T:.2f} T_0$'
    ax1.semilogx(time_pumping / 3600, s_T, linewidth=2.5, color=colors_T[i], label=label_T)
ax1.set_xlabel('Time (hours)', fontsize=11)
ax1.set_ylabel('Drawdown (m)', fontsize=11)
ax1.set_title('Transmissivity Sensitivity ($S$ constant)', fontsize=12, fontweight='bold')
ax1.legend(loc='lower right', fontsize=9)
ax1.grid(True, alpha=0.3, which='both')

# Panel 2: Storativity sensitivity
for i, S_var in enumerate(S_range):
    s_S = theis_drawdown(r_sensitivity, time_pumping, Q, T, S_var)
    label_S = f'$S = {S_var/S:.2f} S_0$'
    ax2.semilogx(time_pumping / 3600, s_S, linewidth=2.5, color=colors_S[i], label=label_S)
ax2.set_xlabel('Time (hours)', fontsize=11)
ax2.set_ylabel('Drawdown (m)', fontsize=11)
ax2.set_title('Storativity Sensitivity ($T$ constant)', fontsize=12, fontweight='bold')
ax2.legend(loc='lower right', fontsize=9)
ax2.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('groundwater_flow_sensitivity.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{groundwater_flow_sensitivity.pdf}
\caption{Sensitivity of time-drawdown response to hydraulic parameters at observation well 50 m
from pumping well. Left panel shows the effect of varying transmissivity from 0.5$T_0$ to 2.0$T_0$
while holding storativity constant, demonstrating that higher transmissivity reduces drawdown
at all times by increasing the aquifer's ability to transmit water. Right panel illustrates
storativity variation from 0.5$S_0$ to 2.0$S_0$ with constant transmissivity, showing that
storativity primarily affects early-time response and the timing of drawdown propagation. Higher
storativity delays the arrival of the pressure transient and reduces early drawdown by providing
more water from storage. At late times ($t > 10$ hours), all storativity curves converge to
parallel straight lines with identical slopes, confirming that late-time Cooper-Jacob analysis
depends only on transmissivity.}
\end{figure}

\section{Specific Capacity and Well Performance}

Well performance is commonly characterized by specific capacity, defined as the pumping rate
divided by drawdown at the well. This metric provides a practical measure of well productivity
that depends on both aquifer properties (transmissivity) and well construction (radius, screen
efficiency). Understanding the relationship between specific capacity, pumping rate, and aquifer
parameters enables optimal wellfield design and prediction of long-term well performance.

\begin{pycode}
# Well performance analysis
Q_range = np.linspace(0.005, 0.03, 10)  # 5 to 30 L/s
r_w = 0.15  # well radius
t_performance = 24 * 3600  # drawdown after 24 hours

# Calculate drawdown at well for different pumping rates
s_well = np.zeros(len(Q_range))
for i, Q_i in enumerate(Q_range):
    s_well[i] = theis_drawdown(r_w, t_performance, Q_i, T, S)

# Specific capacity
specific_capacity = Q_range / s_well  # m^3/s per m = m^2/s

# Well efficiency (actual vs theoretical)
# Theoretical drawdown (no well losses)
# Actual includes well losses: s_actual = s_formation + B*Q + C*Q^2
well_loss_coeff = 50  # s^2/m^5
s_well_loss = well_loss_coeff * Q_range**2
s_actual = s_well + s_well_loss
efficiency = s_well / s_actual * 100

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(13, 10))

# Plot 1: Drawdown vs pumping rate
ax1.plot(Q_range * 1000, s_well, 'b-', linewidth=2.5, label='Formation loss')
ax1.plot(Q_range * 1000, s_actual, 'r--', linewidth=2.5, label='Total (with well loss)')
ax1.fill_between(Q_range * 1000, s_well, s_actual, alpha=0.3, color='orange', label='Well loss')
ax1.set_xlabel('Pumping rate (L/s)', fontsize=11)
ax1.set_ylabel('Drawdown at well (m)', fontsize=11)
ax1.set_title('Well Drawdown vs Pumping Rate', fontsize=12, fontweight='bold')
ax1.legend(loc='upper left', fontsize=9)
ax1.grid(True, alpha=0.3)

# Plot 2: Specific capacity
ax2.plot(Q_range * 1000, specific_capacity * 1000, 'g-', linewidth=2.5, marker='o', markersize=8)
ax2.set_xlabel('Pumping rate (L/s)', fontsize=11)
ax2.set_ylabel('Specific capacity (L/s per m)', fontsize=11)
ax2.set_title('Specific Capacity Decline', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

# Plot 3: Well efficiency
ax3.plot(Q_range * 1000, efficiency, 'm-', linewidth=2.5, marker='s', markersize=8)
ax3.axhline(y=70, color='red', linestyle='--', linewidth=2, alpha=0.7, label='70\\% threshold')
ax3.set_xlabel('Pumping rate (L/s)', fontsize=11)
ax3.set_ylabel('Well efficiency (\\%)', fontsize=11)
ax3.set_title('Well Efficiency vs Pumping Rate', fontsize=12, fontweight='bold')
ax3.legend(loc='upper right', fontsize=9)
ax3.grid(True, alpha=0.3)
ax3.set_ylim(50, 100)

# Plot 4: Step-drawdown test simulation
Q_steps = np.array([0.01, 0.015, 0.02, 0.025]) * 1000  # L/s
step_duration = 3 * 3600  # 3 hours per step
t_step = np.arange(0, len(Q_steps) * step_duration, 60)  # minute intervals
s_step = np.zeros_like(t_step, dtype=float)

for i, (Q_step, t_start) in enumerate(zip(Q_steps, np.arange(0, len(Q_steps) * step_duration, step_duration))):
    idx_step = (t_step >= t_start) & (t_step < t_start + step_duration)
    t_since_start = t_step[idx_step] - t_start + 60  # avoid t=0
    s_step[idx_step] = theis_drawdown(r_w, t_since_start, Q_step / 1000, T, S)
    # Add cumulative effect from previous steps
    for j in range(i):
        t_since_prev = t_step[idx_step] - j * step_duration + 60
        s_step[idx_step] += theis_drawdown(r_w, t_since_prev, Q_steps[j] / 1000, T, S)

ax4.plot(t_step / 3600, s_step, 'b-', linewidth=2.5)
for i, t_start in enumerate(np.arange(0, len(Q_steps) * step_duration, step_duration)):
    ax4.axvline(x=t_start / 3600, color='red', linestyle='--', alpha=0.7)
    if i < len(Q_steps):
        ax4.text(t_start / 3600 + 0.2, ax4.get_ylim()[1] * 0.9, f'$Q={Q_steps[i]:.0f}$ L/s',
                fontsize=9, color='red', fontweight='bold')
ax4.set_xlabel('Time (hours)', fontsize=11)
ax4.set_ylabel('Drawdown at well (m)', fontsize=11)
ax4.set_title('Step-Drawdown Test Simulation', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('groundwater_flow_well_performance.pdf', dpi=150, bbox_inches='tight')
plt.close()

# Calculate optimal pumping rate (max efficiency > 70%)
optimal_Q_idx = np.where(efficiency > 70)[0]
if len(optimal_Q_idx) > 0:
    optimal_Q = Q_range[optimal_Q_idx[-1]] * 1000  # L/s
else:
    optimal_Q = Q_range[0] * 1000
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{groundwater_flow_well_performance.pdf}
\caption{Well performance characteristics for aquifer with $T = \py{f'{T:.2e}'}$ m²/s and
$S = \py{f'{S:.2e}'}$. Top-left panel shows drawdown components: formation loss (blue) increases
linearly with pumping rate following the Theis equation, while well losses (orange shading)
contribute additional quadratic head loss from turbulent flow in the well screen and gravel pack.
Top-right panel displays specific capacity decline from \py{f'{specific_capacity[0]*1000:.1f}'} to
\py{f'{specific_capacity[-1]*1000:.1f}'} L/s/m as pumping rate increases, reflecting non-linear
well losses. Bottom-left panel shows well efficiency decreasing from \py{f'{efficiency[0]:.1f}'}\\%
to \py{f'{efficiency[-1]:.1f}'}\\% with the 70\\% threshold indicated by red dashed line,
suggesting optimal pumping rate near \py{f'{optimal_Q:.1f}'} L/s. Bottom-right panel presents
a step-drawdown test simulation with four pumping steps (10, 15, 20, 25 L/s), each lasting 3 hours,
used to determine formation loss and well loss coefficients.}
\end{figure}

\section{Results Summary}

\begin{pycode}
print(r"\begin{table}[H]")
print(r"\centering")
print(r"\caption{Aquifer Hydraulic Properties and Well Characteristics}")
print(r"\begin{tabular}{lcc}")
print(r"\toprule")
print(r"Parameter & Symbol & Value \\")
print(r"\midrule")
print(f"Hydraulic conductivity & $K$ & {K:.2e} m/s \\\\")
print(f"Aquifer thickness & $b$ & {b:.1f} m \\\\")
print(f"Transmissivity & $T$ & {T:.2e} m$^2$/s \\\\")
print(f"Storativity & $S$ & {S:.2e} (dimensionless) \\\\")
print(f"Pumping rate & $Q$ & {Q*1000:.1f} L/s \\\\")
print(f"Well radius & $r_w$ & {r_w:.2f} m \\\\")
print(r"\midrule")
print(f"Cooper-Jacob slope & $\\Delta s$/log cycle & {slope_cj:.3f} m \\\\")
print(f"Estimated transmissivity & $T_{est}$ & {T_estimated:.2e} m$^2$/s \\\\")
print(f"Estimated storativity & $S_{est}$ & {S_estimated:.2e} \\\\")
print(f"Parameter error (T) & $\\epsilon_T$ & {abs(T_estimated-T)/T*100:.1f}\\% \\\\")
print(f"Parameter error (S) & $\\epsilon_S$ & {abs(S_estimated-S)/S*100:.1f}\\% \\\\")
print(r"\midrule")
print(f"Max drawdown (24 hr) & $s_{max}$ & {max_drawdown:.2f} m \\\\")
print(f"Drawdown at 100 m & $s(r=100)$ & {drawdown_at_100m:.3f} m \\\\")
print(f"Optimal pumping rate & $Q_{opt}$ & {optimal_Q:.1f} L/s \\\\")
print(f"Specific capacity (15 L/s) & SC & {specific_capacity[5]*1000:.2f} L/s/m \\\\")
print(r"\bottomrule")
print(r"\end{tabular}")
print(r"\label{tab:results}")
print(r"\end{table}")
\end{pycode}

\section{Discussion}

\begin{example}[Field Application of Theis Analysis]
Consider a municipal water supply well pumping from a confined sandstone aquifer. Step-drawdown
testing reveals specific capacity of 8.5 L/s/m at the design pumping rate of 15 L/s. Using
the Theis equation with parameters determined from a 24-hour constant-rate test, engineers can:
\begin{itemize}
\item Predict long-term drawdown to ensure adequate saturated thickness
\item Design optimal well spacing in a wellfield to minimize interference
\item Estimate sustainable yield based on recharge and boundary conditions
\item Assess potential for saltwater intrusion in coastal aquifers
\end{itemize}
\end{example}

\begin{remark}[Assumptions and Limitations]
The Theis solution assumes: (1) homogeneous, isotropic aquifer; (2) horizontal flow;
(3) instantaneous release from storage; (4) fully penetrating well; (5) infinite areal extent.
Real aquifers often violate these assumptions, requiring modifications such as:
\begin{itemize}
\item \textbf{Hantush-Jacob} solution for leaky confined aquifers \cite{hantush1955non}
\item \textbf{Neuman} solution for unconfined aquifers with delayed yield \cite{neuman1972theory}
\item \textbf{Boulton} formulation for dual-porosity systems \cite{boulton1954unsteady}
\item \textbf{Anisotropic} formulations with directional hydraulic conductivity
\end{itemize}
\end{remark}

\begin{remark}[Numerical Methods]
The finite difference solution presented here uses iterative relaxation (SOR method) for
steady-state problems. For transient simulations, explicit or implicit time-stepping schemes
solve the time-dependent equation. Implicit methods (Crank-Nicolson, backward Euler) offer
unconditional stability but require solving linear systems at each timestep. Modern groundwater
codes like MODFLOW implement sophisticated preconditioned conjugate gradient solvers for
large 3D domains \cite{wang2000introduction}.
\end{remark}

\subsection{Comparison with Analytical Solutions}

The numerical steady-state solution shows excellent agreement with the Thiem equation beyond
a few grid cells from the well. Near the well, discretization error increases because the
finite grid cannot properly resolve the logarithmic singularity at $r = r_w$. Adaptive mesh
refinement or analytical element methods can improve near-well accuracy.

\section{Conclusions}

This analysis demonstrates comprehensive application of groundwater flow theory to practical
aquifer characterization and well hydraulics:

\begin{enumerate}
\item Pumping test analysis using the Cooper-Jacob method yields transmissivity
$T = \py{f"{T_estimated:.2e}"}$ m²/s and storativity $S = \py{f"{S_estimated:.2e}"}$,
with parameter errors of \py{f"{abs(T_estimated-T)/T*100:.1f}"}\\% and
\py{f"{abs(S_estimated-S)/S*100:.1f}"}\\% respectively

\item Two-dimensional finite difference solution reveals the radial flow field converging
toward the pumping well with maximum drawdown of \py{f"{max_drawdown:.2f}"} m at steady state

\item Sensitivity analysis confirms that transmissivity controls late-time drawdown and
Cooper-Jacob slope, while storativity primarily affects early-time response and timing

\item Well performance analysis indicates optimal pumping rate near \py{f"{optimal_Q:.1f}"} L/s
to maintain well efficiency above 70\\%, balancing production against well losses

\item Numerical and analytical solutions agree within discretization error, validating
both approaches for engineering applications
\end{enumerate}

The methods presented form the foundation for groundwater resource assessment, sustainable
yield determination, and wellfield optimization in water supply and remediation applications.

\bibliographystyle{unsrt}
\bibliography{groundwater_flow}

\end{document}