diff --git a/opm/core/fluid/IncompPropertiesBasic.cpp b/opm/core/fluid/IncompPropertiesBasic.cpp
index 1b93001b..0ce3c88a 100644
--- a/opm/core/fluid/IncompPropertiesBasic.cpp
+++ b/opm/core/fluid/IncompPropertiesBasic.cpp
@@ -115,6 +115,16 @@ namespace Opm
     /// \return Array of P density values.
     const double* IncompPropertiesBasic::density() const
     {
+        // No difference between reservoir and surface densities
+        // modelled by this class.
+	return pvt_.surfaceDensities();
+    }
+
+    /// \return Array of P density values.
+    const double* IncompPropertiesBasic::surfaceDensity() const
+    {
+        // No difference between reservoir and surface densities
+        // modelled by this class.
 	return pvt_.surfaceDensities();
     }
 
diff --git a/opm/core/fluid/IncompPropertiesBasic.hpp b/opm/core/fluid/IncompPropertiesBasic.hpp
index 9784398f..3c8535df 100644
--- a/opm/core/fluid/IncompPropertiesBasic.hpp
+++ b/opm/core/fluid/IncompPropertiesBasic.hpp
@@ -92,9 +92,14 @@ namespace Opm
         /// \return Array of P viscosity values.
         virtual const double* viscosity() const;
 
+        /// Densities of fluid phases at reservoir conditions.
         /// \return Array of P density values.
         virtual const double* density() const;
 
+        /// Densities of fluid phases at surface conditions.
+        /// \return Array of P density values.
+        virtual const double* surfaceDensity() const;
+
         /// \param[in]  n      Number of data points.
         /// \param[in]  s      Array of nP saturation values.
         /// \param[in]  cells  Array of n cell indices to be associated with the s values.
diff --git a/opm/core/fluid/IncompPropertiesFromDeck.cpp b/opm/core/fluid/IncompPropertiesFromDeck.cpp
index 2789fe92..e61e13c8 100644
--- a/opm/core/fluid/IncompPropertiesFromDeck.cpp
+++ b/opm/core/fluid/IncompPropertiesFromDeck.cpp
@@ -90,6 +90,12 @@ namespace Opm
 	return pvt_.reservoirDensities();
     }
 
+    /// \return Array of P density values.
+    const double* IncompPropertiesFromDeck::surfaceDensity() const
+    {
+	return pvt_.surfaceDensities();
+    }
+
     /// \param[in]  n      Number of data points.
     /// \param[in]  s      Array of nP saturation values.
     /// \param[in]  cells  Array of n cell indices to be associated with the s values.
diff --git a/opm/core/fluid/IncompPropertiesFromDeck.hpp b/opm/core/fluid/IncompPropertiesFromDeck.hpp
index 1472f043..d17dd1b7 100644
--- a/opm/core/fluid/IncompPropertiesFromDeck.hpp
+++ b/opm/core/fluid/IncompPropertiesFromDeck.hpp
@@ -78,9 +78,14 @@ namespace Opm
         /// \return Array of P viscosity values.
         virtual const double* viscosity() const;
 
+        /// Densities of fluid phases at reservoir conditions.
         /// \return Array of P density values.
         virtual const double* density() const;
 
+        /// Densities of fluid phases at surface conditions.
+        /// \return Array of P density values.
+        virtual const double* surfaceDensity() const;
+
         /// \param[in]  n      Number of data points.
         /// \param[in]  s      Array of nP saturation values.
         /// \param[in]  cells  Array of n cell indices to be associated with the s values.
diff --git a/opm/core/fluid/IncompPropertiesInterface.hpp b/opm/core/fluid/IncompPropertiesInterface.hpp
index f3c77caa..a5025e8a 100644
--- a/opm/core/fluid/IncompPropertiesInterface.hpp
+++ b/opm/core/fluid/IncompPropertiesInterface.hpp
@@ -62,9 +62,22 @@ namespace Opm
         /// \return Array of P viscosity values.
         virtual const double* viscosity() const = 0;
 
+        /// Densities of fluid phases at surface conditions.
         /// \return Array of P density values.
         virtual const double* density() const = 0;
 
+        /// Densities of fluid phases at surface conditions.
+        /// Note: a reasonable question to ask is why there can be
+        /// different densities at surface and reservoir conditions,
+        /// when the phases are assumed incompressible. The answer is
+        /// that even if we approximate the phases as being
+        /// incompressible during simulation, the density difference
+        /// between surface and reservoir may be larger. For accurate
+        /// reporting and using data given in terms of surface values,
+        /// we need to handle this difference.
+        /// \return Array of P density values.
+        virtual const double* surfaceDensity() const = 0;
+
         /// \param[in]  n      Number of data points.
         /// \param[in]  s      Array of nP saturation values.
         /// \param[in]  cells  Array of n cell indices to be associated with the s values.