Traditional AHG fitting

The AHG relationships

The idea of hydraulic geometry was introduced by Leopold and Maddock (1953) and the concept of ‘at-a-station hydraulic geometry’ (AHG) provides three power laws that relate water-surface width (TW), average depth (Y), and average velocity (V) to the stream flow (Q) in a given river cross section:

\[ TW = a \cdot Q^b \] \[ Y = c \cdot Q^f \] \[ V = k \cdot Q^m \]

A set of AHG relations applies to within-bank flows at a specific cross section and assumes the channel characteristics do not change sensibly with discharge. Under these assumptions, continuity dictates that no water is gained or lost in the system such that:

\[ Q = TW·Y·V \]

and therefore:

\[ b + f + m = a·c·k = 1 \]

In practice, the enforcement of continuity is determined by the quality of the data, and is not enforced in the estimation of the AHG relationships.

Traditonal Fitting Technique

While the AHG power laws provide an elegant theory for representing a mass-constrained hydraulic system, the data available to determine these relations are rarely as clean. In most cases, hydraulic data is reported as a table(s) relating values of Q to values of one or more hydraulic variables. Thus, the relationships require the power law relationships to be fit,

AHG relations are historically fit using OLS on the logarithmic transformation of the variables, where the exponential of the intercept provides the power law coefficient, and the slope of the linear model provides the exponent. Several studies have highlighted the potential problem of using this approach given that when the log transformed values are back-transformed, the estimates are effectively medians, and not means of the estimates, resulting in a general low bias.

Example:

data <- select(nwis, date, Q = Q_cms, Y = Y_m, TW = TW_m, V = V_ms)

fit <- lm(log(data$Y) ~ log(data$Q))
  
(coef <- exp(fit$coefficients[1]))
#> (Intercept) 
#>   0.2018376

(exp  <- fit$coefficients[2])
#> log(data$Q) 
#>   0.4794202

Full Fits

ols <- function(X, Y, name = NA){
  fit <- lm(log(Y) ~ log(X))
  data.frame(coef = exp(fit$coefficients[1]),
             exp  = fit$coefficients[2],
             name = name,
             row.names = NULL)
}

(x <- bind_rows(
        ols(data$Q, data$Y, "Y"),
        ols(data$Q, data$TW, "TW"),
        ols(data$Q, data$V, "V")))
#>         coef       exp name
#> 1  0.2018376 0.4794202    Y
#> 2 22.6153841 0.1145234   TW
#> 3  0.2178592 0.4056932    V

Continuity check

sum(x$exp) 
#> [1] 0.9996368
prod(x$coef)
#> [1] 0.9944476

System Error