It seems that the flux divergence has an inverted sign than it should.
I believe that line 1350 of oggm.core.flowline.py should read:
ds['climatic_mb'].data[i, :] - ds['dhdt'].data[i, :]
instead of:
ds['dhdt'].data[i, :] - ds['climatic_mb'].data[i, :]
per Cogley et al. (2010)'s Glossary of Glacier Mass Balance and Related Terms.