Bug in GOTM water balance?

[shajar007]

When using GOTM with lake support version 5.4.0 , it is possible to set parameters so that GOTM will calculate lake level. The parameters are water_balance_method = 3 and zeta_method=3.
I guess these parameters are rarely used because many lakes have constant lake level. However, at least in lake Kinneret, calculating lake level is important. By calculating daily water balance out of the GOTM output file it can be shown that GOTM fails to calculate lake level correctly when using these parameters. If lake level is supplied to the model externally using parameters: water_balance_method = 2, zeta_method=2, GOTM calculates the water balance correctly using a residual stream.
Water balance calculation:
Daily water balance = Inflow + (evaporation + precipitation) * surface_area + Residual_stream – Δ_lake_level * surface_area.
Where the variables are from the GOTM output.nc file:
Inflow=Qlayer [m3/s] * 86400
Evaporation = evap [m/s] * 86400
Precipitation = percip [m/s] * 86400
surface_area = Af [m2]
Residual_stream = Qres [m3/s] * 86400
Δ_lake_level = number of layers * (h(t) - h(t-1)) [m/day].

(86400 is converting seconds into a day)
Daily water balance should be zero.
When calculating the water balance with parameters (water_balance_method = 2, zeta_method=2) the daily water balance is exactly zero when there is no change in lake level, and around zero when there is lake level change. This suggests that changes in lake volume is calculated different than Δ_lake_level * surface_area – maybe it is calculated hourly where I use the output file which is daily – resulting in a small difference. The daily balance is +/- 400 m3/day in my lake Kinneret model.
When calculating the water balance with parameters (water_balance_method = 3, zeta_method=3) the daily water balance is off by +/- 2*10^6 m3/day ! The water balance is obviously off with these parameters resulting in the inability to calculate lake level.
This issue can probably be reproduced with any lake model with these water balance parameters. Attached is an R code to calculate water balance, assuming output is daily.
Am I missing a parameter? Thanks

GOTM_water_balance.zip

[knutaros]

Hi Xia, thanks for reporting. Please see some first comments below.

On 1/12/22 08:50, shajar007 wrote: When using GOTM with lake support version 5.4.0 , it is possible to set parameters so that GOTM will calculate lake level. The parameters are water_balance_method = 3 and zeta_method=3. I guess these parameters are rarely used because many lakes have constant lake level. However, at least in lake Kinneret, calculating lake level is important. By calculating daily water balance out of the GOTM output file it can be shown that GOTM fails to calculate lake level correctly when using these parameters. If lake level is supplied to the model externally using parameters: water_balance_method = 2, zeta_method=2, GOTM calculates the water balance correctly using a residual stream. Water balance calculation: Daily water balance = Inflow + (evaporation + precipitation) * surface_area + Residual_stream – Δ_lake_level * surface_area. Where the variables are from the GOTM output.nc file: Inflow=Qlayer [m3/s] * 86400 Evaporation = evap [m/s] * 86400 Precipitation = percip [m/s] * 86400 surface_area = Af [m2] Residual_stream = Qres [m3/s] * 86400 Δ_lake_level = number of layers * (h(t) - h(t-1)) [m/day]. (86400 is converting seconds into a day) Daily water balance should be zero. When calculating the water balance with parameters (water_balance_method = 2, zeta_method=2) the daily water balance is exactly zero when there is no change in lake level, and around zero when there is lake level change. This suggests that changes in lake volume is calculated different than Δ_lake_level * surface_area – maybe it is calculated hourly where I use the output file which is daily – resulting in a small difference. The daily balance is +/- 400 m3/day in my lake Kinneret model.

Your formula above does not consider the time-dependent surface_area. Maybe you want to try to calculate the time-dependent lake volume as the sum of the layer volumes. To do so please try registering Vc for output: +++ b/src/gotm/register_all_variables.F90 @@ -354,6 +354,7 @@     call fm%register('ga', '', 'coordinate scaling', standard_name='??', dimensions=(/id_dim_z/), data1d=ga(1:nlev),category='column_structure')     if (lake) then        call fm%register('Af', 'm^2', 'hypsograph at grid interfaces', standard_name='??', dimensions=(/id_dim_z/), data1d=Af(1:nlev), category='column_structure') +      call fm%register('Vc', 'm^3', 'volume of layers', standard_name='??', dimensions=(/id_dim_z/), data1d=Vc(1:nlev), category='column_structure')        call fm%register('int_flow', 'm3', 'stream: integrated', data0d=int_flows, category='meanflow/waterbalance')        call fm%register('int_water_balance', 'm3', 'integrated total water balance', data0d=int_water_balance, category='meanflow/waterbalance')     else

When calculating the water balance with parameters (water_balance_method = 3, zeta_method=3) the daily water balance is off by +/- 2*10^6 m3/day ! The water balance is obviously off with these parameters resulting in the inability to calculate lake level.

based on your inflow and net_precip data between two instants of time, can you diagnose whether the change of model sea level is completely wrong (sign or orders of magnitude)? Apart from that: is the model sea level always below your uppermost input hypsography level? Cheers, Knut

…

This issue can probably be reproduced with any lake model with these water balance parameters. Attached is an R code to calculate water balance, assuming output is daily. Am I missing a parameter? Thanks GOTM_water_balance.zip <https://github.com/gotm-model/code/files/7852543/GOTM_water_balance.zip> — Reply to this email directly, view it on GitHub <#30>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AC2RV6UFANW7IUH36F7DAZLUVUXFXANCNFSM5LYIS7QQ>. Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[shajar007]

Hi Knut,
I know my calculation is simplified, but for a lake area of 165,000,000 m^2, lake level makes a small difference, much smaller than the differences between the 2 water balance methods.
Regarding registering Vc for output - I understand that I have to change something in the Fortran code and compile it? I don't know how to do it - I will need guidance.
The change between 2 consecutive days is not completely wrong in terms of lake level - the mistake is about 1 mm per day, but it amounts to 1.65 meter in 4 years. The daily volume difference have positive and negative values, tending for the negative side which amounts to a lake level decrease. (so the lake level is always below hypsography).
The mean absolute daily bias in water balance is 493398 m^3, but the mean bias is 286730 m^3 - i.e. the bias is both negative and positive, more to the positive side. Just to compare, when using external lake level and residual stream, the mean absolute bias is 165 M^3 and the mean bias is -3 m^3.
I attach my calculation in a csv file for a 4 year run.
water_balance1.csv
Thanks! Shajar

[shajar007]

Update
I run the model with hourly and daily output – again for both water balance methods.
Calculating water balance from the hourly output showed good results (up to +/- 400 m3 mismatch in water balance for both methods), however lake level had the same error as before.
The water balance in the hourly output did not add-up to the daily output! Most notably the evaporation in the hourly output (after aggregating to daily) is about 1.6 times the evaporation in the daily output. The residual stream compensates the difference in the evaporation.
My guess is that evaporation in the hourly output is correct. But aggregation to a day(?) for the evaporation has a bug. In water_balance_method=2 the residual stream compensates for the evaporation error where in the water_balance_method=3 the error is apparent in the water balance calculation described above.
I hope this help pinpointing the bug.
Regards, Shajar

[TobiasKAndersen]

Hi Shajar,
You are running GOTM coupled to WET, right? From my understanding of the GOTM output, if you want to get daily values you should change output to daily and the time step for integration (dt). If you have set time step to hourly and output to daily, you will get daily output at a specific time step (and not integrated for the entire day). So I think you should stick with the hourly calculations.
I hope this clears up some misunderstanding.
Cheers, Tobias

[shajar007]

Hmmm. I'll reinspect and update accordingly.
Thanks! Shajar

[shajar007]

I rechecked my calculations according to Tobias comment, I think my calculations are correct.
The rates are given per second. In case of hourly output I multiply by 3600 to get hourly rate, multiply by surface area for evaporation and rain to get water volume in m^3.
In case of daily output I multiply by 86400 to get daily rate and do the same as above to get water volume.
Aggregation of hourly volume into days should approximately equal the daily volume. The evaporation and residual stream are significantly different between daily and hourly output.
Of course this is an approximated calculation because it assumes rate is constant for the time step of the output (day or hour) and lake surface is taken from last time step. However these should not make big differences.
Shajar

[shajar007]

Hi @knutaros , @TobiasKAndersen
I further investigated the problem and found that probably the bug is not in the water balance but rather in the daily output of evaporation and residual stream.

1. when parameters are water_balance_method = 3 and zeta_method=3, I run once with daily output and once with hourly output and compare. I found that the daily evaporation is taken from evaporation at 12:00:00 multiplied by 24. This makes a big difference since evaporation at 12:00 noon is not representative of the whole day! So practically - evaporation in the daily output is useless.
2. when parameters are water_balance_method = 2 and zeta_method=2 , I check the same as above. In addition to the evaporation issue, the residual stream also has great difference between hourly and daily output. Here the hourly residual stream, after aggregation to daily, is about 1000 times the daily. obviously this is wrong, because the daily output should be aggregation of an hourly data.
3. other variables - like inflows, also doesn't aggregate the hourly data to be exactly as the daily data, but differences are small.
   Hope you find the time to look into it.
   Regards, Shajar