mcolvin / Bluff-Lake-Project

Bluff Lake Project Repository
0 stars 0 forks source link

DO Model #7

Closed VictoriaStarnes closed 4 years ago

VictoriaStarnes commented 4 years ago

Bluff-Lake-Project/_analysis/DOModel.R

I have been trying to copy the equation from Dr. Miranda's 2001 paper on Eagle Lake (emailed as an attachment). The code should only be in lines 34-44. I think I've either misread the equation, or I've made an error in the code. However, I can't figure out what I've done wrong. Would you care to take a look at the paper and my code? Thanks!

mcolvin commented 4 years ago

@VictoriaStarnes not sure if it emailed this to you last night without tagging you in it.

The code threw an error when I ran it

For respiration: image

For DO Model 2

image

It would be helpful to annotate your code more fully. See line 43 where I added an annotation.

You are going to need to dust off the old calc/precalc book. Equation 1 is an integral that sums the area under the curve for dDO/dt (equation 2 in Miranda et al 2001) for a dt=1 minute. Take a look at the integrate function in R, something like eq_2<-function(DO_t){DO_t-SR*Z^-1 + WR +DE*Z^-1} Or it might look like this depending on how integrate takes the state variable DO eq_2<-function(DO_t){-SR*Z^-1 + WR +DE*Z^-1}

Might mess with both and feed them through the integrate function below.

do_dawn=do_dusk-integrate(eq_2,dusk,dawn)

That will do a quick simple numeric integration of the DO dynamics from dusk to dawn, rather than multiplying the function by 600 minutes.

I will let you mess with it a wee bit as a challenge to push your learning, if you are still spinning your wheels by Thursday let me know and I will work through it with you.

VictoriaStarnes commented 4 years ago

You can ignore the respiration section. It has now been deleted. I have a separate script for water column respiration.

I pulled the code from GitHub and did not see any of the annotations that you made. However, I added more descriptions to my code.

I can't get the integrate function to work. I've worked on it for a few hours (4) now and I don't think that I will be able to figure it out. No matter what I try I always get the errors "the integral is probably divergent" or "non-numeric argument to binary operator." Searches for these errors have not been helpful. I would love to walk through it whenever you have time. Thanks

mcolvin commented 4 years ago

@VictoriaStarnes Here is some standalone code that seems to give reasonable results

DO_dusk<-10 
parameters <- c( 
    tempC = 25, 
    Z = 2) # depth 
DO_fun<-function(t,x,parms)
    {
    DO<-x
    tempC<-parms["tempC"]
    Z<-parms["Z"]
    WR<-0.07203515/60 #lake average water respiration in mg/L/min
    # DO Saturation at a given temperature        
    DOsat<-4.09+(10.5*exp(-0.0371*tempC)) 
    SR<- ((((0.287*tempC-2.5)*44.661)*0.7)*0.001)/60
    # diffusive exchange
    DE<-(DO*DOsat)*(2.2*10^-5)*(0.03^-1)*
        60*10^-2 
    # below is eq 2. but i am not sure the do term should be ther
    #dDO<- DO - SR*Z^-1 + WR +DE*Z^-1
    # I think this is the right ode
    dDO<- -1*(SR*Z^-1+WR+DE*Z^-1)
    return(list(dDO))
    }
solution<- deSolve::ode(
    y=DO_dusk, 
    times=c(0:(10*60)), 
    func=DO_fun, 
    parms=parameters, 
    method="rk4")
plot(solution,ylab="DO mg/L",las=1,main="")

image