R-based model using `AGWRC=f(AGWS)` and HSPF equations

[rburghol]

• Daily Timestep
• AGWS computed each day as $AGWS_{t+1} = AGWS_t - AGWO$
• Must compute $AGWRC = f(AGWS)$ -- see below
• Using Q to get AGWRC regression equations, then plugging in agws

Equations for $AGWRC = f(AGWS)$

• Start with regression equation $AGWRC = f(Q)$
• Substitute in $Q(in) = AGWS * (1.0 - AGWRC)$
• Must derive AGWS from Q using equation above, and since AGWS is in inches may as well make that conversion

Algorithm

• Notes:
  • Final model must be in inches, AGWS = inches, AGWO = inches/day
  • First life is converting
• To run the model: start with initial AGWS and function of the relationship between AGWRC=f(AGWS)
• loop thru timesteps (could loop thru elements of a dataframe that had a time column, or just a specified number of steps)
  • $AGWRC = f(AGWS)$
  • $AGWO = f(AGWS, AGWRC)$
  • $AGWS = AGWS - AGWO$
  • stash each variable in timeseries dataframe

Tasks

• Function to convert Q to inches (Ilona)
• Create function for regression and storing regression values (Check Units) (Will)
• Redo regression with Q in inches (Will)
• Case study to check that formulas work with real/updated numbers
• Test uniroot speed (Ben)
  • Change uniroot limits (0.75, 0.999)

[rburghol]

@ilonah22 @willprokopik @benhdye
@julieshortridge
My suggestions for today:

• 1 person take the lead on the iteration timestep code loop.
  • include a function to calculate AGWRC = f(AGWS), which includes a default AGWRC parameter like calc_agwrc = function(AGWS, m, x, agwrc_default). Initially, this function should just return the agwrc_default back, but when the team completes the transformation from Q to AGWS that code can be substituted in.
  • test your code output for some sample starting AGWS values, and compare them to our existing projection equations/BF event flows.
• 2 people work as a team to handle the mathematics converting the AGWRC = f(Q), into AGWRC = f(AGWS).
  • I would simply start with replacing Q by the AGWS * (1.0 - AGWRC) term. Then solve for AGWRC.
  • if AGWRC gets zeroed out, well then you may have to come up with something more clever :).
  • Note: you may need to re-calculate your regression equation using flow as expressed in watershed inches, otherwise your coefficient/slope/intercept will probably not work.
  • if you are successful in solving for AGWRC, then verify your code by doing some manual calculations and see what happens.
  • Update the calc_agwrc = function(AGWS, m, x, agwrc_default) to use your new relatiinship
• Once both of the above are completed, put them together.

[ilonah22]

Work from today's meeting:

Sorry for putting this up a bit late, I had to rewrite my work so it was legible.

@rburghol, We initially tried to follow the outline in your previous comment, but found that, to the best of our knowledge, it won't be possible to substitute AGWS for Q without already having an AGWRC value. Also, @willprokopik and @benhdye, I noticed an error that I made earlier, but fixed it so the equations may look a little different than what you remember.

We determined that AGWO and Q should be interchangeable for our purposes. We were able to get AGWRC=f(AGWS), but that required using multiple days of storage data, which may not be helpful if the final goal is one AGWS input. We hit a roadblock and weren't really sure where to go from from here. So if you or @julieshortridge have a chance to look at this before Tuesday, hopefully we can figure something out, unless there was some other method you already had in mind.

[rburghol]

@ilonah22 @benhdye @willprokopik - thanks for starting this. Here is how I see it:

• To make for easier reading I will use $C = AGWRC; S = AGWS$
• a: $C = mQ + b$ ; From regression we get this
• b: $Q = S * (1 - C)$ ; from HSPF
• $C = m * (S * (1 - C)) + b$; subbing b into a
• → $C = m * (S - SC) + b$
• → $C = mS - mSC + b$
• → $C + mSC = mS + b$
• → $C * (1 + mS) = mS + b$
• $C = (mS + b) / (1 + mS)$

Now, this isn't quite as simple as all that since I think the equation is actually $C = m*ln(Q) + b$, which will make the math harder, but the principle should remain the same.

Did I get this right @julieshortridge

[willprokopik]

@rburghol I took a look at the work @ilonah22 provided, and I also looked into using this equation: $C = m \ln(Q) + b$ and here is what I found:

From Ilona's work I plugged $S$ = 2.0; $m$ = 0.1; $b$ = 0.8 into:
$C = \frac{mS + b}{1 + mS}$

$Q = S(1 - C)$

$C_{\text{check}} = mQ + b$
and I got $C_{\text{check}}$ = 0.833333333333

I then rewrote this $C = m \ln(Q) + b$ so that it could be solved: $g(C) = C - \left[m \ln\big(S(1 - C)\big) + b\right] = 0$ (I think I did this correctly? Let me know if this doesn't look right). I created a code chunk that could test out the same values that I plugged into Ilona's work ($S$ = 2.0; $m$ = 0.1; $b$ = 0.8):

solve_agwrc_log <- function(S, m, b) {
  
  g <- function(C) {
    C - (m * log(S * (1 - C)) + b)
  }
  
  res <- tryCatch(
    uniroot(g, lower = 0.001, upper = 0.999)$root,
    error = function(e) NA_real_
  )
  
  return(res)
}

and it returned 0.7360972
This suggests that the two types of calculations produce meaningfully different values of $C$ for the same $S$... $Q = S(1 − C) = 2.0(1 − 0.736)$ = 0.528
vs
$Q = S(1 - C) = 2.0(1 - 0.833)$ = 0.333

[rburghol]

Thanks @willprokopik -- that is sensible. My first question is, was I missing something with my algebra (copied from my above post):

• a: $C = mQ + b$ ; From regression we get this
• b: $Q = S * (1 - C)$ ; from HSPF
• $C = m * (S * (1 - C)) + b$; subbing b into a
• → $C = m * (S - SC) + b$
• → $C = mS -mSC + b$
• → $C + mSC = mS + b$
• → $C * (1 + mS) = mS + b$
• $C = (mS + b) / (1 + mS)$

And then, can you re-do my above (@ilonah22 @benhdye @willprokopik ), but using the correct equation for i.e., $C = m*ln(Q) + b$

• a: $C = m * ln(Q) + b$ ; From regression we get this
• b: $ln(Q) = ln(S * (1 - C))$ ; from HSPF
• $C = m * ln(S * (1 - C)) + b$; subbing b into a
• → $C = m * ln((S - SC)) + b$
• → $C = m * (ln(S) + ln(1 - C)) + b$
• → $C = m * ln(S) + m*ln(1 - C) + b$
• → $C = m * ln(S) + m*ln(1 - C) + b$
• This seems hard, options:
  - (maybe) Google recs Lambert W function
  - Assume the baseflow variation is linear instead of log-linear
  - This will require an analysis of inflection point (i.e., where does linear make a reasonable fit)
  - Use a lookup table : @COBrogan and I have talked about this, generally as a less than desirable approach, but, it may actually be OK -- the real questions in my mind are "does it mimic behavior in the real world?" and/or "can be it be implemented in a computationally efficient manner?"
  - See below for approximation of $ln(1 - C) \approx -C$
• Or, since since C almost always < 0.1, when $C &lt;&lt; 1$ supposedly Taylor series expansion says $ln(1 - C)$ can be approximated as $ln(1 - C) \approx -C$
• → $C - m*ln(1 - C) = m * ln(S) + b$
• → $C + mC = m * ln(S) + b$;
• → $C (1 + m) = m * ln(S) + b$
• → $C = (m * ln(S) + b) / (1 + m) $
• -- note @COBrogan I got this from the google LLM so who knows if it's true. But, shown below, it's close at the range of C = 0.05 to 0.1, which while it still may be unreasonable for our use case, since 1/1000 th difference in AGWRC could lead to different baseflow outcomes, could be worth a chance.

log(0.9)
[1] -0.1053605
 log(0.95)
[1] -0.05129329

[benhdye]

I think we can use uniroot like Will did to solve in the log case too

solve_agwrc_log <- function(S, m, b) {
  g <- function(C) {
    C - (m * log(S * (1 - C)) + b)
  }
  uniroot(g, lower = 1e-6, upper = 1 - 1e-6, tol = 1e-10)$root
}

S <- 2.0
m <- 0.1
b <- 0.8

C_solution <- solve_agwrc_log(S, m, b)
Q      <- S * (1 - C_solution)

cat(sprintf("C     = %.6f\n", C_solution))
cat(sprintf("Q     = %.6f\n", Q))
cat(sprintf("Check = %.6f\n", m * log(Q) + b))

> cat(sprintf("C     = %.6f\n", C_solution))
C     = 0.736097
> cat(sprintf("Q     = %.6f\n", Q))
Q     = 0.527805
> cat(sprintf("Check = %.6f\n", m * log(Q) + b))
Check = 0.736097

Check matches the calculated value for C so Uniroot found a calculated solution for C even in the log.

[rburghol]

Hey @willprokopik @benhdye -- appreciate the inquiry on this. Perhaps uniroot is good, though to be sure we may not have access to that in fortran/python, but if it's worth the effort we can find a way.

That said, the 0.8 case that you have proposed doesn't offer us actual validation. You need to solve a case from our actual AGWRC regression equations, and demonstrate that whatever proposed solution you come up with yields the correct numbers.

[julieshortridge]

Hi everyone, sorry to be late to this- I was tied up with other projects last week but have tried to get caught up this morning. Here are my comments on some questions/issues that have come up so far:

• Rob, I don't think you're mistaken anywhere in your algebra from the simplified regression (C=mQ+b) or the actual regression (C=m*ln(Q)+b). I checked Will's equations and don't see any issues there either.
• From Will: "This suggests that the two types of calculations produce meaningfully different values of C for the same S": I think this is to be expected because one is solving for Q using the logarithmic regression equation and one is not
• I don't think there is a way to algebraically solve for C when there are terms both within and outside the ln term. Given the complexity of both the Lambert W and Taylor series expansion (and the fact that both require making certain assumptions about the form and value of C), I think that the analyst's suggestion to use uniroot to solve for C is a logical approach.
• If we ultimately have to utilize a uniroot type solver in a non-R language, I think this should be feasible. The function itself is not complicated - it's just searching across an interval to find the value of C where C - m*ln(S(1-C)) + b = 0. If this ultimately proved necessary, I think the main issue would be constraining the input values of the function (m, S, and b) so it stays solvable.
• One comment that tripped me up from Rob was was "since C almost always < 0.1". Based on our two starting equations, C = m*ln(Q) + b AND Q = S(1-C), isn't it the opposite: that 1-C is generally << 1.0?
• I agree with the comment that the b=0.8 does not provide good validation - shouldn't we be using our actual regression coefficients here:
  https://github.com/HARPgroup/baseflow_storage/blob/ih_data/data/QC_coefficients.csv

Let me know if I'm missing something.

[rburghol]

Thanks for the review @julieshortridge and good catch on my last round of math:

• 1 - C << 1.0 I was confusing myself thinking of the C as k, which is actually 1 - C, but that's not the same thing. Drat.