- Normal depth of a canal with Manning's equation
- Stage-discharge rating curve
- Tidal constituents from a sea level record
- Water density and pressure at depth
- 100-year flood with a gamma law
- Theis drawdown of a pumping test and the Cooper-Jacob check
- Green-Ampt infiltration and ponding time
- Arrival of a contaminant at a well with the Ogata-Banks solution
- Wavelength and shoaling of a swell on a beach
- Drawdown of a well near a river with images and complex potential
Normal depth of a canal with Manning's equation
For a trapezoidal canal, compute the discharge at 1.5 m depth, solve for the normal depth of the 40 m3/s design flood and check the Froude number.
#1. Trapezoidal canal: bottom width 4 m, side slope 2H:1V, Manning n 0.030, bed slope 0.0008; -> "[Q]" "(4+2*y)*y*((4+2*y)*y/(4+2*y*sqrt(5)))^(2/3)*sqrt(0.0008)/0.030"; #2. Discharge at a flow depth of 1.5 m; -> "[Q15]" (math eval [Q] "{\"y\": 1.5}"); #3. Normal depth for the design flood of 40 m3/s; -> "[yn]" (num solve (concat [Q] "-40") "y" 0.1 5); #4. Velocity and Froude number at that depth, hydraulic depth A/T; -> "[A]" (math calc "(4+2*[yn])*[yn]"); -> "[Fr]" (math calc "40/[A]/sqrt(_g0*[A]/(4+4*[yn]))"); concat "Q at 1.5 m = " (math round [Q15] 2) " m3/s, normal depth for 40 m3/s = " (math round [yn] 3) " m, V = " (math round (math calc "40/[A]") 2) " m/s, Froude = " (math round [Fr] 3)
Q at 1.5 m = 9.77 m3/s, normal depth for 40 m3/s = 2.978 m, V = 1.35 m/s, Froude = 0.316
- Get the critical depth by solving Fr = 1
- Draw the rating curve with
num table
Stage-discharge rating curve
Fit the power law Q = a (h - h0)^b to the gaugings of a river station, then read the discharge of a flood peak above the measured range.
#1. Gaugings of a river station: stage h in m, measured discharge Q in m3/s; -> "[h]" "[0.45, 0.62, 0.80, 1.05, 1.30, 1.62, 2.10]"; -> "[q]" "[0.76, 2.22, 4.30, 8.28, 13.17, 19.78, 32.82]"; #2. Rating curve Q = a*h-h0^b fitted by non-linear least squares; -> "[fit]" (stat fit_model "a*(x-h0)^b" "x" [h] [q] "{\"a\": 10, \"h0\": 0.2, \"b\": 1.5}"); json load "rx_rating" [fit]; #3. Discharge of the flood peak read at 2.85 m on the staff gauge; -> "[Qp]" (math eval (json select "rx_rating" /expression) "{\"x\": 2.85}"); concat "a = " (math round (json select "rx_rating" /parameters/a) 3) ", h0 = " (math round (json select "rx_rating" /parameters/h0) 3) " m, b = " (math round (json select "rx_rating" /parameters/b) 3) ", r2 = " (math round (json select "rx_rating" /r2) 4) ", Q at 2.85 m = " (math round [Qp] 1) " m3/s"
a = 12.425, h0 = 0.275 m, b = 1.612, r2 = 0.9998, Q at 2.85 m = 57.1 m3/s
- Weight the gaugings by their uncertainty
- Integrate a hydrograph with
num interpolate_integralto get the flood volume
Tidal constituents from a sea level record
Take 30 days of hourly sea level, compute its amplitude spectrum, pick the tidal lines and classify the tide with the form factor.
#1. 30 days of hourly sea level in m: mean level and four constituents M2, S2, K1, O1; -> "[eta]" "2.5+1.2*cos(2*pi*t/12.4206)+0.4*cos(2*pi*t/12+0.7)+0.15*cos(2*pi*t/23.9345+1.1)+0.10*cos(2*pi*t/25.8193+2.3)"; -> "[x]" (matrix col (num table [eta] "t" 0 719 1) 2); #2. Amplitude spectrum, sampling 24 values per day; -> "[sp]" (dsp spectrum [x] 24); json load "rx_tide" [sp]; #3. Spectral lines higher than 5 cm; -> "[pk]" (dsp peaks (json select "rx_tide" /amplitude) "{\"height\": 0.05}"); json load "rx_peaks" [pk]; #4. Form factor K1+O1 / M2+S2 to classify the tide; -> "[F]" (math calc (concat "(" (json select "rx_peaks" /values[0]) "+" (json select "rx_peaks" /values[1]) ")/(" (json select "rx_peaks" /values[2]) "+" (json select "rx_peaks" /values[3]) ")")); concat "lines in cycles/day: " (vector scale (json select "rx_peaks" /indices) (math calc "24/720")) ", amplitudes: " (json select "rx_peaks" /values) ", form factor = " (math round [F] 3)
lines in cycles/day: [0.9333333333333333,1,1.9333333333333333,2], amplitudes: [ 0.10338664735793467, 0.1525429251899293, 1.1976164507455016, 0.41290282434961517 ], form factor = 0.159
- Remove the surge with
dsp filterbefore the analysis - Fit the exact constituent periods with
stat fit_model
Water density and pressure at depth
Interpolate the density of water at the lake temperature from a table, then compute the pressure at 35 m depth and the mass of the reservoir.
#1. Density of pure water in kg/m3 from 0 to 40 degC; -> "[T]" "[0, 4, 10, 15, 20, 25, 30, 40]"; -> "[rho]" "[999.84, 999.97, 999.70, 999.10, 998.21, 997.05, 995.65, 992.22]"; #2. Cubic spline at the water temperature of the lake, 17.3 degC; -> "[r]" (num interpolate [T] [rho] 17.3 "cubic"); #3. Absolute pressure at 35 m depth and mass of the 1.2 hm3 reservoir; -> "[p]" (math calc "[r]*_g0*35+_atm"); concat "rho = " (math round [r] 2) " kg/m3, pressure at 35 m = " (math round (unit convert [p] "Pa" "bar") 3) " bar = " (math round (unit convert [p] "Pa" "psi") 1) " psi, mass = " (math round (math calc "[r]*1.2e6/1e9") 4) " Mt"
rho = 998.73 kg/m3, pressure at 35 m = 4.441 bar = 64.4 psi, mass = 1.1985 Mt
- Add the salinity with the UNESCO equation for sea water
- Find the temperature of maximum density
100-year flood with a gamma law
Twenty annual maximum discharges are fitted by a gamma law by maximum likelihood: the shape solves ln k - digamma(k) = ln(mean) - ln(geometric mean), its standard error uses the trigamma function, and the 100-year flood is the 0.99 quantile from the incomplete gamma function.
#1. Annual maximum discharge of a river over 20 years, m3/s; -> "[Q]" "[412, 535, 298, 640, 377, 455, 820, 510, 365, 590, 433, 702, 285, 486, 560, 395, 678, 470, 350, 615]"; -> "[n]" 20; -> "[m]" (stat mean [Q]); #2. Gamma law by maximum likelihood: ln k - digamma of k = ln mean - ln geometric mean; -> "[s]" (math calc (concat "ln([m])-ln(" (stat geo_mean [Q]) ")")); -> "[k]" (num solve (concat "ln(k)-digamma(k)-" [s]) "k" 0.5 200); -> "[theta]" (math calc "[m]/[k]"); #3. Standard error of k with the trigamma function; -> "[sek]" (math calc "sqrt([k]/([n]*([k]*polygamma(1, [k])-1)))"); #4. Log-likelihood with lgamma, and the 100 year flood: P of k, q/theta = 0.99; -> "[ll]" (math calc (concat "([k]-1)*[n]*ln(" (stat geo_mean [Q]) ")-[n]*[m]/[theta]-[n]*[k]*ln([theta])-[n]*lgamma([k])")); -> "[q100]" (num solve (concat "gamma_inc(" [k] ", q/" [theta] ")-0.99") "q" 500 3000); concat "k = " (math round [k] 2) " +- " (math round [sek] 2) ", theta = " (math round [theta] 2) " m3/s, log-likelihood = " (math round [ll] 2) ", Q100 = " (math round [q100] 0) " m3/s"
k = 12.91 +- 4.03, theta = 38.63 m3/s, log-likelihood = -126.51, Q100 = 877 m3/s
- Compare with a Gumbel fit and the same return period.
- The log-likelihood allows to compare laws with the AIC.
Theis drawdown of a pumping test and the Cooper-Jacob check
In a confined aquifer pumped at 20 L/s, the script computes the Theis argument u, the well function W(u) as an improper integral (the exponential integral E1), the drawdown at a piezometer after one day, and compares with the Cooper-Jacob straight line used on semi-log field plots.
#1. Pumping test in a confined aquifer: Q = 20 L/s, T = 5e-3 m2/s, S = 2e-4, piezometer at 50 m, after 1 day; -> "[Q]" 0.020; -> "[T]" 5e-3; -> "[S]" 2e-4; -> "[r]" 50; -> "[t]" 86400; #2. Theis argument u = r^2 S / 4 T t; -> "[u]" (math calc "[r]^2*[S]/(4*[T]*[t])"); #3. Well function W u = exponential integral, from u to infinity of exp -x / x; -> "[W]" (num integrate "exp(-x)/x" "x" [u] "inf"); -> "[s]" (math calc "[Q]/(4*pi*[T])*[W]"); #4. Cooper-Jacob straight line approximation, valid when u is small; -> "[sj]" (math calc (concat "[Q]/(4*pi*[T])*" (math log (math calc "2.25*[T]*[t]/([r]^2*[S])")))); concat "u = " (math round_sig [u] 3) ", W(u) = " (math round [W] 4) ", Theis drawdown = " (math round [s] 3) " m, Cooper-Jacob = " (math round [sj] 3) " m, gap = " (math round (math calc "100*abs([sj]-[s])/[s]") 2) " %, drawdown per log cycle = " (math round (math calc "2.303*[Q]/(4*pi*[T])") 3) " m"
u = 0.000289, W(u) = 7.5709, Theis drawdown = 2.41 m, Cooper-Jacob = 2.41 m, gap = 0.02 %, drawdown per log cycle = 0.733 m
- Fit T and S to measured drawdowns with
stat fit_modelon the Cooper-Jacob line. - Add the recovery after the pump stops with the superposition of an image well.
Green-Ampt infiltration and ponding time
For a sandy loam with the Rawls Green-Ampt parameters, the script solves the implicit equation of the cumulative infiltration after two hours, checks it with log1p, gets the infiltration capacity and the ponding time under a 3 cm/h storm.
#1. Green-Ampt parameters of a sandy loam: Ks = 1.09 cm/h, suction 11.01 cm, moisture deficit 0.288; -> "[K]" 1.09; -> "[psi]" 11.01; -> "[dth]" 0.288; -> "[pd]" (math calc "[psi]*[dth]"); #2. Cumulative infiltration after 2 h of ponding: F - psi dtheta ln 1 + F/psi dtheta = K t; -> "[F]" (num solve (concat "x-" [pd] "*ln(1+x/" [pd] ")-" [K] "*2") "x" 0.01 50); #3. Check the equation with log1p, accurate for small F, and the infiltration rate; -> "[check]" (math calc (concat "[F]-[pd]*" (math log1p (math calc "[F]/[pd]")) "-[K]*2")); -> "[f]" (math calc "[K]*(1+[pd]/[F])"); #4. Ponding time under a rain of 3 cm/h; -> "[Fp]" (math calc "[K]*[pd]/(3-[K])"); concat "F(2 h) = " (math round [F] 3) " cm, residual = " (math round_sig [check] 2) ", infiltration rate = " (math round [f] 3) " cm/h, ponding after " (math round (math calc "[Fp]/3*60") 1) " min of a 3 cm/h rain with " (math round [Fp] 3) " cm infiltrated"
F(2 h) = 5.293 cm, residual = -0.00000000000000044, infiltration rate = 1.743 cm/h, ponding after 36.2 min of a 3 cm/h rain with 1.81 cm infiltrated
- Loop on the time to build the infiltration curve and plot it with
plot xy. - Compare with the Horton law fitted with
stat fit_model.
Arrival of a contaminant at a well with the Ogata-Banks solution
A dissolved pollutant enters a sandy aquifer. The script evaluates the Ogata-Banks solution of the advection-dispersion equation at a well 100 m away, then the arrival time of 10 percent of the source concentration with the inverse error function, compared with pure advection.
#1. Contaminant entering a sandy aquifer at x = 0: pore velocity 0.5 m/d, longitudinal dispersivity 5 m, well at 100 m; -> "[v]" 0.5; -> "[D]" (math calc "5*[v]"); -> "[x]" 100; #2. Ogata-Banks solution at t = 200 d, both terms; -> "[t]" 200; -> "[C]" (math calc (concat "0.5*(" (math erfc (math calc "([x]-[v]*[t])/(2*sqrt([D]*[t]))")) "+exp([v]*[x]/[D])*" (math erfc (math calc "([x]+[v]*[t])/(2*sqrt([D]*[t]))")) ")")); #3. Arrival of 10 percent of the source concentration: erfc y = 0.2, y = erf_inv 0.8; -> "[y]" (math erf_inv 0.8); -> "[t10]" (num solve (concat "(" [x] "-" [v] "*t)/(2*sqrt(" [D] "*t))-" [y]) "t" 1 400); #4. Advection only would give x/v; concat "C/C0 at 200 d = " (math round [C] 4) ", erfc-1(0.2) = " (math round [y] 5) ", 10 percent arrival after " (math round [t10] 1) " d, 50 percent after about " (math round (math calc "[x]/[v]") 0) " d, so dispersion brings the front " (math round (math calc "[x]/[v]-[t10]") 1) " d earlier"
C/C0 at 200 d = 0.5616, erfc-1(0.2) = 0.90619, 10 percent arrival after 133.7 d, 50 percent after about 200 d, so dispersion brings the front 66.3 d earlier
- Add the retardation factor R = 1 + rho Kd / n: divide v and D by R.
- Draw the breakthrough curve with
plot functionfrom 0 to 400 days.
Wavelength and shoaling of a swell on a beach
A swell of 8 s reaches a depth of 5 m. The script solves the linear dispersion relation with tanh for the wave number, gets the wavelength, the celerity and the group velocity with sinh, the shoaling coefficient and the wave height, and classifies the depth.
#1. Swell of period 8 s and height 1.5 m in deep water, reaching a beach at 5 m depth; -> "[T]" 8; -> "[h]" 5; -> "[H0]" 1.5; -> "[g]" 9.80665; -> "[w]" (math calc "2*pi/[T]"); -> "[L0]" (math calc "[g]*[T]^2/(2*pi)"); #2. Dispersion relation w^2 = g k tanh k h solved for the wave number; -> "[k]" (num solve (concat [w] "^2-" [g] "*x*tanh(x*" [h] ")") "x" 0.01 2); #3. Celerity, group velocity with sinh, shoaling coefficient; -> "[c]" (math calc (concat [w] "/[k]")); -> "[n]" (math calc (concat "0.5*(1+2*[k]*[h]/" (math sinh (math calc "2*[k]*[h]")) ")")); -> "[Ks]" (math calc "sqrt(([L0]/[T]/2)/([n]*[c]))"); #4. Depth class from kh and its tanh; concat "L0 = " (math round [L0] 2) " m, L = " (math round (math calc "2*pi/[k]") 2) " m, tanh(kh) = " (math round (math tanh (math calc "[k]*[h]")) 4) ", c = " (math round [c] 3) " m/s, cg = " (math round (math calc "[n]*[c]") 3) " m/s, shoaling Ks = " (math round [Ks] 4) ", H = " (math round (math calc "[Ks]*[H0]") 3) " m, " (if (< (math calc "[k]*[h]") (math calc "pi/10")) {"shallow water"} {"intermediate depth"})
L0 = 99.89 m, L = 53.07 m, tanh(kh) = 0.5313, c = 6.634 m/s, cg = 5.969 m/s, shoaling Ks = 1.0227, H = 1.534 m, intermediate depth
- Add the refraction with Snell law on the wave angle.
- Find the breaking depth where H = 0.78 h with
num solve.
Drawdown of a well near a river with images and complex potential
A well pumps 100 m from a river that keeps a constant head. The method of images adds a recharge well on the other side; the complex potential Q/2pi log((z - zw)/(z - zi)) gives in one complex logarithm the drawdown (real part) and the stream function (imaginary part) at a piezometer, and the drawdown in the well.
#1. Well 100 m from a straight river, steady pumping Q = 15 L/s, transmissivity 4e-3 m2/s, well radius 0.15 m; -> "[Q]" 0.015; -> "[T]" 4e-3; -> "[zw]" "100"; -> "[zi]" "-100"; #2. The river is a constant head: an image recharge well on the other side, complex potential Q/2pi log z - zw / z - zi; -> "[P]" "60+40i"; -> "[Om]" (complex log (complex div (complex sub [P] [zw]) (complex sub [P] [zi]))); #3. Drawdown at the point P: real part, stream function: imaginary part; -> "[s]" (math calc (concat "-[Q]/(2*pi*[T])*" (complex re [Om]))); -> "[psi]" (math calc (concat "[Q]/(2*pi)*" (complex im [Om]))); #4. Drawdown at the well screen: ln 2d / rw; -> "[sw]" (math calc "[Q]/(2*pi*[T])*ln(200/0.15)"); concat "log ratio = " [Om] ", drawdown at P = " (math round [s] 3) " m, stream function = " (math round_sig [psi] 4) " m3/s, drawdown in the well = " (math round [sw] 3) " m"
log ratio = -1.0700330817481354+2.1112158270654806i, drawdown at P = 0.639 m, stream function = 0.005040 m3/s, drawdown in the well = 4.294 m
- Map the drawdown on a grid with
num tableand draw it withplot function. - Find the stagnation point on the axis: the capture zone of the well.
