Hydrograph as a mathematical function. See https://tonyladson.wordpress.com/2015/07/21/hydrograph-as-a-function/

Hydro.R

Hydro <- function(tt, t.peak = 1, Qmin = 1, Qmax = 10, beta = 5) {
  Qmin + (Qmax - Qmin)*( (tt/t.peak) * (exp(1 - tt/t.peak)))^beta
  
}



tt.seq <- seq(0,4,0.01)

op <- par(oma = c(1, 2, 0, 0))
my.beta <- 5
plot(tt.seq, Hydro(tt.seq, beta = my.beta), type = 'l',
     las = 1, 
     xlab = 'Time',
     ylab = 'Discharge')
par(op)

op <- par(mfrow = c(1, 2))
my.beta <- 1
plot(tt.seq, Hydro(tt.seq, beta = my.beta), type = 'l',
     las = 1, 
     xlab = 'Time',
     ylab = 'Discharge')
legend('topright', legend = bquote(beta == .(my.beta)), bty = 'n')

my.beta <- 10
plot(tt.seq, Hydro(tt.seq, beta = my.beta), type = 'l',
     las = 1, 
     xlab = 'Time',
     ylab = '')
legend('topright', legend = bquote(beta == .(my.beta)), bty = 'n')
par(op)


# Inflection points
library(Deriv)

my.ex <- expression(Qmin + (Qmax - Qmin)*( (tt/t.peak) * (exp(1 - tt/t.peak)))^beta)

Deriv(my.ex, 'tt')  # First derivative

my.ex2 <- Deriv(Deriv(my.ex, 'tt'), 'tt') # Second derivative

# Convert second derivative to function

Hydro2D <- function(tt, t.peak, Qmin , Qmax , beta ){ eval( my.ex2[[1]])}

# To location the inflection points
# Find where Hydro2D is equal to zero

# Inflection point after peak
uniroot(Hydro2D, interval = c(2, 10), Qmin = 1, Qmax = 10, beta = 5, t.peak = 2 )

#$root
#[1] 2.894449

# check
2 * (1 + (1/sqrt(5)))
#[1] 2.894427 ok

# Inflection point before peak
uniroot(Hydro2D, interval = c(1, 2), Qmin = 1, Qmax = 10, beta = 5, t.peak = 2 )
# [1] 1.105573
2 * (1 - (1/sqrt(5)))
#[1] 1.105573 ok

# Plot the inflection points
op <- par(oma = c(1, 2, 0, 0))

plot(tt.seq, Hydro(tt.seq, Qmin = 1, Qmax = 10, beta = 5, t.peak = 2), type = 'l',
     las = 1, 
     xlab = 'Time',
     ylab = 'Discharge')

x <- uniroot(Hydro2D, interval = c(1, 2), Qmin = 1, Qmax = 10, beta = 5, t.peak = 2 )
points(x$root, Hydro(x$root, Qmin = 1, Qmax = 10, beta = 5, t.peak = 2 ), col = 'red')

x <- uniroot(Hydro2D, interval = c(2, 4), Qmin = 1, Qmax = 10, beta = 5, t.peak = 2 )
points(x$root, Hydro(x$root, Qmin = 1, Qmax = 10, beta = 5, t.peak = 2 ), col = 'red')

legend('topright','inflection points', pch = 1, col = 'red', bty = 'n', cex = 0.8)

par(op)