Hydrographs can be produced with different approaches. A straightforward approach is to use the storage cascade of Nash. The Nash cascade is a sequence of $N$ linear stores in series, each having a residence time of $K$ time units. The resulting mathematical form of the Nash-cascade is:

$$h(t) = (t/K)^{N-1}*exp\left(\frac{-t/K}{K*\Gamma(N)}\right)$$

where $h(t)$ is the unit hydrograph in terms of a unit discharge and $\Gamma$ is the gamma function $\gamma(N)=(N-1)!$ for integer values $N$. $N$ can also have rational values.

| Nash.R

n <- 7 # number of storages in sequence
k <- 0.3 # exponential decay in each store were k = 1/n with n time steps
fc = 2.5
x <- seq(0.1,3,by=0.05)
y <- 1/(k*gamma(n))*(x/k)^(n-1)*exp(-x/k)*fc # Nash Cascade
plot(x,y,col=2, xlab= " ", ylab="Q", xlim=range(0:3), ylim=range(0:3))
curve(1/(k *gamma(n))*(x/k )^(n-1)*exp(-x/k)*fc, add = TRUE, col="red")

First, you need to install the package berryfunctions. Then you can use functions to produce unit hydrographs.

| Berry

library(berryFunctions)
Time <- 0:100
plot(Time, unitHydrograph(n=3,  k=2, t=Time), type="l", las=1,
     main="Unit Hydrograph - linear storage cascade")
lines(Time, unitHydrograph(n=2,  k=1, t=Time), col=2)
lines(Time, unitHydrograph(n=5,k=5, t=Time), col=4)
text(c(12, 20, 50), c(0.1, 0.04, 0.025), c("n=2, k=3","n=2, k=8","n=5.5, k=8"),
     col=c(1,2,4), adj=0)

The time of concentration can be calculated with the formula of Kerby-Kirpich:

$$t_c = K \cdot {(L \cdot N)}^{0.467} \cdot S^{-0.235}$$

where: $t_{ov}$: overland flow time of concentration, in minutes $K$: a units conversion coefficient, in which K = 0.828 for traditional units and K = 1.44 for SI units $L$: the overland-flow length, in feet or meters as dictated by K $N$: a dimensionless retardance coefficient $S$: the dimensionless slope of terrain conveying the overland flow

| Kirpich.R

K = 1.44 # for SI units
L = 8000 # the overland-flow length, in meters as dictated by K
N = 3    # retardance coefficient
S = 0.010 # the dimensionless slope of terrain conveying the overland flow
 
tm = K*(L*N)^0.467*S^-0.235
tm
th = tm/60
th

The unit hydrograph can also be estimated using the Advection Dispersion equation.

$$Q_(t) = \frac{x}{\sqrt{4*pi*D*t^3}}*exp\left[-\frac{(x-v_c*t)^2}{4*D*t}\right]$$

with the dispersion coefficient as

$$ D = v_c*\alpha = v_c * f_D * x $$

with $Q(t)$: discharge in $[m^3/s]$ $x$: Distance in meters $D$: Dispersion coefficient $t$: time in hours $v_c$ vecolicity in $[m/s]$ $f_D$ Scale-Factor $[m]$

| Ade.R

c<-1
D<-0.5
x<-12
fc = 12
t<-seq(0,24,0.5)
u<-fc*x/(4*pi*D*t^3)^(1/2)*exp(-(x-c*t)^2/(4*D*t))
plot(t,u,xlab="time", ylab="Discharge",
  xlim=c(0, 24.0), ylim=c(0, 2.0),type="l", col="blue", axes=FALSE)
axis(1, seq(0,24.0,2))
axis(2, seq(0,2.0,0.5))
grid (NULL,NULL, lty = 3, col = "grey")

Make a table (array of the dimensionless unit hydrograph) that you get from Link to tabulated unit hydrograph

Calculate the concentration time, and multiply the dimensionless time with the concentration time. Then multiply the dimension-less flow (from 0 to 1) with the peak flow. This will give you the full hydrograph.

| UnitHydrograph.R

```{r uhdata}
ut = c(0.00,0.10,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90,1.00,1.10,1.20,1.30,1.40,1.50,1.60,1.70,1.80,1.90,2.00,2.20,2.40,2.60,2.80,3.00,3.20,3.40,3.60,3.80,4.00,4.50,5.00)
uq = c(0.00,0.03,0.10,0.19,0.31,0.47,0.66,0.82,0.93,0.99,1.00,0.99,0.93,0.86,0.78,0.68,0.56,0.46,0.39,0.33,0.28,0.21,0.15,0.11,0.08,0.06,0.04,0.03,0.02,0.02,0.01,0.01,0.00)
 
plot(ut,uq, type="l")
 
K = 1.44 # for SI units
L = 6500 # the overland-flow length, in meters as dictated by K
N = 1.00 # retardance coefficient
S = 0.01 # the dimensionless slope of terrain conveying the overland flow
qp= 1.41 # here the design flow is 1.41 m³/s
 
tm = K*(L*N)^0.467*S^-0.235 # concentration time in minutes
tm                          # show it
th = tm/60                  # concentration time in hours
th                          # show it
 
ur = ut*th                  # hydrograph time in real time units, here hours
qr = uq*qp                  # peak flow to be determined from hydrograph analysis
 
plot(ur,qr, type="l")