Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Thermocline depth calculation #118

Open
vikifc opened this issue Oct 26, 2023 · 4 comments
Open

Thermocline depth calculation #118

vikifc opened this issue Oct 26, 2023 · 4 comments

Comments

@vikifc
Copy link

vikifc commented Oct 26, 2023

HI! I am finding a problem when calculating the thermocline depth of a temperature profile in R. Here it goes:

wtr = c(13.9, 13.8, 13.7, 13.7, 13.5, 13.3, 13.1, 13.0, 12.9, 12.8)
water.density(wtr)
[1] 999.2871 999.3007 999.3143 999.3143 999.3411 999.3674 999.3931 999.4059 999.4184 999.4309
depths = c(0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5)
thermo.depth(wtr, depths, Smin = 0.5, seasonal = TRUE, mixed.cutoff = 1)
[1] 1.991029

I don't understand why it finds a thermocline at 1.99 m with the Smin=0.5 condition. I also tried with really high and low Smin numbers and the result is the same, do someone know how to fix that?

I would really appreciate the help, thanks!

@hdugan
Copy link

hdugan commented Oct 27, 2023

Smin is a bit confusing, because it is easily overridden in the code if it's set too high or too low.

What are you trying to accomplish? Would the mixed.cutoff parameter be better suited to your goal?

@vikifc
Copy link
Author

vikifc commented Oct 27, 2023

I think mixed.cutoff considers the maximum difference between de highest temperature and the lowest temperature in the whole profile, right?
I am trying to calculate the depth of the thermocline, where both conditions are met: 1) the density gradient is at least 0.5 kg/m3/m, and 2) the temperature gradient is at least 1ºC/m. At 1.99 m none of that is true.

@hdugan
Copy link

hdugan commented Oct 29, 2023

Yeah I seem the problem. Without recoding the current function, you could apply your own density cutoff. Something like the following:

library(LakeMetabolizer)
library(tidyverse)
wtr = c(13.9, 13.8, 13.7, 13.7, 13.5, 13.3, 13.1, 13.0, 12.9, 12.8)
depths = c(0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5)

# Derivative function
drv <- function(x, y) c(NA, (diff(y) /diff(x))) 

# Create your own density cutoff 
data.frame(wtr = wtr, depths = depths) |> 
  mutate(dens = water.density(wtr)) |> 
  mutate(dens.der = drv(depths, dens)) |> 
  mutate(thermo.depth = thermo.depth(wtr,depths)) |> 
  mutate(thermo.depth = if_else(max(dens.der, na.rm = T) < 0.5, NA, thermo.depth[1]))

Using the data.frame structure probably makes more sense if you have multiple profiles and group them accordingly by date. But hopefully this is useful.

@vikifc
Copy link
Author

vikifc commented Oct 30, 2023

Yes, I think I'll go for a solution like that. Thank you for your help!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants