Part IV: Applied Earth Science | Chapter 3

Hydrogeology

Darcy's law, aquifer hydraulics, well equations, and contaminant transport

4.1 Darcy's Law and Hydraulic Conductivity

Henry Darcy (1856) established experimentally that the volumetric flow rate through a porous medium is proportional to the hydraulic gradient. For one-dimensional flow:

$$Q = -KA\frac{dh}{dx}$$

where $Q$ is the volumetric flow rate (m$^3$/s), $K$ is hydraulic conductivity (m/s), $A$ is cross-sectional area, and $dh/dx$ is the hydraulic gradient. The specific discharge (Darcy velocity) is:

$$q = \frac{Q}{A} = -K\frac{dh}{dx}$$

Derivation 1: Darcy's Law from Navier-Stokes

Darcy's law can be derived from the Navier-Stokes equations for slow viscous flow through a porous medium. At low Reynolds number (creeping flow), inertial terms are negligible:

$$\nabla P = \mu \nabla^2 \mathbf{v} - \rho g \hat{z}$$

Averaging over a representative elementary volume (REV) of the porous medium yields:

$$\mathbf{q} = -\frac{k}{\mu}(\nabla P + \rho g \hat{z}) = -\frac{k \rho g}{\mu}\nabla h = -K\nabla h$$

where $k$ is intrinsic permeability (m$^2$), $\mu$ is dynamic viscosity, and$h = z + P/(\rho g)$ is the hydraulic head. The hydraulic conductivity is:

$$K = \frac{k\rho g}{\mu}$$

This shows that $K$ depends on both the medium (through $k$) and the fluid (through $\rho$ and $\mu$). The actual pore velocity is faster than the Darcy velocity by the factor $1/n$ where $n$ is porosity:$v_{\text{pore}} = q/n$.

4.2 Aquifer Types and Groundwater Flow

Derivation 2: The Groundwater Flow Equation

Combining Darcy's law with the conservation of mass gives the groundwater flow equation. For a confined aquifer of thickness $b$, transmissivity $T = Kb$, and storativity $S$:

$$T\nabla^2 h = S\frac{\partial h}{\partial t} - W$$

where $W$ is a source/sink term (recharge or pumping). For an unconfined aquifer, the Dupuit approximation (horizontal flow, hydrostatic pressure) gives:

$$\frac{\partial}{\partial x}\left(Kh\frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y}\left(Kh\frac{\partial h}{\partial y}\right) = S_y\frac{\partial h}{\partial t} - W$$

where $S_y$ is the specific yield (drainable porosity, typically 0.1-0.3). This is nonlinear because the saturated thickness $h$ appears in the coefficients.

Aquifer types include: Confined (bounded above and below by aquitards, water under pressure), Unconfined (water table at atmospheric pressure),Perched (local saturated zone above main water table), andArtesian (confined with potentiometric surface above ground).

4.3 Well Hydraulics: The Theis Equation

Derivation 3: Transient Radial Flow to a Well

C.V. Theis (1935) solved the problem of transient flow to a pumping well in a confined aquifer, by analogy with heat conduction. The governing equation in radial coordinates is:

$$T\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial s}{\partial r}\right) = S\frac{\partial s}{\partial t}$$

where $s = h_0 - h$ is the drawdown. With boundary conditions of constant pumping rate $Q$ and initially uniform head, the solution is:

$$s(r, t) = \frac{Q}{4\pi T}W(u), \quad u = \frac{r^2 S}{4Tt}$$

where $W(u)$ is the well function (exponential integral):

$$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$$

For small $u$ (large time or small distance), the Cooper-Jacob approximation applies:

$$s \approx \frac{Q}{4\pi T}\left(-0.5772 - \ln u\right) = \frac{2.3Q}{4\pi T}\log_{10}\frac{2.25Tt}{r^2 S}$$

This logarithmic approximation is the basis of the Cooper-Jacob straight-line method for pumping test analysis. A plot of $s$ vs $\log t$ gives a straight line with slope $\Delta s = 2.3Q/(4\pi T)$ per log cycle.

4.4 Contaminant Transport

Derivation 4: The Advection-Dispersion Equation

Contaminant transport in groundwater involves advection (transport with flowing water), dispersion (mechanical mixing and molecular diffusion), and reactions. The 1D advection-dispersion equation (ADE) is:

$$\frac{\partial C}{\partial t} = D_L \frac{\partial^2 C}{\partial x^2} - v\frac{\partial C}{\partial x} - \lambda C$$

where $C$ is concentration, $v = q/n$ is pore velocity, $D_L = \alpha_L v + D_m$is the longitudinal dispersion coefficient ($\alpha_L$ is dispersivity, $D_m$ is molecular diffusion), and $\lambda$ is a first-order decay constant.

For a continuous source at $x = 0$ with concentration $C_0$ starting at $t = 0$(Ogata-Banks solution, no decay):

$$C(x,t) = \frac{C_0}{2}\left[\text{erfc}\left(\frac{x - vt}{2\sqrt{D_L t}}\right) + \exp\left(\frac{vx}{D_L}\right)\text{erfc}\left(\frac{x + vt}{2\sqrt{D_L t}}\right)\right]$$

Derivation 5: Retardation by Sorption

Many contaminants are retarded by sorption onto aquifer solids. For linear, equilibrium sorption ($S = K_d C$ where $K_d$ is the distribution coefficient):

$$R\frac{\partial C}{\partial t} = D_L \frac{\partial^2 C}{\partial x^2} - v\frac{\partial C}{\partial x}$$

where the retardation factor is:

$$R = 1 + \frac{\rho_b K_d}{n}$$

with $\rho_b$ being the bulk density and $n$ the porosity. The contaminant moves at velocity $v_c = v/R$, slower than the groundwater. For organic contaminants, $K_d$ is often estimated from the organic carbon partition coefficient: $K_d = K_{oc} f_{oc}$.

4.5 Applications in Water Resources

Hydrogeology is critical for water supply, environmental remediation, and geotechnical engineering. Over 2 billion people depend on groundwater as their primary water source. Key challenges include sustainable extraction rates, saltwater intrusion in coastal aquifers, contamination from industrial and agricultural sources, and managing transboundary aquifer resources.

The Ogallala Aquifer (High Plains, USA) illustrates the sustainability challenge: withdrawals for irrigation far exceed natural recharge, and water levels have declined by more than 30 m in parts of Kansas and Texas since 1950. Similar overdraft situations exist in the North China Plain, the Arabian Peninsula, and the Indo-Gangetic Basin.

Managed Aquifer Recharge

Managed aquifer recharge (MAR) deliberately replenishes groundwater by spreading surface water in infiltration basins, injecting treated water through wells, or enhancing natural recharge. MAR is practiced in over 100 countries and is critical for water banking (storing wet-season surplus for dry-season use), seawater intrusion barriers, and aquifer storage and recovery (ASR) systems. The Orange County (California) Groundwater Replenishment System is the world's largest, injecting 400,000 m$^3$/day of purified recycled water.

Groundwater-Surface Water Interaction

The interaction between groundwater and surface water (rivers, lakes, wetlands) is critical for ecosystem health and water management. Gaining streams receive groundwater (baseflow), while losing streams recharge the aquifer. The fraction of streamflow from groundwater (the baseflow index) ranges from ~0.1 in arid regions to ~0.9 in humid regions with permeable geology. Temperature and chemical tracers (e.g., radon-222, electrical conductivity) are used to quantify these fluxes.

4.6 Groundwater Recharge and the Water Balance

The water balance equation for a catchment relates precipitation $P$, evapotranspiration$ET$, runoff $R_{\text{off}}$, and change in storage $\Delta S$:

$$P = ET + R_{\text{off}} + \Delta S$$

Groundwater recharge is the portion of precipitation that percolates past the root zone to reach the water table. In humid temperate regions, recharge is typically 15-30% of precipitation; in arid regions, it can be less than 1%. Recharge estimation methods include chloride mass balance, water table fluctuation, and lysimeter measurements.

Saltwater Intrusion

In coastal aquifers, freshwater floats on denser saltwater. The Ghyben-Herzberg relation gives the depth to the freshwater-saltwater interface:

$$z_s = \frac{\rho_f}{\rho_s - \rho_f} h_f \approx 40 h_f$$

where $h_f$ is the freshwater head above sea level, $\rho_f = 1000$ kg/m$^3$, and $\rho_s = 1025$ kg/m$^3$. This means that for every meter of freshwater head above sea level, the interface extends ~40 m below sea level. Over-pumping that lowers the water table by 1 m can cause the interface to rise by 40 m, contaminating wells with saltwater.

Capture Zone Analysis

The capture zone of a pumping well is the region from which groundwater will eventually reach the well. For uniform flow with velocity $U$ and a well pumping at rate $Q$, the capture zone width far upgradient approaches:

$$y_{\text{max}} = \frac{Q}{2bU}$$

where $b$ is the aquifer thickness. The stagnation point (where groundwater velocity is zero) is located at distance $x_0 = -Q/(2\pi bU)$ downgradient of the well. Capture zone delineation is essential for wellhead protection and remediation system design.

Computational Lab: Hydrogeology

Darcy Flow, Theis Well Function, and Contaminant Transport

Python

script.py361 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# HYDROGEOLOGY SIMULATION
# ============================================================
print("HYDROGEOLOGY SIMULATION")
print("=" * 70)

# ============================================================
# 1. Darcy's Law: Hydraulic Conductivity of Materials
# ============================================================
print()
print("1. Hydraulic Properties of Geological Materials")
print("-" * 70)

materials = [
    ("Gravel", 1e-1, 0.30, 0.25),
    ("Coarse sand", 1e-2, 0.35, 0.27),
    ("Fine sand", 1e-4, 0.40, 0.21),
    ("Silt", 1e-6, 0.45, 0.08),
    ("Clay", 1e-9, 0.50, 0.03),
    ("Sandstone", 1e-5, 0.15, 0.10),
    ("Limestone", 1e-5, 0.10, 0.05),
    ("Fractured granite", 1e-6, 0.02, 0.01),
    ("Unfractured granite", 1e-12, 0.005, 0.001),
]

print(f"  {'Material':22} {'K (m/s)':14} {'Porosity':12} {'Spec. yield':14} {'Class'}")
print("  " + "-" * 65)

for name, K, n, Sy in materials:
    cls = "Aquifer" if K > 1e-7 else ("Aquitard" if K > 1e-10 else "Aquiclude")
    print(f"  {name:20} {K:12.1e} {n:10.2f} {Sy:12.2f} {cls}")

# Darcy flow example
print()
print("  Darcy flow example:")
K_ex = 1e-4  # m/s
i_ex = 0.01  # hydraulic gradient
A_ex = 100   # m^2 cross-section
n_ex = 0.30

q = K_ex * i_ex
Q = q * A_ex
v_pore = q / n_ex
print(f"    K = {K_ex:.1e} m/s, gradient = {i_ex}, area = {A_ex} m^2")
print(f"    Darcy velocity q = {q:.2e} m/s = {q*86400:.2f} m/day")
print(f"    Flow rate Q = {Q:.4f} m^3/s = {Q*86400:.1f} m^3/day")
print(f"    Pore velocity = {v_pore:.2e} m/s = {v_pore*86400:.2f} m/day")

# ============================================================
# 2. Theis Well Function
# ============================================================
print()
print()
print("2. Theis Well Function W(u)")
print("-" * 55)

def well_function(u):
    # Series approximation for W(u) = -Ei(-u)
    if u < 1:
        W = -0.5772 - np.log(u) + u - u**2/(2*2) + u**3/(3*6) - u**4/(4*24) + u**5/(5*120)
    else:
        # Asymptotic expansion for large u
        W = np.exp(-u)/u * (1 - 1/u + 2/u**2 - 6/u**3)
    return max(W, 0)

print(f"  {'u':14} {'W(u)':14} {'Cooper-Jacob':14} {'Difference (%)'}")
print("  " + "-" * 50)

for u in [1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 0.05, 0.1, 0.5, 1.0, 2.0, 5.0]:
    Wu = well_function(u)
    CJ = -0.5772 - np.log(u)  # Cooper-Jacob approx
    diff = abs(Wu - CJ) / Wu * 100 if Wu > 0 else 0
    print(f"  {u:12.1e} {Wu:12.4f} {CJ:12.4f} {diff:10.1f}")

# ============================================================
# 3. Pumping Test Analysis
# ============================================================
print()
print()
print("3. Pumping Test: Drawdown vs Time and Distance")
print("-" * 65)

T = 500  # m^2/day transmissivity
S = 1e-4  # storativity (confined)
Q_pump = 1000  # m^3/day

print(f"  T = {T} m^2/day, S = {S:.0e}, Q = {Q_pump} m^3/day")
print()

# Drawdown at r = 50 m over time
print(f"  Drawdown at r = 50 m:")
print(f"  {'Time (min)':14} {'Time (day)':14} {'u':14} {'W(u)':12} {'s (m)'}")
print("  " + "-" * 55)

r = 50.0
for t_min in [1, 5, 10, 30, 60, 120, 360, 720, 1440, 2880, 7200]:
    t_day = t_min / 1440.0
    u = r**2 * S / (4 * T * t_day)
    Wu = well_function(u)
    s = Q_pump / (4 * np.pi * T) * Wu
    print(f"  {t_min:12} {t_day:12.4f} {u:12.2e} {Wu:10.4f} {s:8.3f}")

# Drawdown at t = 1 day vs distance
print()
print(f"  Drawdown at t = 1 day:")
print(f"  {'r (m)':12} {'u':14} {'W(u)':12} {'s (m)'}")
print("  " + "-" * 45)

for r_val in [1, 5, 10, 25, 50, 100, 200, 500, 1000, 2000]:
    u = r_val**2 * S / (4 * T * 1.0)
    Wu = well_function(u)
    s = Q_pump / (4 * np.pi * T) * Wu
    print(f"  {r_val:10} {u:12.2e} {Wu:10.4f} {s:8.3f}")

# ============================================================
# 4. Contaminant Transport (1D ADE)
# ============================================================
print()
print()
print("4. 1D Contaminant Transport (Advection-Dispersion)")
print("-" * 65)

v_flow = 0.5  # m/day pore velocity
alpha_L = 10.0  # m dispersivity
D_L = alpha_L * v_flow  # m^2/day
C0 = 100.0  # initial concentration (mg/L)

print(f"  Pore velocity = {v_flow} m/day, dispersivity = {alpha_L} m")
print(f"  Dispersion coefficient = {D_L} m^2/day")
print()

def erfc(x):
    # Approximation of complementary error function
    t = 1.0 / (1.0 + 0.3275911 * abs(x))
    poly = t * (0.254829592 + t * (-0.284496736 + t * (1.421413741 +
           t * (-1.453152027 + t * 1.061405429))))
    result = poly * np.exp(-x**2)
    return result if x >= 0 else 2.0 - result

# Concentration at various x and t
print(f"  C/C0 at different distances and times (no retardation):")
print(f"  {'x (m)':10}", end="")
for t_day in [100, 200, 365, 500, 730]:
    print(f"{'t=' + str(t_day) + 'd':12}", end="")
print()
print("  " + "-" * 60)

for x in [0, 25, 50, 100, 150, 200, 250, 300, 400, 500]:
    print(f"  {x:8}", end="")
    for t_day in [100, 200, 365, 500, 730]:
        if t_day > 0:
            arg1 = (x - v_flow * t_day) / (2 * np.sqrt(D_L * t_day))
            # Simplified (neglect second term for Peclet >> 1)
            C_ratio = 0.5 * erfc(arg1)
            print(f"  {C_ratio:10.3f}", end="")
        else:
            print(f"  {'---':10}", end="")
    print()

# ============================================================
# 5. Retardation Factors
# ============================================================
print()
print()
print("5. Retardation Factors for Common Contaminants")
print("-" * 65)

rho_b = 1800  # kg/m^3 bulk density
n_por = 0.30
f_oc = 0.005  # fraction organic carbon

contaminants = [
    ("Chloride", 0, 0),
    ("Nitrate", 0, 0),
    ("Benzene", 83, 0),
    ("Toluene", 190, 0),
    ("TCE", 126, 0),
    ("PCE", 364, 0),
    ("Naphthalene", 1300, 0),
    ("Cadmium", 0, 5.0),
    ("Lead", 0, 500),
    ("Chromium(VI)", 0, 0.1),
]

print(f"  {'Contaminant':16} {'Koc (L/kg)':14} {'Kd (L/kg)':12} {'R':8} {'v_c (m/day)'}")
print("  " + "-" * 55)

for name, Koc, Kd_direct in contaminants:
    if Koc > 0:
        Kd = Koc * f_oc / 1000  # L/kg -> m^3/kg (divide by 1000)
    elif Kd_direct > 0:
        Kd = Kd_direct / 1000
    else:
        Kd = 0
    R = 1 + rho_b * Kd / n_por
    v_c = v_flow / R
    Kd_Lkg = Kd * 1000
    Koc_str = f"{Koc}" if Koc > 0 else "---"
    print(f"  {name:14} {Koc_str:12} {Kd_Lkg:10.3f} {R:6.1f} {v_c:8.4f}")

# ============================================================
# PLOTS
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 11))

# Plot 1: Theis drawdown cone
ax = axes[0, 0]
r_arr = np.linspace(1, 1000, 500)
times = [0.1, 0.5, 1.0, 5.0, 10.0, 30.0]

for t_day in times:
    s_arr = []
    for r_val in r_arr:
        u = r_val**2 * S / (4 * T * t_day)
        Wu = well_function(u)
        s_val = Q_pump / (4 * np.pi * T) * Wu
        s_arr.append(s_val)
    ax.plot(r_arr, s_arr, linewidth=2, label=f't = {t_day} day')

ax.set_xlabel('Distance from well (m)')
ax.set_ylabel('Drawdown (m)')
ax.set_title('Theis Drawdown Cone of Depression')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_xlim(1, 1000)
ax.set_xscale('log')

# Plot 2: Cooper-Jacob time-drawdown
ax = axes[0, 1]
t_arr = np.logspace(-3, 2, 500)  # days
r_obs = [10, 50, 100, 200]

for r_val in r_obs:
    s_arr = []
    for t_val in t_arr:
        u = r_val**2 * S / (4 * T * t_val)
        Wu = well_function(u)
        s_val = Q_pump / (4 * np.pi * T) * Wu
        s_arr.append(s_val)
    ax.semilogx(t_arr, s_arr, linewidth=2, label=f'r = {r_val} m')

# Cooper-Jacob straight line (for r=50)
t_CJ = np.logspace(-1, 2, 100)
s_CJ = 2.3 * Q_pump / (4 * np.pi * T) * np.log10(2.25 * T * t_CJ / (50**2 * S))
s_CJ = np.maximum(s_CJ, 0)
ax.semilogx(t_CJ, s_CJ, 'k--', linewidth=1, label='Cooper-Jacob (r=50m)')

ax.set_xlabel('Time (days)')
ax.set_ylabel('Drawdown (m)')
ax.set_title('Time-Drawdown Curves')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3, which='both')

# Plot 3: Contaminant breakthrough curves
ax = axes[1, 0]
x_obs = 200  # m
t_plot = np.linspace(1, 1500, 500)

for R_val, label in [(1.0, 'Conservative (R=1)'), (2.0, 'R=2'),
                      (5.0, 'R=5'), (10.0, 'R=10')]:
    v_eff = v_flow / R_val
    D_eff = D_L / R_val
    C_plot = []
    for t_val in t_plot:
        if t_val > 0:
            arg = (x_obs - v_eff * t_val) / (2 * np.sqrt(D_eff * t_val))
            C_val = 0.5 * erfc(arg) * C0
        else:
            C_val = 0
        C_plot.append(C_val)
    ax.plot(t_plot, C_plot, linewidth=2, label=label)

ax.axhline(y=50, color='red', linestyle=':', alpha=0.5, label='MCL (50 mg/L)')
ax.set_xlabel('Time (days)')
ax.set_ylabel('Concentration (mg/L)')
ax.set_title(f'Breakthrough Curves at x = {x_obs} m')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# Plot 4: Contaminant plume snapshot
ax = axes[1, 1]
x_plume = np.linspace(0, 500, 500)
t_snapshot = 365  # 1 year

for R_val, label, color in [(1.0, 'R=1', 'blue'), (3.0, 'R=3', 'red'),
                              (10.0, 'R=10', 'green')]:
    v_eff = v_flow / R_val
    D_eff = D_L / R_val
    C_plume = []
    for x_val in x_plume:
        arg = (x_val - v_eff * t_snapshot) / (2 * np.sqrt(D_eff * t_snapshot))
        C_val = 0.5 * erfc(arg) * C0
        C_plume.append(C_val)
    ax.plot(x_plume, C_plume, color=color, linewidth=2, label=label)
    # Mark front position
    front = v_eff * t_snapshot
    ax.axvline(x=front, color=color, linestyle=':', alpha=0.3)

ax.set_xlabel('Distance from source (m)')
ax.set_ylabel('Concentration (mg/L)')
ax.set_title(f'Contaminant Plume at t = {t_snapshot} days')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
plt.close()
print()
print()
print("6. Ghyben-Herzberg Saltwater Interface")
print("-" * 55)

rho_fresh = 1000.0
rho_salt = 1025.0
ratio_GH = rho_fresh / (rho_salt - rho_fresh)

print(f"  Ghyben-Herzberg ratio: {ratio_GH:.0f}")
print(f"  (1 m head above sea level => {ratio_GH:.0f} m interface below)")
print()
print(f"  {'Head above SL (m)':22} {'Interface depth (m)':22} {'Safe well depth (m)'}")
print("  " + "-" * 60)

for h_above in [0.5, 1.0, 2.0, 3.0, 5.0, 7.0, 10.0, 15.0, 20.0]:
    z_interface = ratio_GH * h_above
    safe_depth = z_interface * 0.5  # rule of thumb: pump from upper half
    print(f"  {h_above:20.1f} {z_interface:20.0f} {safe_depth:10.0f}")

# Capture zone width
print()
print()
print("7. Well Capture Zone Dimensions")
print("-" * 55)

b_aq = 20.0  # m aquifer thickness
K_cap = 1e-4  # m/s
i_cap = 0.005  # hydraulic gradient

U_cap = K_cap * i_cap  # Darcy velocity

print(f"  K = {K_cap:.1e} m/s, gradient = {i_cap}, b = {b_aq} m")
print(f"  Ambient flow: U = {U_cap:.2e} m/s = {U_cap*86400:.2f} m/day")
print()
print(f"  {'Q_pump (m3/day)':18} {'y_max (m)':14} {'x_stag (m)':14} {'Capture area (m2)'}")
print("  " + "-" * 55)

for Q_cap in [10, 50, 100, 200, 500, 1000, 2000]:
    Q_ms = Q_cap / 86400  # m^3/s
    y_max = Q_ms / (2 * b_aq * U_cap)
    x_stag = Q_ms / (2 * np.pi * b_aq * U_cap)
    area = np.pi * y_max * x_stag * 3  # rough capture zone area
    print(f"  {Q_cap:16} {y_max:12.0f} {x_stag:12.0f} {area:12.0f}")

print()
print()
print("[Chart saved: Theis drawdown, time-drawdown, breakthrough, plume snapshot]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Saltwater Intrusion and Aquifer Sustainability Analysis

Python

script.py73 lines

import numpy as np

print("SALTWATER INTRUSION AND AQUIFER SUSTAINABILITY")
print("=" * 70)

# Ghyben-Herzberg analysis
rho_f = 1000.0
rho_s = 1025.0
ratio = rho_f / (rho_s - rho_f)

print()
print("1. Ghyben-Herzberg Interface Position")
print("-" * 55)
print(f"  Density ratio factor: {ratio:.0f}")
print()
print(f"  {'Pumping scenario':25} {'Head (m)':12} {'Interface (m)':14} {'Safe depth (m)'}")
print("  " + "-" * 60)

scenarios = [
    ("Natural (no pumping)", 3.0),
    ("Light pumping", 2.0),
    ("Moderate pumping", 1.0),
    ("Heavy pumping", 0.5),
    ("Over-pumping", 0.2),
    ("Critical", 0.1),
]

for name, h in scenarios:
    z_int = ratio * h
    safe = z_int * 0.4
    print(f"  {name:23} {h:10.1f} {z_int:12.0f} {safe:10.0f}")

# Aquifer sustainability
print()
print()
print("2. Aquifer Sustainability: Recharge vs Withdrawal")
print("-" * 60)

area = 1000e6  # 1000 km^2 in m^2
precip = 0.6   # m/yr
recharge_frac = 0.15

recharge_rate = precip * recharge_frac * area  # m^3/yr
recharge_Mm3 = recharge_rate / 1e6

print(f"  Area = 1000 km^2, Precipitation = {precip*1000:.0f} mm/yr")
print(f"  Recharge fraction = {recharge_frac*100:.0f}%")
print(f"  Sustainable yield = {recharge_Mm3:.0f} Mm^3/yr = {recharge_rate/86400/365:.0f} m^3/s")
print()
print(f"  {'Withdrawal (Mm3/yr)':24} {'Ratio W/R':14} {'Sustainability':18} {'Depletion time'}")
print("  " + "-" * 65)

storage = 5000  # Mm^3 total storage

for W in [10, 30, 60, 90, 120, 150, 200, 300]:
    ratio_wr = W / recharge_Mm3
    if ratio_wr < 0.5:
        sust = "Sustainable"
    elif ratio_wr < 1.0:
        sust = "Marginal"
    else:
        sust = "UNSUSTAINABLE"
    net_loss = max(0, W - recharge_Mm3)
    if net_loss > 0:
        depl_yr = storage / net_loss
        depl = f"{depl_yr:.0f} yr"
    else:
        depl = "indefinite"
    print(f"  {W:22} {ratio_wr:12.2f} {sust:16} {depl}")

print()
print("Complete: sustainability analysis finished.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Isotope Geology Next: Natural Hazards →