diff --git a/docs/whats-new.rst b/docs/whats-new.rst
index 56f7f0a25..358d83229 100644
--- a/docs/whats-new.rst
+++ b/docs/whats-new.rst
@@ -8,7 +8,8 @@ v1.6.x (not released)
 
 Enhancements
 ~~~~~~~~~~~~
-
+- Resolved inverted sign of flowline diagnostics flux divergence (:pull:`1815`).
+  By `Brandon Tober <https://github.com/btobers>`_
 - Added ``tasks.compute_fl_diagnostics_quantiles``, this task is designed to
   calculate quantiles from multiple fl_diagnostic files. It enables users to
   compute metrics such as the median flowline across various GCM projections
diff --git a/oggm/core/flowline.py b/oggm/core/flowline.py
index b66e8b36d..91785eec2 100644
--- a/oggm/core/flowline.py
+++ b/oggm/core/flowline.py
@@ -1346,8 +1346,9 @@ def run_until_and_store(self, y1,
                                 val * fac_sec)
                             surface_h_previous[fl_id] = fl.surface_h
                         if 'flux_divergence' in ovars_fl and (yr > self.y0):
-                            # calculated after the formula dhdt = mb + flux_div
-                            val = ds['dhdt'].data[i, :] - ds['climatic_mb'].data[i, :]
+                            # calculated after the formula mb = dhdt + flux_div
+                            # Cogley et al. (2010): Glossary of Glacier Mass Balance
+                            val = ds['climatic_mb'].data[i, :] - ds['dhdt'].data[i, :]
                             # special treatment for retreating: If the glacier
                             # terminus is getting ice free it means the
                             # climatic mass balance is more negative than the
diff --git a/oggm/tests/test_models.py b/oggm/tests/test_models.py
index d8557a89f..30ac39650 100644
--- a/oggm/tests/test_models.py
+++ b/oggm/tests/test_models.py
@@ -1823,7 +1823,7 @@ def test_run(self, class_case_dir):
                     np.isclose(model.fls[0].thick, 0.),
                     dhdt_ref[-1] < 0.)
                 flux_divergence_ref.append(
-                    (dhdt_ref[-1] - climatic_mb_ref[-1]) *
+                    (climatic_mb_ref[-1] - dhdt_ref[-1]) *
                     np.where(has_become_ice_free, 0.1, 1.))
             if int(yr) == 500:
                 secfortest = model.fls[0].section