diff --git a/landice/output_processing_li/plot_globalStats.py b/landice/output_processing_li/plot_globalStats.py
index 617004818..9800efdf5 100755
--- a/landice/output_processing_li/plot_globalStats.py
+++ b/landice/output_processing_li/plot_globalStats.py
@@ -14,7 +14,10 @@
 import textwrap
 
 rhoi = 910.0
-rhosw = 1028.
+rhosw = 1028.0
+rhofw = 1000.0
+Aocn = 3.625e14  # m^2 (Gregory et al. 2019)
+m_to_mm = 1000.0
 
 print("** Gathering information.  (Invoke with --help for more details. All arguments are optional)")
 parser = ArgumentParser(description=__doc__)
@@ -188,10 +191,16 @@
 
 
 def VAF2seaLevel(vol):
-    return vol / scaleVol / 3.62e14 * rhoi / rhosw * 1000.
+    """
+    This function accounts for the fact that ice, when melted,
+    takes on freshwater density.  This density correction is only
+    applied to the volume above flotation, when in reality, it should
+    be applied for all melted ice.
+    """
+    return vol / scaleVol / Aocn * rhoi / rhofw * m_to_mm
 
 def seaLevel2VAF(vol):
-    return vol * scaleVol * 3.62e14 * rhosw / rhoi / 1000.
+    return vol * scaleVol * Aocn * rhofw / rhoi / m_to_mm
 
 def addSeaLevAx(axName):
     seaLevAx = axName.secondary_yaxis('right', functions=(VAF2seaLevel, seaLevel2VAF))