Describe the bug
I've found a bug with the way WRF Hydro handles level pool physics for lakes/reservoirs. Fortran by default uses 4 byte floats, which only reliably stores around 7 decimal points. This is in contrast with "double precision" 8 byte floats which have at least 15 decimal points of precision. From the level pool code in WRF Hydro, it appears the entire model is using 4 byte floats: https://github.com/NCAR/wrf_hydro_nwm_public/blob/main/src/Routing/Reservoirs/Level_Pool/module_levelpool.F90
The is problematic when the lake surface area is large and the routing timestep is small. The specific problem is calculating the change in height for the lake. From the source code:
dh1 = ((It - discharge)/sap)*dt
dt = routing time step
sap = surface area of the lake in m^2
If the routing timestep dt is small (let's say <10 seconds) and the surface area sap is large (let's say 150km^2 = 150,000,000m^2), you can quite evidently see that dh1 will push towards the edge of the floating point precision, and then the lake height is calculated (H = H + dh), it will cause issues where H would not change even if dh > 0.
An example of this is in a Google forum thread I posted a few days ago, where elevation for the most would not change unless the inflow was so high that it essentially forced dh to be significantly higher than the floating point imprecision. It would otherwise stay flat even when we expected a delta in lake elevation.
Related to https://groups.google.com/a/ucar.edu/g/wrf-hydro_users/c/aYQuEEzdhYA
I'm not sure what solution is best to fix this. You could switch everything over for the lakes to use 8 byte floating point numbers instead of 4 byte floating point numbers to increase the precision. The downside of that is it would be different than the rest of WRF Hydro, and it would take more RAM and compute probably. Modifying the algorithm to work with millimeters instead of meters could also work, but might be messy and still might be problematic with even larger lakes. Please advise and let me know when a fix can be implemented for this.
To Reproduce
Steps to reproduce the behavior:
- Enable lake outputs and look specifically at large lakes with surface area >100km^2
- Set the routing timesteps DTRT_CH and DTRT_TER to a low value like 2. Note that this is a realistic setting when the grid cell sizes are ~30meters: https://wrf-hydro.readthedocs.io/en/latest/model-physics.html#id7
Expected behavior
Water surface elevation should constantly change and not remain constant
Screenshots
See https://groups.google.com/a/ucar.edu/g/wrf-hydro_users/c/aYQuEEzdhYA
Additional context
N/A
Describe the bug
I've found a bug with the way WRF Hydro handles level pool physics for lakes/reservoirs. Fortran by default uses 4 byte floats, which only reliably stores around 7 decimal points. This is in contrast with "double precision" 8 byte floats which have at least 15 decimal points of precision. From the level pool code in WRF Hydro, it appears the entire model is using 4 byte floats: https://github.com/NCAR/wrf_hydro_nwm_public/blob/main/src/Routing/Reservoirs/Level_Pool/module_levelpool.F90
The is problematic when the lake surface area is large and the routing timestep is small. The specific problem is calculating the change in height for the lake. From the source code:
dh1 = ((It - discharge)/sap)*dt
dt = routing time step
sap = surface area of the lake in m^2
If the routing timestep dt is small (let's say <10 seconds) and the surface area sap is large (let's say 150km^2 = 150,000,000m^2), you can quite evidently see that dh1 will push towards the edge of the floating point precision, and then the lake height is calculated (H = H + dh), it will cause issues where H would not change even if dh > 0.
An example of this is in a Google forum thread I posted a few days ago, where elevation for the most would not change unless the inflow was so high that it essentially forced dh to be significantly higher than the floating point imprecision. It would otherwise stay flat even when we expected a delta in lake elevation.
Related to https://groups.google.com/a/ucar.edu/g/wrf-hydro_users/c/aYQuEEzdhYA
I'm not sure what solution is best to fix this. You could switch everything over for the lakes to use 8 byte floating point numbers instead of 4 byte floating point numbers to increase the precision. The downside of that is it would be different than the rest of WRF Hydro, and it would take more RAM and compute probably. Modifying the algorithm to work with millimeters instead of meters could also work, but might be messy and still might be problematic with even larger lakes. Please advise and let me know when a fix can be implemented for this.
To Reproduce
Steps to reproduce the behavior:
Expected behavior
Water surface elevation should constantly change and not remain constant
Screenshots
See https://groups.google.com/a/ucar.edu/g/wrf-hydro_users/c/aYQuEEzdhYA
Additional context
N/A