Hydrology: units, equations and validation

Domain hub · History · Examples · API reference

Conventions and model limits

Lengths and depths are metres, times seconds and rates m/s; divide a rain rate in mm/h by 3.6e6 and convert days with units.convert(value, "day", "s"). The water balance accepts any common depth rate, such as mm/day, used consistently. Evaporation rates are kg/(m^2 s), which equals mm/s of liquid water; multiply by 86400 for mm/day. Vapour pressures are pascals, and saturation_slope is the derivative of earthsciencekit.atmosphere.saturation_vapor_pressure() with the same constant latent heat LATENT_HEAT_VAPORIZATION = 2.5e6 J/kg; the von Karman constant VON_KARMAN = 0.41 sets the neutral aerodynamic resistance. Both live in earthsciencekit.constants.

Darcy fluxes are positive down the head gradient, and hydraulic_head adds the elevation to the gauge pressure head. Aquifers are homogeneous and isotropic: Dupuit’s water table has no recharge and a horizontal base, and the Theis and Cooper-Jacob drawdowns assume a fully penetrating well of negligible radius in an infinite confined aquifer. Cooper-Jacob is the late-time limit, within 1% of Theis for u < 0.01.

In the unsaturated zone heights and fluxes are positive upward and pressure heads negative under suction. richards_steady_head solves the steady Richards equation for Gardner’s exponential conductivity above a water table at height zero; it rejects downward fluxes at or beyond the saturated conductivity and heights above the limit of capillary supply. Green-Ampt and Horton rates are infiltration capacities of a uniform soil; Green-Ampt times start when the surface ponds, and the suction head is a positive magnitude.

Channel formulas are for steady uniform flow in SI units: Chézy’s coefficient is in m^0.5/s and Manning’s n in s/m^(1/3). The unit hydrograph is linear and time-invariant: unit_hydrograph_convolution takes excess rainfall depths per interval and ordinates in m^3/s per metre of excess at the same interval, and the Nash shape integrates to one per second. muskingum_coefficients returns [C0, C1, C2] on the last axis; muskingum_routing starts from a steady state and needs scalar parameters. Stream orders follow Strahler (1952). TOPMODEL indices use the upslope area per unit contour width in metres, and nonpositive local deficits mark saturated ground.

Input contracts

Algebraic models accept finite real scalars or nonempty arrays. Booleans, complex inputs, nonnumeric values and NaN/Inf are rejected. Physical bounds apply element by element; NumPy broadcasting applies between arguments, except for the vector inputs of unit_hydrograph_convolution and muskingum_routing. No unit metadata are carried in raw numbers. Bounds protect mathematical domains and obvious physical contradictions, but do not establish model applicability.

Equation and output map

water_balance

earthsciencekit.hydrology.water_balance()

Output: m/s. Perrault (1674): a closed catchment with no groundwater exchange across its divide; any common depth rate (such as mm/day) may be used throughout.

precipitation - evapotranspiration - runoff

runoff_coefficient

earthsciencekit.hydrology.runoff_coefficient()

Output: 1. Both depths over the same area and period; long periods make the storage change negligible.

runoff / precipitation

hydraulic_radius

earthsciencekit.hydrology.hydraulic_radius()

Output: m. Rectangular section; tends to the depth for a wide channel.

width * depth / (width + 2 * depth)

chezy_velocity

earthsciencekit.hydrology.chezy_velocity()

Output: m/s. Chezy (1775): steady uniform turbulent flow; the coefficient is not a constant but depends on roughness and depth.

chezy_coefficient * np.sqrt(hydraulic_radius * slope)

dalton_evaporation

earthsciencekit.hydrology.dalton_evaporation()

Output: kg/(m^2 s). Dalton (1802), a bulk mass-transfer law; the surface temperature it needs is rarely measured.

transfer_coefficient * (saturation_vapor_pressure - vapor_pressure)

hydraulic_head

earthsciencekit.hydrology.hydraulic_head()

Output: m. Slow flow with negligible velocity head and uniform fluid density.

elevation + pressure / (density * gravity)

darcy_flux

earthsciencekit.hydrology.darcy_flux()

Output: m/s. Darcy (1856): laminar flow with a pore Reynolds number below about 1-10 in a saturated isotropic medium.

-hydraulic_conductivity * head_gradient

seepage_velocity

earthsciencekit.hydrology.seepage_velocity()

Output: m/s. Average linear velocity; dispersion spreads tracers about it.

darcy_flux / porosity

dupuit_head

earthsciencekit.hydrology.dupuit_head()

Output: m. Dupuit (1863): no recharge, a horizontal base, and slopes small enough that the head is uniform on each vertical.

np.sqrt(upstream_head**2 - (upstream_head**2 - downstream_head**2) * distance / length)

dupuit_discharge

earthsciencekit.hydrology.dupuit_discharge()

Output: m^2/s. Dupuit (1863), the Dupuit-Forchheimer discharge formula, exact for the discharge even where the head profile is not.

hydraulic_conductivity * (upstream_head**2 - downstream_head**2) / (2 * length)

manning_velocity

earthsciencekit.hydrology.manning_velocity()

Output: m/s. Manning (1891), SI form; fully rough turbulent flow, steady and uniform.

hydraulic_radius**(2 / 3) * np.sqrt(slope) / roughness

chezy_from_manning

earthsciencekit.hydrology.chezy_from_manning()

Output: m^0.5/s. Manning (1891), SI units.

hydraulic_radius**(1 / 6) / roughness

green_ampt_rate

earthsciencekit.hydrology.green_ampt_rate()

Output: m/s. Green and Ampt (1911): a piston-like front in a uniform soil with a negligible ponding depth; the rate tends to the saturated conductivity.

saturated_conductivity * (1 + suction_head * moisture_deficit / cumulative_infiltration)

green_ampt_time

earthsciencekit.hydrology.green_ampt_time()

Output: s. Green and Ampt (1911), ponded from time zero; early on the depth grows as the square root of time.

(cumulative_infiltration - suction_head * moisture_deficit * np.log(1 + cumulative_infiltration / (suction_head * moisture_deficit))) / saturated_conductivity

gardner_conductivity

earthsciencekit.hydrology.gardner_conductivity()

Output: m/s. Gardner (1958), a one-parameter fit without hysteresis; it makes the steady Richards equation solvable in closed form.

saturated_conductivity * np.exp(sorptive_number * pressure_head)

unsaturated_darcy_flux

earthsciencekit.hydrology.unsaturated_darcy_flux()

Output: m/s. Richards (1931): the total head is pressure_head + z and air stays at atmospheric pressure.

-conductivity * (pressure_head_gradient + 1)

richards_steady_head

earthsciencekit.hydrology.richards_steady_head()

Output: m. Richards (1931) steady equation flux = -K(psi) (dpsi/dz + 1) with Gardner’s (1958) conductivity; an upward flux can only be sustained below np.log(1 + saturated_conductivity / flux) / sorptive_number.

np.log(((flux + saturated_conductivity) * np.exp(-sorptive_number * height) - flux) / saturated_conductivity) / sorptive_number

nash_unit_hydrograph

earthsciencekit.hydrology.nash_unit_hydrograph()

Output: 1/s. Nash (1957), a gamma density with mean lag reservoirs * storage_constant; a smooth synthetic shape for Sherman’s unit hydrograph, linear and time-invariant.

time**(reservoirs - 1) * np.exp(-time / storage_constant) / (storage_constant**reservoirs * scipy.special.gamma(reservoirs))

unit_hydrograph_convolution

earthsciencekit.hydrology.unit_hydrograph_convolution()

Output: m^3/s. Sherman (1932): linear superposition and time invariance; baseflow and losses are removed beforehand.

np.convolve(excess_rainfall, unit_hydrograph)

horton_infiltration_rate

earthsciencekit.hydrology.horton_infiltration_rate()

Output: m/s. Horton (1933, 1940), an empirical law for rain at or above capacity throughout.

final_rate + (initial_rate - final_rate) * np.exp(-decay_constant * time)

horton_cumulative_infiltration

earthsciencekit.hydrology.horton_cumulative_infiltration()

Output: m. Horton (1933, 1940), with rain at or above capacity from time zero.

final_rate * time + (initial_rate - final_rate) * (1 - np.exp(-decay_constant * time)) / decay_constant

well_function

earthsciencekit.hydrology.well_function()

Output: 1. The exponential integral E1, -0.5772 - ln(u) + u - ... for small u.

scipy.special.exp1(u)

theis_drawdown

earthsciencekit.hydrology.theis_drawdown()

Output: m. Theis (1935): constant rate, a fully penetrating well of negligible radius in an infinite, homogeneous, isotropic confined aquifer with water released instantly from storage.

pumping_rate / (4 * np.pi * transmissivity) * scipy.special.exp1(radius**2 * storativity / (4 * transmissivity * time))

muskingum_coefficients

earthsciencekit.hydrology.muskingum_coefficients()

Output: 1. McCarthy (1938), storage K (X I + (1 - X) O); all coefficients are nonnegative when 2 K X <= time_step <= 2 K (1 - X).

[time_step - 2 K X, time_step + 2 K X, 2 K (1 - X) - time_step] / (2 K (1 - X) + time_step)

muskingum_routing

earthsciencekit.hydrology.muskingum_routing()

Output: m^3/s. McCarthy (1938): O[j+1] = C0 I[j+1] + C1 I[j] + C2 O[j]; choose the time step so that no coefficient is negative.

O[j + 1] = C0 * I[j + 1] + C1 * I[j] + C2 * O[j]

stream_number

earthsciencekit.hydrology.stream_number()

Output: 1. Horton’s law of stream numbers (1945), an empirical average; Strahler’s (1952) ordering is the usual convention.

bifurcation_ratio**(highest_order - order)

stream_length

earthsciencekit.hydrology.stream_length()

Output: m. Horton’s law of stream lengths (1945), an empirical average over a network.

first_order_length * length_ratio**(order - 1)

cooper_jacob_drawdown

earthsciencekit.hydrology.cooper_jacob_drawdown()

Output: m. Cooper and Jacob (1946): the first two terms of the Theis series, within 1% for u < 0.01 and meaningless (it goes negative) for u above about 0.56.

pumping_rate / (4 * np.pi * transmissivity) * (-np.euler_gamma - np.log(radius**2 * storativity / (4 * transmissivity * time)))

cooper_jacob_transmissivity

earthsciencekit.hydrology.cooper_jacob_transmissivity()

Output: m^2/s. Cooper and Jacob (1946), from late-time data where u < 0.01; the storativity follows from the zero-drawdown intercept t0 as 2.25 * transmissivity * t0 / radius**2.

np.log(10) * pumping_rate / (4 * np.pi * drawdown_per_log_cycle)

saturation_slope

earthsciencekit.hydrology.saturation_slope()

Output: Pa/K. Derivative of earthsciencekit.atmosphere.saturation_vapor_pressure() (constant latent heat 2.5e6 J/kg).

saturation_vapor_pressure(temperature) * LATENT_HEAT_VAPORIZATION / (R_VAPOR * temperature**2)

psychrometric_constant

earthsciencekit.hydrology.psychrometric_constant()

Output: Pa/K. Dry-air heat capacity and constant latent heat; equal transfer coefficients for heat and vapour.

CP_DRY_AIR * pressure / (EPSILON_WATER * LATENT_HEAT_VAPORIZATION)

penman_evaporation

earthsciencekit.hydrology.penman_evaporation()

Output: kg/(m^2 s). Penman (1948): the saturation curve is linearized between air and surface temperature, and heat and vapour share one transfer coefficient; a saturated surface (open water or wet grass).

(saturation_slope(temperature) * net_radiation / LATENT_HEAT_VAPORIZATION + psychrometric_constant(pressure) * aerodynamic_evaporation) / (saturation_slope(temperature) + psychrometric_constant(pressure))

aerodynamic_resistance

earthsciencekit.hydrology.aerodynamic_resistance()

Output: s/m. Neutral stability; the measurement height above the displacement height must exceed both roughness lengths. For FAO grass (0.12 m) at 2 m it is about 208 / wind_speed.

np.log((measurement_height - displacement_height) / momentum_roughness) * np.log((measurement_height - displacement_height) / vapor_roughness) / (VON_KARMAN**2 * wind_speed)

bulk_surface_resistance

earthsciencekit.hydrology.bulk_surface_resistance()

Output: s/m. Monteith (1965) big-leaf canopy with parallel leaves; FAO grass has about 70 s/m.

stomatal_resistance / active_leaf_area_index

penman_monteith_evaporation

earthsciencekit.hydrology.penman_monteith_evaporation()

Output: kg/(m^2 s). Monteith (1965): a big-leaf canopy at one temperature, neutral stability; zero surface resistance recovers Penman with aerodynamic_evaporation = air_density * CP_DRY_AIR * vapor_pressure_deficit / (psychrometric_constant * LATENT_HEAT_VAPORIZATION * aerodynamic_resistance).

(Delta * available_energy + air_density * CP_DRY_AIR * vapor_pressure_deficit / aerodynamic_resistance) / (LATENT_HEAT_VAPORIZATION * (Delta + gamma * (1 + surface_resistance / aerodynamic_resistance)))

topographic_index

earthsciencekit.hydrology.topographic_index()

Output: 1. Beven and Kirkby (1979): steady recharge, a water table parallel to the surface and transmissivity falling exponentially with deficit.

np.log(specific_area / np.tan(slope))

topmodel_deficit

earthsciencekit.hydrology.topmodel_deficit()

Output: m. Beven and Kirkby (1979) with uniform transmissivity at saturation and quasi-steady drainage.

mean_deficit + scaling_depth * (mean_index - topographic_index)

topmodel_baseflow

earthsciencekit.hydrology.topmodel_baseflow()

Output: m/s. Beven and Kirkby (1979): exponential transmissivity profile, water table parallel to the surface.

saturated_transmissivity * np.exp(-mean_index) * np.exp(-mean_deficit / scaling_depth)

Infiltration under steady rain

earthsciencekit.hydrology.simulate_infiltration() integrates the cumulative infiltration F of a Green-Ampt soil under rain of rate r, dF/dt = min(r, K (1 + psi dtheta / F)), from F = 0 at the first sample. The surface ponds when F = K psi dtheta / (r - K), at t_p = K psi dtheta / (r (r - K)); rain no faster than K never ponds. It returns the cumulative infiltration, the infiltration rate and the infiltration-excess runoff rate r - f. Surface storage, redistribution and a water table are neglected.

Independent validation

Every algebraic API has a fixed-value check, broadcasting or bound checks and invalid-input contracts. The hydrology tests also check that the water balance closes over a cycle, that Manning’s formula is Chézy’s with R**(1/6) / n and the hydraulic radius tends to the depth, that every section of Dupuit’s water table carries the Dupuit discharge, that the steady Richards profile carries a constant flux, is hydrostatic without flow and tends to unit gradient, that the Green-Ampt time inverts the rate and grows as the square of the depth early on, that the Nash hydrograph integrates to one and the storm hydrograph conserves volume and superposes, that Horton’s cumulative infiltration integrates the rate, that Theis tends to Cooper-Jacob and its log-cycle slope returns the transmissivity, that Muskingum coefficients sum to one and routing conserves volume while delaying and lowering the peak, that Penman-Monteith with zero surface resistance is Penman and saturation_slope is the derivative of the saturation curve, that TOPMODEL deficits average to the mean and baseflow falls by e per scaling depth, and that simulate_infiltration matches the Mein-Larson closed form after ponding. See Validation policy for the exact evidence and Numerical methods and sampled results for numerical conventions.