Explicitly include GPP, ER, Gas exchange in model outputs

[jsadler2]

One possible approach to adding process guidance to the vanilla deep learning model is to explicitly represent the DO mass balance by outputting GPP, ER, and Gas exchange:

Baseline LSTM (#40):

Proposed:

This is based on the Odum equation for DO mass balance:

[amcarter]

A couple more ideas for this version of the model:

1. We could add process guidance based on physical relationships between stream slope, width and depth with gas exchange (K)
2. Both GPP and ER could be constrained to have a degree of autocorrelation (unless there is high precipitation possibly)

[jsadler2]

1. We could add process guidance based on physical relationships between stream slope, width and depth with gas exchange (K)

This could be interesting. I guess there are ways to estimate K with just stream slope, width and depth? One question I have here, is that K actually doesn't appear in any of these equations as they are now. I guess we'd have to figure out how to do that. But then we'd also need a estimate of DO at that time step, right? And then my head starts to hurt because that's what we are trying to estimate in the first place, so it seems like it'd be on both sides of the equation.

[jsadler2]

2. Both GPP and ER could be constrained to have a degree of autocorrelation (unless there is high precipitation possibly)

I like this idea.

[lekoenig]

Both GPP and ER could be constrained to have a degree of autocorrelation (unless there is high precipitation possibly)

Thanks, Alice - I like this idea. I had wondered if we could impose a similar idea in the model even if we didn't include GPP and ER explicitly - i.e., could we assume that the daily DO amplitude does not change much from one day to the next (again, aside from a high precipitation event), and penalize the model accordingly.

[lekoenig]

^ Following up on my comment here, this approach of incorporating DO amplitude (and maybe the autocorrelation or change in DO amplitude from yesterday) in the model loss function seems like it might be more straightforward than predicting GPP, ER, and K. But I'm uncertain whether we try that before trying to incorporate the DO mass balance into the PGDL model (if at all).

[amcarter]

Also, an option for farther down the road would be to take this approach but then add on a second layer of neural network after the GPP, ER and K are predicted that included a subset of the original predictors. This would allow for non-deterministic relationships between metabolism and oxygen concentrations that arise because of variation in stream temperatures, differences in groundwater inputs, or other types of process error.

[jsadler2]

I think your other comment (#47) accounts for this by including h_t in the calculation of DO_min and DO_max

[amcarter]

One more thought - could we also estimate water temperature? Meaning that the estimated variables would be:

If so, then we could calculate DO saturation, and then have equations like the ones above to calculate DO max/mean/min:

I'm still don't fully understand the h, W, and b terms - so I apologize if these equations don't make sense!

[jsadler2]

Also if we are estimating water temp, we'd probably have to add that to the loss function:

Which also means we'd need water temperature values (either observations or estimates) for calculating the loss.

[amcarter]

I think it makes more sense to keep DO max and DO min separate, because the range doesn't preserve any information about the magnitude of DO, just the change each day.

[jsadler2]

Ah right. Good point.

[jsadler2]

@amcarter, @lekoenig

I was looking at these equations more and thinking about implementing them.

I'm wondering if for this one (and the DO_max) if it makes sense to have the h_t*W term? To me that term is like a fudge factor. It gives the model another degree of freedom. But I actually don't know if we need it. Is simplifying it to this accurate?

[jsadler2]

mmm. One thing that I just thought of is that DO_sat may not actually be the ceiling since streams can be super saturated. Is that right? So maybe we'd have to have some wiggle room either in this calculation or the calculation of DO_sat. Maybe you were thinking that that would already be the case.

[lekoenig]

Based on a suggestion from Bob, we could also edit the model equations at the top of the page such that we move from y_t to DO using DO-min rather than DO-mean. In other words, we predict DO-min as a function of ER, K, and maybe temperature, and then use DO-min combined with the weights and biases to predict DO-max and DO-mean in the subsequent equations.

One question I still have that maybe @jsadler2 can help me with - is the equation {DO-mean_t = GPP_t + ER_t + Gas_t} deterministic, or does that line represent relationships that the model is able to learn?

[jsadler2]

Based on a suggestion from Bob, we could also edit the model equations at the top of the page such that we move from y_t to DO using DO-min rather than DO-mean. In other words, we predict DO-min as a function of ER, K, and maybe temperature, and then use DO-min combined with the weights and biases to predict DO-max and DO-mean in the subsequent equations.

@lekoenig - I was thinking that the changes that @amcarter suggested (#47 (comment)) implements that idea

[jsadler2]

is the equation {DO-mean_t = GPP_t + ER_t + Gas_t} deterministic, or does that line represent relationships that the model is able to learn?

That would be deterministic. As in the model wouldn't learn anything about that equation. It would be learning about how to arrive at GPP, ERR and Gas, because those are dependent on the W_out and b_out which would be learned in training.

Does that answer your question?

[lekoenig]

Yep, that makes sense. Thanks!

[jsadler2]

I've started implementing this (#56) and run into a few questions that I'm hoping to get some feedback on:

Here are the equations that I'm working toward implementing (I made some adjustments that I thought made sense in the do_mean calculation):

And here is the code so far (from here):

        self.rnn_layer = layers.LSTM(
            hidden_size,
            return_sequences=True,
            recurrent_dropout=recurrent_dropout,
            dropout=dropout,
        )
        self.metab_out = layers.Dense(3)
        self.do_range_multiplier = layers.Dense(1)
        self.do_mean_wgt = layers.Dense(1)

    def call(self, inputs, DO_sat):
        h = self.rnn_layer(inputs)
        metab = self.metab_out(h)
        GPP = metab[:, :, 0]
        ER = metab[:, :, 1]
        K = metab[:, :, 2]
        DO_min = DO_sat - ER + K
        DO_max = DO_min + GPP - K
        DO_mean = DO_min + self.do_range_multiplier(DO_max - DO_min) + tf.squeeze(self.do_mean_wgt(h))
        return tf.stack((DO_min, DO_mean, DO_max, GPP, ER, K), axis=2)

My questions

1. Using DO_sat - I'm thinking that DO_sat will come directly from the data in the Appling data release. This is instead of being calculated by the stream temp as suggested above. I was thinking that using the DO_sat estimate from the data release would save us a step, but we don't have to go that route. Do you think that makes sense?

2. Gas exchage term - In the Odum equation, the gas exchange is represented as K * DO_deficit. I can't think of a simple way to know the DO_deficit at time t ... after all that is what we are kind of solving for (dO/dt is on the other side of the equation). I'm assuming this trickiness is what the MCMC is used for in the Appling model. Is that right?

   I don't think the code above is the right way to set the gas exchange term up. As is, I was thinking we'd train K based on the K in the data release, but using K like in the Odum equation assumes we know the DO deficit. What I think we maybe should do is more like it is in the equations - look at the entire gas exhange term as a whole and do not worry about K specifically. So something like this:

    def call(self, inputs, DO_sat):
        h = self.rnn_layer(inputs)
        metab = self.metab_out(h)
        GPP = metab[:, :, 0]
        ER = metab[:, :, 1]
        Gas= metab[:, :, 2]
        DO_min = DO_sat - ER + Gas
        DO_max = DO_min + GPP - Gas
        DO_mean = DO_min + self.do_range_multiplier(DO_max - DO_min) + tf.squeeze(self.do_mean_wgt(h))
        return tf.stack((DO_min, DO_mean, DO_max, GPP, ER, Gas), axis=2)

I have one question about this approach: if we take this approach, should we somehow incorporate K explicitly? We have estimates of K_600 in the Appling data release. Should we use those to somehow guide the model in estimating the gas exchange term? Or do we just let the model have "free reign" on that term?

[jsadler2]

It would be great to hear your thoughts on this, @lekoenig, @galengorski , @amcarter. Hopefully I explained it well enough. Also, I'm at the very fringe of my understanding of the bio/physical process, so please correct any of my misunderstanding

[amcarter]

Jeff, I'm not certain on these questions, but here are some thoughts:

1. Using DO_sat - I'm thinking that DO_sat will come directly from the data in the Appling data release.

That makes sense to me, except I am wondering where it will come from in the streams where we are trying to make DO predictions and do not have metabolism data (and therefore DO_sat) available as an input. We could calculate it from stream temperature and air pressure, but I thought our plan was to rely on the model to estimate stream temperature rather than limiting prediction to basins where we have that data.

2. Gas exchange term

I agree that these equations don't seem quite right. It would be circular to try and calculate the gas exchange based on oxygen concentrations at time t - as you say farther down it makes more sense to estimate the actual flux of gas (in mg O2/m2/d) at each time point. If we go that route, I think we should calculate this at each time point using the K600 and O2 concentrations in the data release rather than giving it free reign. I still unsure about how we are incorporating that term though, even in the second set of equations. I'll need to spend a little more time thinking that through to convince myself - then I'll get back to you. Apologies for missing this for a few days.

[galengorski]

I too am at the fringe of my understanding here, but a few questions come to mind:

1. Using DO_sat - I'm thinking that DO_sat will come directly from the data in the Appling data release.

I agree with the point made by @amcarter about how this approach would present difficulties in extending the analysis to areas where we don't have metabolism estimates. However if the goals of the study are to investigate if PGDL improves performance, that doesn't necessarily have to be done at every location. That being said, if we could estimate DO_sat using stream temperature and air pressure that would increase the extensibility of the method. Is there a possibility of predicting DO_sat? Is there an advantage to directly using estimates of DO_sat from the data release rather than estimating them explicitly using the model?

2. Gas exchange term

I'm probably missing something here but isn't K already constrained in the second set of equations above? If you're predicting the gas exchange using the gas exchange from the data release, I thought the gas exchange from the data release uses K600 (fit using their model) and measured water temperature. Equations 4 and 5 from Appling et al.:

[amcarter]

I spent some time working through this this afternoon - the original equations we were using were much more of a back of the envelope ballpark of O max and min. I think that it is still doable to estimate these from GPP and ER, but we will need to make some simplifying assumptions. Below is a walkthrough my current thinking of how we would need to approach this - I worked through these equations with Bob this afternoon and he is also tentatively on board. I will try to test this on some data and maybe have some plots to share tomorrow

[amcarter]

[amcarter]

Following up on our conversation this morning, these are the updated equations for DO min and max

here, both Osat and K2 are calculated as a function of a temperature T - I've tried a few different ones below. Note that this formulation does require an estimate of depth (z) unless we want to train the model on k2, which we could calculate using the depth in the data release.

First, here are some comparisons for how these equations work for predicting DO max and DO min using:

1. Daily mean temperature (meanT)
2. Daily max and daily min temperature (deltaT)
3. Temperature at the time of maximum light (Lmax, for DO max only)

The good news is that it doesn't look like the differences in temperature choice are very important, and using the daily mean water temperature (or possibly even air temperature) would suffice.

Second, in the equation for DO max above I've decided we should keep the GPP max calculation, but that means we need to add a timestep correction factor on ER and K2, which is different from what we talked about today in the meeting. That does mean we'd have to model within day light based on latitude and day of year, and that the estimates would change a little depending on the timestep we use to calculate it (but probably not by much). After fixing the ER timestep issue, this correction seems to work the best of all the options. Below I compare it to correcting to the average daily GPP based on daylength and to an uncorrected estimates for all three temperature methods above.

[jsadler2]

Thanks for this update, @amcarter. I think I'm understanding this, but I want to make sure before implementing in code.

estimating k2

I think it makes sense to use k2 because we can calculate this with k600, z, and T from the metabolism data release. My thinking is that

• the model will estimate these three values
• then we'll use k600 and the estimated water temp with this equation to get K2
• 
• then to get k2 we do K2 x z

estimating light, L

For the light term, L_t/\sum{L}, you said we can estimate that using the day of year and latitude. And then I assume that L_t is the light at the maximum time of day. Is that right? If so, can you point me to the equations for estimating this? Also, just to check, this ratio is a ratio of the max light to the total amount of light?

no GPP in calculating DO min

I'm curious about the equation for calculating DO_min. It makes sense mathematically how you derived that, but (bio)logically, I'm curious why GPP doesn't factor in. If there are two stream locations with the same ER rate, and one has a higher GPP rate, shouldn't the one with the higher GPP rate have a higher DO_min, since its DO_max will be higher? Not trying to poke holes, just wanting to understand.

time step correction factor

Is the "timestep correction factor" \delta t? And are we including that because we shouldn't be counting ER for the entire 24 hours when GPP is only happening in the daylight hours (so \delta t would be less than 1 since we are working on a daily time scale)?

[amcarter]

@jsadler2, here are my thoughts on these - and I'm happy to talk more about them if you'd like.

estimating k2

Yes, this looks exactly right.

estimating light, L
I assume that L_t is the light at the maximum time of day

Yes, this is the daily max light. L_t/\sum{L} is the max light value for each day divided by the cumulative light for that day. These can be estimated using functions from streamMetabolizer. I'm not actually sure what these functions are referencing under the hood, but basically they are adjusting the times so that noon matches solar noon, then calculating a daily curve light curve based on time of year. Here, datetime is a vector of 'yyyy-mm-dd HH:MM:SS' in the local time zone.

library(streamMetabolizer)  
solar.time = calc_solar_time(datetime, longitude)  
light = calc_light(solar.time, latitude, longitude)

no GPP in calculating DO min
I'm curious why GPP doesn't factor in

Yeah, it was kind of hard for me to wrap my head around this too - but the GPP is just setting the amplitude of the DO curve each day, which is why it's much easier to estimate accurately even if the DO sensor is miscalibrated. The ER determines the location of that curve relative to DO saturation. So two sites with the same ER but different GPP would in theory have the same minimum but different maximums, assuming temperature and depth are all the same.

time step correction factor
Is the "timestep correction factor" \delta t?

I don't think that timestep correction factor is actually a great word to describe this, so sorry for making it confusing. It's just that when we are solving the odum equation for DO, the dt gets moved to the other side of the equation and is multiplied by each of the fluxes (GPP, ER, D). The thing I was getting wrong in the first place was thinking I could just cancel it out because it was part of every term. The problem is that the L_max/sum(L) term is actually accounting for the timestep in GPP, meaning that it can't be canceled from the ER and K terms. This difference in how we deal with the timesteps in each term is because our estimates for GPP, ER, and D are single values calculated for each day, but the assumption in the model is that they are distributed differently throughout the day. ER is assumed to be constant over the 24 hrs, so at a single timestep, the amount of ER happening is ER/(number of timesteps per day), and similarly with diffusion (that's why the delta T does cancel in the DO min equation). GPP however is distributed unevenly throughout the day and is assumed to be proportional to the light availability. That's why it is zero at night and during the timestep we're interested in it is the fraction of total light that occurs in that timestep.
tl;dr - delta T should be equal to 1/(number of timesteps per day) or 1/48 for 30 minute resolution data because the L_max/sum(L) is accounting for this in the GPP.

[jsadler2]

Thanks for those explanations, @amcarter. That's very helpful.

I'm still wondering about how to estimate both parts of the light term.

From what I gather streamMetabolizer estimates instantaneous light at a given time. But I'm not clear on how we'd get L_t from that since, if I'm understanding the \delta t thing right, we need light over the time period t. Equally puzzling to me is how we'd get the cumulative light, \sum(L)

I think this will be a good topic for tomorrow's conversation.