Profiles / Space, Earth & energy
HO

Hydrologist & oceanographer

Open channel flow, rating curves, tidal spectra and water properties.

10 recipes on this page · open in the interactive view

Recipes
  1. Normal depth of a canal with Manning's equation
  2. Stage-discharge rating curve
  3. Tidal constituents from a sea level record
  4. Water density and pressure at depth
  5. 100-year flood with a gamma law
  6. Theis drawdown of a pumping test and the Cooper-Jacob check
  7. Green-Ampt infiltration and ponding time
  8. Arrival of a contaminant at a well with the Ogata-Banks solution
  9. Wavelength and shoaling of a swell on a beach
  10. Drawdown of a well near a river with images and complex potential
executed on MentDB

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.

Pipeline
Script
admin
#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)
Result
mentdb
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
Go further
  • Get the critical depth by solving Fr = 1
  • Draw the rating curve with num table
executed on MentDB

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.

Pipeline
Script
admin
#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"
Result
mentdb
a = 12.425, h0 = 0.275 m, b = 1.612, r2 = 0.9998, Q at 2.85 m = 57.1 m3/s
Go further
executed on MentDB

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.

Pipeline
Script
admin
#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)
Result
mentdb
lines in cycles/day: [0.9333333333333333,1,1.9333333333333333,2], amplitudes: [
  0.10338664735793467,
  0.1525429251899293,
  1.1976164507455016,
  0.41290282434961517
], form factor = 0.159
Go further
executed on MentDB

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.

Pipeline
Script
admin
#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"
Result
mentdb
rho = 998.73 kg/m3, pressure at 35 m = 4.441 bar = 64.4 psi, mass = 1.1985 Mt
Go further
  • Add the salinity with the UNESCO equation for sea water
  • Find the temperature of maximum density
executed on MentDB

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.

Pipeline
Script
admin
#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"
Result
mentdb
k = 12.91 +- 4.03, theta = 38.63 m3/s, log-likelihood = -126.51, Q100 = 877 m3/s
Go further
  • Compare with a Gumbel fit and the same return period.
  • The log-likelihood allows to compare laws with the AIC.
executed on MentDB

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.

Pipeline
Script
admin
#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"
Result
mentdb
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
Go further
  • Fit T and S to measured drawdowns with stat fit_model on the Cooper-Jacob line.
  • Add the recovery after the pump stops with the superposition of an image well.
executed on MentDB

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.

Pipeline
Script
admin
#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"
Result
mentdb
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
Go further
  • 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.
executed on MentDB

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.

Pipeline
Script
admin
#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"
Result
mentdb
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
Go further
  • Add the retardation factor R = 1 + rho Kd / n: divide v and D by R.
  • Draw the breakthrough curve with plot function from 0 to 400 days.
executed on MentDB

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.

Pipeline
Script
admin
#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"})
Result
mentdb
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
Go further
  • Add the refraction with Snell law on the wave angle.
  • Find the breaking depth where H = 0.78 h with num solve.
executed on MentDB

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.

Pipeline
Script
admin
#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"
Result
mentdb
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
Go further
  • Map the drawdown on a grid with num table and draw it with plot function.
  • Find the stagnation point on the axis: the capture zone of the well.

Also in Space, Earth & energy

AS
Astronomer
Kepler laws, parallax distances, magnitudes and stellar luminosity, computed with units and CODATA constants.
10 recipes
SM
Space mission engineer
Orbits, transfers and rockets: velocities, periods and delta-v budgets in a few lines.
10 recipes
AE
Aeronautical engineer
Standard atmosphere, lift and drag polar, glide performance, with knots, feet and nautical miles.
10 recipes
ME
Meteorologist
Humidity, heat and cold indices, sea level pressure and rainfall extremes.
10 recipes
CE
Climate & environment scientist
CO2 trends, temperature anomalies and carbon footprints in CO2 equivalent.
10 recipes
GS
Geologist & seismologist
Seismograms, magnitudes, earthquake statistics and radiometric ages.
10 recipes
RE
Renewable energy engineer
Solar and wind yield, battery sizing and the levelized cost of electricity.
10 recipes