Level Pool Mass Balance Bug

[hellkite500]

The level pool formulation doesn't correctly conserve mass when exceeding max conditionions.

Expected Behavior

When a reservoir is below max, and the current conditions would cause it to exceed max, it should fill the available storage and release the excess water.

Current Behavior

When current conditions cause the reservoir to exceed max, storage becomes max and all inflow is passed to outflow. This creates water in the reservoir!

Possible Solution

292 PR

When the reservoir exceeds max, the discharge should be discharge = qi1 - (((maxh - Htmp) * sap) / dt)

When the reservoir was not full, but would exceed, maxh-Htmp = available storage and only the excess inflow goes to outflow.

When the reservoir is already full, maxh-htmp = 0, and all inflow goes to outflow.

[hellkite500]

As has been realized, this fix DOES NOT address what happens to orifice + weir outflows in the overtop condition. The proposed solution doesn't account for these, but could be modified to do so.

discharge = qi1 - (((maxh - Htmp) * sap) / dt) + (orifice + weir)
H = <recompute based on Htmp, discharge, qi1>

Also need to appropriately recompute H (it may not longer be maxh since the extra volume is removed from the reservoir.) This would likely mean that we maintain NEARLY max elevation, minus the extent of the weir to move flow, so by the next time step we are below max again, then will cross over with the new inflow, move the excess and however much the weir can handle, and the elevation is once again below max.

[jmccreight]

@jdmattern-noaa attached previously attached this information to #292 which add further detail:

Modification on the Level Pool Method

This is an example of how the new fix would impact the outflow simulated by the Level Pool (LP) method. We added here results for the lake_id ==16944276. Note that graphs may not clearly represent the difference. For this reservoir, the water elevations are equal before and after the LP fix. The outflows are equal at every timestep except at 12300 seconds. This is the point where LP calculates the water elevation to exceed its maximum value. The fixed LP calculates 218.43 cms less outflow than the unfixed LP, which is the amount of water needed to fill the reservoir to the max elevation. Therefore, 512.60 cms should be the outflow because it is the excess water that spills over the top at that timestep as opposed to releasing the full inflow of 731.03 cms.

[hellkite500]

As has been realized, this fix DOES NOT address what happens to orifice + weir outflows in the overtop condition. The proposed solution doesn't account for these, but could be modified to do so.

discharge = qi1 - (((maxh - Htmp) * sap) / dt) + (orifice + weir)
H = <recompute based on Htmp, discharge, qi1>

Also need to appropriately recompute H (it may not longer be maxh since the extra volume is removed from the reservoir.) This would likely mean that we maintain NEARLY max elevation, minus the extent of the weir to move flow, so by the next time step we are below max again, then will cross over with the new inflow, move the excess and however much the weir can handle, and the elevation is once again below max.

@laurareads @dnyates Reference to latest Level Pool formulation discussion.

[laurareads]

Stay tuned for modifications according to @nels comment. We will pass unit tests on a small domain and then send off to you for comments.