diff --git a/dumux/material/binarycoefficients/h2o_heavyoil.hh b/dumux/material/binarycoefficients/h2o_heavyoil.hh
index cead97c2f5147da6159fed9bb45bf574bd07a896..1dae2afa265f21a969f1b65f45c1552c9c1e8bb1 100644
--- a/dumux/material/binarycoefficients/h2o_heavyoil.hh
+++ b/dumux/material/binarycoefficients/h2o_heavyoil.hh
@@ -41,10 +41,8 @@ public:
     /*!
      * \brief Henry coefficient \f$[N/m^2]\f$  for heavy oil in liquid water.
      *
-     * See:
-     *
+     * \todo values copied from TCE, please improve it
      */
-
     template <class Scalar>
     static Scalar henryOilInWater(Scalar temperature)
     {
@@ -57,10 +55,8 @@ public:
     /*!
      * \brief Henry coefficient \f$[N/m^2]\f$  for water in liquid heavy oil.
      *
-     * See:
-     *
+     * \todo arbitrary value, please improve it
      */
-
     template <class Scalar>
     static Scalar henryWaterInOil(Scalar temperature)
     {
@@ -72,17 +68,18 @@ public:
     /*!
      * \brief Binary diffusion coefficient [m^2/s] for molecular water and heavy oil.
      *
+     * \todo value is just an order of magnitude, please improve it
      */
     template <class Scalar>
     static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
     {
-        return 1e-6; // [m^2/s] This is just an order of magnitude. Please improve it!
+        return 1e-6; // [m^2/s] TODO: This is just an order of magnitude. Please improve it!
     }
 
     /*!
      * \brief Diffusion coefficient [m^2/s] for tce in liquid water.
      *
-     * \todo
+     * \todo value is just an order of magnitude, please improve it
      */
     template <class Scalar>
     static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
diff --git a/dumux/material/fluidsystems/1p.hh b/dumux/material/fluidsystems/1p.hh
index 9131015e2e5ea8a282d1453358f97fb49667aaac..6d664b1fcd6dddb43d3a372a6bc57068c36fdbd0 100644
--- a/dumux/material/fluidsystems/1p.hh
+++ b/dumux/material/fluidsystems/1p.hh
@@ -265,11 +265,13 @@ public:
      * \brief Calculate the fugacity coefficient \f$\mathrm{[Pa]}\f$ of an individual
      *        component in a fluid phase
      *
-     * The fugacity coefficient \f$\mathrm{\phi_\kappa}\f$ is connected to the
-     * fugacity \f$\mathrm{f_\kappa}\f$ and the component's molarity
-     * \f$\mathrm{x_\kappa}\f$ by means of the relation
+     * The fugacity coefficient \f$\mathrm{\phi_\kappa_\alpha}\f$ is connected to the
+     * fugacity \f$\mathrm{f^\kappa_\alpha}\f$ and the component's mole
+     * fraction \f$\mathrm{x^\kappa_\alpha}\f$ by means of the relation
      *
-     * \f[ f_\kappa = \phi_\kappa * x_{\kappa} \f]
+     * \f[
+     f^\kappa_\alpha = \phi^\kappa_\alpha\;x^\kappa_\alpha\;p_\alpha
+     * \f]
      */
     template <class FluidState>
     static Scalar fugacityCoefficient(const FluidState &fluidState,
diff --git a/dumux/material/fluidsystems/2pimmiscible.hh b/dumux/material/fluidsystems/2pimmiscible.hh
index bd7751f1d4aa81684db91e2472bf4a12d38f0c0b..7c5f3b837ad5de64d8a1bec9bd2e0e63721120e1 100644
--- a/dumux/material/fluidsystems/2pimmiscible.hh
+++ b/dumux/material/fluidsystems/2pimmiscible.hh
@@ -293,17 +293,20 @@ public:
 
     using Base::fugacityCoefficient;
     /*!
-     * \brief Calculate the fugacity coefficient \f$\mathrm{[Pa]}\f$ of an individual
+     * \brief Calculate the fugacity coefficient \f$\mathrm{[-]}\f$ of an individual
      *        component in a fluid phase
+     *
+     * The fugacity coefficient \f$\mathrm{\phi_\kappa_\alpha}\f$ is connected to the
+     * fugacity \f$\mathrm{f^\kappa_\alpha}\f$ and the component's mole
+     * fraction \f$\mathrm{x^\kappa_\alpha}\f$ by means of the relation
+     *
+     * \f[
+     f^\kappa_\alpha = \phi^\kappa_\alpha\;x^\kappa_\alpha\;p_\alpha
+     * \f]
+     *
      * \param fluidState The fluid state of the two-phase model
      * \param phaseIdx Index of the fluid phase
      * \param compIdx index of the component
-     *
-     * The fugacity coefficient \f$\mathrm{\phi_\kappa}\f$ is connected to the
-     * fugacity \f$\mathrm{f_\kappa}\f$ and the component's molarity
-     * \f$\mathrm{x_\kappa}\f$ by means of the relation
-     *
-     * \f[ f_\kappa = \phi_\kappa * x_{\kappa} \f]
      */
     template <class FluidState>
     static Scalar fugacityCoefficient(const FluidState &fluidState,
diff --git a/dumux/material/fluidsystems/2pliquidvapor.hh b/dumux/material/fluidsystems/2pliquidvapor.hh
index 869110fd324c1097c4592883cb40290e3c65c9f5..7950293f7426eb0d2e9def6c044353426cfbc288 100644
--- a/dumux/material/fluidsystems/2pliquidvapor.hh
+++ b/dumux/material/fluidsystems/2pliquidvapor.hh
@@ -377,7 +377,7 @@ public:
     }
 
     /*!
-     * \brief Calculate the fugacity coefficient [Pa] of an individual
+     * \brief Calculate the fugacity coefficient [-] of an individual
      *        component in a fluid phase
      *
      * The fugacity coefficient \f$\phi^\kappa_\alpha\f$ of
diff --git a/dumux/material/fluidsystems/CMakeLists.txt b/dumux/material/fluidsystems/CMakeLists.txt
index 43b27a7c2bebc65d416487c4803425dba1cc8f82..82c2e6159217cfda43c80af2fb524ac3f0f40175 100644
--- a/dumux/material/fluidsystems/CMakeLists.txt
+++ b/dumux/material/fluidsystems/CMakeLists.txt
@@ -10,6 +10,7 @@ install(FILES
   h2oair.hh
   h2oairmesitylene.hh
   h2oairxylene.hh
+  h2oheavyoil.hh
   h2oheavyoilfluidsystem.hh
   h2on2.hh
   h2on2kinetic.hh
diff --git a/dumux/material/fluidsystems/base.hh b/dumux/material/fluidsystems/base.hh
index d7df3931919eef749c31a03c4c9c95d8ea993ddf..10a0285b1365422bee7022ab5e832f0b81d452ca 100644
--- a/dumux/material/fluidsystems/base.hh
+++ b/dumux/material/fluidsystems/base.hh
@@ -61,18 +61,21 @@ public:
     }
 
     /*!
-     * \brief Calculate the fugacity coefficient \f$\mathrm{[Pa]}\f$ of an individual
+     * \brief Calculate the fugacity coefficient \f$\mathrm{[-]}\f$ of an individual
      *        component in a fluid phase
+     *
+     * The fugacity coefficient \f$\mathrm{\phi_\kappa_\alpha}\f$ is connected to the
+     * fugacity \f$\mathrm{f^\kappa_\alpha}\f$ and the component's mole
+     * fraction \f$\mathrm{x^\kappa_\alpha}\f$ by means of the relation
+     *
+     * \f[
+     f^\kappa_\alpha = \phi^\kappa_\alpha\;x^\kappa_\alpha\;p_\alpha
+     * \f]
+     *
      * \param fluidState The fluid state
      * \param paramCache mutable parameters
      * \param phaseIdx Index of the fluid phase
      * \param compIdx Index of the component
-     *
-     * The fugacity coefficient \f$\mathrm{\phi_\kappa}\f$ is connected to the
-     * fugacity \f$\mathrm{f_\kappa}\f$ and the component's molarity
-     * \f$\mathrm{x_\kappa}\f$ by means of the relation
-     *
-     * \f[ f_\kappa = \phi_\kappa * x_{\kappa} \f]
      */
     template <class FluidState>
     static Scalar fugacityCoefficient(const FluidState &fluidState,
diff --git a/dumux/material/fluidsystems/brineair.hh b/dumux/material/fluidsystems/brineair.hh
index 2e80056cdc0b58974232821f383812f0fe5596aa..7071f69f0e740e357a67a10e25c92f9765b07ef8 100644
--- a/dumux/material/fluidsystems/brineair.hh
+++ b/dumux/material/fluidsystems/brineair.hh
@@ -436,7 +436,7 @@ public:
      * where \f$\mathrm{p_\alpha}\f$ is the pressure of the fluid phase.
      *
      * For liquids with very low miscibility this boils down to the
-     * inverse Henry constant for the solutes and the saturated vapor pressure
+     * Henry constant for the solutes and the saturated vapor pressure
      * both divided by phase pressure.
      */
     template <class FluidState>
diff --git a/dumux/material/fluidsystems/h2oair.hh b/dumux/material/fluidsystems/h2oair.hh
index 7b3066d221fb939ec245e67310264d33b919d810..b2bd16f5a9db42e382db3f9ff07a612b70e3b964 100644
--- a/dumux/material/fluidsystems/h2oair.hh
+++ b/dumux/material/fluidsystems/h2oair.hh
@@ -547,7 +547,7 @@ public:
      * where \f$p_\alpha\f$ is the pressure of the fluid phase.
      *
      * For liquids with very low miscibility this boils down to the
-     * inverse Henry constant for the solutes and the saturated vapor pressure
+     * Henry constant for the solutes and the saturated vapor pressure
      * both divided by phase pressure.
      */
     template <class FluidState>
diff --git a/dumux/material/fluidsystems/h2oairmesitylene.hh b/dumux/material/fluidsystems/h2oairmesitylene.hh
index a0eb42bb7a85750ace4b373c7c27613b39a821e6..4e13346599f27758c3d0092a93c2f655c92b71f5 100644
--- a/dumux/material/fluidsystems/h2oairmesitylene.hh
+++ b/dumux/material/fluidsystems/h2oairmesitylene.hh
@@ -428,9 +428,9 @@ public:
      *
      * In this case, things are actually pretty simple. We have an ideal
      * solution. Thus, the fugacity coefficient is 1 in the gas phase
-     * (fugacity equals the partial pressure of the component in the gas phase
-     * respectively in the liquid phases it is the inverse of the
-     * Henry coefficients scaled by pressure
+     * (fugacity equals the partial pressure of the component in the gas phase)
+     * respectively in the liquid phases it is the Henry coefficients divided
+     * by pressure.
      * \param fluidState The fluid state
      * \param phaseIdx The index of the phase
      * \param compIdx The index of the component
diff --git a/dumux/material/fluidsystems/h2oairxylene.hh b/dumux/material/fluidsystems/h2oairxylene.hh
index b65f5440aac3552f74f9fe66f60059250664d508..edb77e44f488b5534fec39f3b3cfc763ef15e733 100644
--- a/dumux/material/fluidsystems/h2oairxylene.hh
+++ b/dumux/material/fluidsystems/h2oairxylene.hh
@@ -430,9 +430,9 @@ public:
      *
      * In this case, things are actually pretty simple. We have an ideal
      * solution. Thus, the fugacity coefficient is 1 in the gas phase
-     * (fugacity equals the partial pressure of the component in the gas phase
-     * respectively in the liquid phases it is the inverse of the
-     * Henry coefficients scaled by pressure
+     * (fugacity equals the partial pressure of the component in the gas phase)
+     * respectively in the liquid phases it is the Henry coefficients
+     * divided by pressure.
      */
     template <class FluidState>
     static Scalar fugacityCoefficient(const FluidState &fluidState,
diff --git a/dumux/material/fluidsystems/h2oheavyoil.hh b/dumux/material/fluidsystems/h2oheavyoil.hh
new file mode 100644
index 0000000000000000000000000000000000000000..bd811f3b7e3d3bc77f21be7a5bbcd36787813919
--- /dev/null
+++ b/dumux/material/fluidsystems/h2oheavyoil.hh
@@ -0,0 +1,491 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A fluid system with water, gas and NAPL as phases and
+ *        \f$H_2O\f$ and \f$NAPL (contaminant)\f$ as components.
+ * It uses heavyoil Properties, but allows for a scaling of density
+ * thus enabling an artifical DNAPL if desired
+ */
+#ifndef DUMUX_H2O_HEAVYOIL_FLUID_SYSTEM_HH
+#define DUMUX_H2O_HEAVYOIL_FLUID_SYSTEM_HH
+
+#include <dumux/material/idealgas.hh>
+#include <dumux/material/components/h2o.hh>
+#include <dumux/material/components/tabulatedcomponent.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/heavyoil.hh>
+
+#include <dumux/material/binarycoefficients/h2o_heavyoil.hh>
+
+#include <dumux/material/fluidsystems/base.hh>
+
+namespace Dumux
+{
+namespace FluidSystems
+{
+
+/*!
+ * \ingroup Fluidsystems
+ * \brief A compositional fluid with water and heavy oil
+ *        components in both, the liquid and the gas phase.
+ */
+template <class TypeTag, class Scalar,
+          class H2OType = Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar> > >
+class H2OHeavyOil
+    : public BaseFluidSystem<Scalar, H2OHeavyOil<TypeTag, Scalar, H2OType> >
+{
+    typedef H2OHeavyOil<TypeTag, Scalar, H2OType> ThisType;
+    typedef BaseFluidSystem<Scalar, ThisType> Base;
+
+public:
+    typedef Dumux::HeavyOil<Scalar> HeavyOil;
+    typedef H2OType H2O;
+
+
+    static const int numPhases = 3;
+    static const int numComponents = 2;
+
+    static const int wPhaseIdx = 0; // index of the water phase
+    static const int nPhaseIdx = 1; // index of the NAPL phase
+    static const int gPhaseIdx = 2; // index of the gas phase
+
+    static const int H2OIdx = 0;
+    static const int NAPLIdx = 1;
+
+    // export component indices to indicate the main component
+    // of the corresponding phase at atmospheric pressure 1 bar
+    // and room temperature 20°C:
+    static const int wCompIdx = H2OIdx;
+    static const int nCompIdx = NAPLIdx;
+
+    /*!
+     * \brief Initialize the fluid system's static parameters generically
+     *
+     * If a tabulated H2O component is used, we do our best to create
+     * tables that always work.
+     */
+    static void init()
+    {
+        init(/*tempMin=*/273.15,
+             /*tempMax=*/623.15,
+             /*numTemp=*/100,
+             /*pMin=*/0.0,
+             /*pMax=*/20e6,
+             /*numP=*/200);
+    }
+
+    /*!
+     * \brief Initialize the fluid system's static parameters using
+     *        problem specific temperature and pressure ranges
+     *
+     * \param tempMin The minimum temperature used for tabulation of water [K]
+     * \param tempMax The maximum temperature used for tabulation of water [K]
+     * \param nTemp The number of ticks on the temperature axis of the  table of water
+     * \param pressMin The minimum pressure used for tabulation of water [Pa]
+     * \param pressMax The maximum pressure used for tabulation of water [Pa]
+     * \param nPress The number of ticks on the pressure axis of the  table of water
+     */
+    static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
+                     Scalar pressMin, Scalar pressMax, unsigned nPress)
+    {
+        if (H2O::isTabulated) {
+            std::cout << "Initializing tables for the H2O fluid properties ("
+                      << nTemp*nPress
+                      << " entries).\n";
+
+            H2O::init(tempMin, tempMax, nTemp,
+                      pressMin, pressMax, nPress);
+        }
+    }
+
+
+    /*!
+     * \brief Return whether a phase is liquid
+     *
+     * \param phaseIdx The index of the fluid phase to consider
+     */
+    static bool isLiquid(int phaseIdx)
+    {
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+        return phaseIdx != gPhaseIdx;
+    }
+
+    static bool isIdealGas(int phaseIdx)
+    { return phaseIdx == gPhaseIdx && H2O::gasIsIdeal() && HeavyOil::gasIsIdeal(); }
+
+    /*!
+     * \brief Returns true if and only if a fluid phase is assumed to
+     *        be an ideal mixture.
+     *
+     * We define an ideal mixture as a fluid phase where the fugacity
+     * coefficients of all components times the pressure of the phase
+     * are indepent on the fluid composition. This assumtion is true
+     * if Henry's law and Raoult's law apply. If you are unsure what
+     * this function should return, it is safe to return false. The
+     * only damage done will be (slightly) increased computation times
+     * in some cases.
+     *
+     * \param phaseIdx The index of the fluid phase to consider
+     */
+    static bool isIdealMixture(int phaseIdx)
+    {
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+        // we assume Henry's and Raoult's laws for the water phase and
+        // and no interaction between gas molecules of different
+        // components, so all phases are ideal mixtures!
+        return true;
+    }
+
+    /*!
+     * \brief Returns true if and only if a fluid phase is assumed to
+     *        be compressible.
+     *
+     * Compressible means that the partial derivative of the density
+     * to the fluid pressure is always larger than zero.
+     *
+     * \param phaseIdx The index of the fluid phase to consider
+     */
+    static bool isCompressible(int phaseIdx)
+    {
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+        // gases are always compressible
+        if (phaseIdx == gPhaseIdx)
+            return true;
+        else if (phaseIdx == wPhaseIdx)
+            // the water component decides for the water phase...
+            return H2O::liquidIsCompressible();
+
+        // the NAPL component decides for the napl phase...
+        return HeavyOil::liquidIsCompressible();
+    }
+
+    /*!
+     * \brief Return the human readable name of a phase (used in indices)
+     */
+    static std::string phaseName(int phaseIdx)
+    {
+        switch (phaseIdx) {
+        case wPhaseIdx: return "w";
+        case nPhaseIdx: return "n";
+        case gPhaseIdx: return "g";;
+        };
+        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+    }
+
+    /*!
+     * \brief Return the human readable name of a component (used in indices)
+     */
+    static std::string componentName(int compIdx)
+    {
+        switch (compIdx) {
+        case H2OIdx: return H2O::name();
+        case NAPLIdx: return HeavyOil::name();
+        };
+        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+    }
+
+    /*!
+     * \brief Return the molar mass of a component in [kg/mol].
+     */
+    static Scalar molarMass(int compIdx)
+    {
+        switch (compIdx) {
+        case H2OIdx: return H2O::molarMass();
+        case NAPLIdx: return HeavyOil::molarMass();
+        };
+        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+    }
+
+    /*!
+     * \brief Given all mole fractions in a phase, return the phase
+     *        density [kg/m^3].
+     */
+    using Base::density;
+    template <class FluidState>
+    static Scalar density(const FluidState &fluidState, int phaseIdx)
+    {
+        if (phaseIdx == wPhaseIdx) {
+            // See: doctoral thesis of Steffen Ochs 2007
+            // Steam injection into saturated porous media : process analysis including experimental and numerical investigations
+            // http://elib.uni-stuttgart.de/bitstream/11682/271/1/Diss_Ochs_OPUS.pdf
+            Scalar rholH2O = H2O::liquidDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+            Scalar clH2O = rholH2O/H2O::molarMass();
+
+            // this assumes each dissolved molecule displaces exactly one
+            // water molecule in the liquid
+            return
+                clH2O*(H2O::molarMass()*fluidState.moleFraction(wPhaseIdx, H2OIdx)
+                       +
+                       HeavyOil::molarMass()*fluidState.moleFraction(wPhaseIdx, NAPLIdx));
+        }
+        else if (phaseIdx == nPhaseIdx) {
+            // assume pure NAPL for the NAPL phase
+            Scalar pressure = HeavyOil::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
+            return HeavyOil::liquidDensity(fluidState.temperature(phaseIdx), pressure);
+        }
+
+        assert (phaseIdx == gPhaseIdx);
+        Scalar pH2O =
+            fluidState.moleFraction(gPhaseIdx, H2OIdx)  *
+            fluidState.pressure(gPhaseIdx);
+        Scalar pNAPL =
+            fluidState.moleFraction(gPhaseIdx, NAPLIdx)  *
+            fluidState.pressure(gPhaseIdx);
+        return
+            H2O::gasDensity(fluidState.temperature(phaseIdx), pH2O) +
+            HeavyOil::gasDensity(fluidState.temperature(phaseIdx), pNAPL);
+    }
+
+    /*!
+     * \brief Return the viscosity of a phase.
+     */
+    using Base::viscosity;
+    template <class FluidState>
+    static Scalar viscosity(const FluidState &fluidState,
+                            int phaseIdx)
+    {
+        if (phaseIdx == wPhaseIdx) {
+            // assume pure water viscosity
+            return H2O::liquidViscosity(fluidState.temperature(phaseIdx),
+                                        fluidState.pressure(phaseIdx));
+        }
+        else if (phaseIdx == nPhaseIdx) {
+            // assume pure NAPL viscosity
+            return HeavyOil::liquidViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+        }
+
+        assert (phaseIdx == gPhaseIdx);
+
+        /* Wilke method. See:
+         *
+         * See: R. Reid, et al.: The Properties of Gases and Liquids,
+         * 4th edition, McGraw-Hill, 1987, 407-410
+         * 5th edition, McGraw-Hill, 20001, p. 9.21/22
+         *
+         * in this case, we use a simplified version in order to avoid
+         * computationally costly evaluation of sqrt and pow functions and
+         * divisions
+         * -- compare e.g. with Promo Class p. 32/33
+         */
+        const Scalar mu[numComponents] = {
+            H2O::gasViscosity(fluidState.temperature(phaseIdx), H2O::vaporPressure(fluidState.temperature(phaseIdx))),
+            HeavyOil::gasViscosity(fluidState.temperature(phaseIdx), HeavyOil::vaporPressure(fluidState.temperature(phaseIdx)))
+        };
+
+        return mu[H2OIdx]*fluidState.moleFraction(gPhaseIdx, H2OIdx)
+                + mu[NAPLIdx]*fluidState.moleFraction(gPhaseIdx, NAPLIdx);
+    }
+
+
+    using Base::binaryDiffusionCoefficient;
+    template <class FluidState>
+    static Scalar diffusionCoefficient(const FluidState &fluidState, int phaseIdx)
+    {
+        const Scalar T = fluidState.temperature(phaseIdx);
+        const Scalar p = fluidState.pressure(phaseIdx);
+
+        // liquid phase
+        if (phaseIdx == wPhaseIdx)
+            return BinaryCoeff::H2O_HeavyOil::liquidDiffCoeff(T, p);
+
+        // gas phase
+        else if (phaseIdx == gPhaseIdx)
+            return BinaryCoeff::H2O_HeavyOil::gasDiffCoeff(T, p);
+
+        else
+            DUNE_THROW(Dune::InvalidStateException, "non-existent diffusion coefficient for phase index " << phaseIdx);
+    }
+
+     /*!
+     * \brief Henry coefficients
+     */
+    template <class FluidState>
+    static Scalar henryCoefficient(const FluidState &fluidState,
+                                   int phaseIdx,
+                                   int compIdx)
+    {
+        assert(0 <= phaseIdx  && phaseIdx < numPhases);
+        assert(0 <= compIdx  && compIdx < numComponents);
+
+        const Scalar T = fluidState.temperature(phaseIdx);
+
+        if (compIdx == NAPLIdx && phaseIdx == wPhaseIdx)
+            return Dumux::BinaryCoeff::H2O_HeavyOil::henryOilInWater(T);
+
+        else if (phaseIdx == nPhaseIdx && compIdx == H2OIdx)
+            return Dumux::BinaryCoeff::H2O_HeavyOil::henryWaterInOil(T);
+
+        else
+            DUNE_THROW(Dune::InvalidStateException, "non-existent henry coefficient for phase index " << phaseIdx
+                                                     << " and component index " << compIdx);
+    }
+
+     /*!
+     * \brief Partial pressures in the gas phase, taken from saturation vapor pressures
+     */
+    template <class FluidState>
+    static Scalar partialPressureGas(const FluidState &fluidState, int phaseIdx,
+                                      int compIdx)
+    {
+        assert(0 <= compIdx  && compIdx < numComponents);
+
+        const Scalar T = fluidState.temperature(phaseIdx);
+        if (compIdx == NAPLIdx)
+            return HeavyOil::vaporPressure(T);
+        else if (compIdx == H2OIdx)
+            return H2O::vaporPressure(T);
+        else
+            DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
+    }
+
+     /*!
+     * \brief Inverse vapor pressures, taken from inverse saturation vapor pressures
+     */
+    template <class FluidState>
+    static Scalar inverseVaporPressureCurve(const FluidState &fluidState,
+                                            int phaseIdx,
+                                            int compIdx)
+    {
+        assert(0 <= compIdx  && compIdx < numComponents);
+
+        const Scalar pressure = fluidState.pressure(phaseIdx);
+        if (compIdx == NAPLIdx)
+            return HeavyOil::vaporTemperature(pressure);
+        else if (compIdx == H2OIdx)
+            return H2O::vaporTemperature(pressure);
+        else
+            DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
+    }
+
+
+
+    /*!
+     * \brief Given all mole fractions in a phase, return the specific
+     *        phase enthalpy [J/kg].
+     *
+     *  \todo This system neglects the contribution of gas-molecules in the liquid phase.
+     *        This contribution is probably not big. Somebody would have to find out the enthalpy of solution for this system. ...
+     */
+    using Base::enthalpy;
+    template <class FluidState>
+    static Scalar enthalpy(const FluidState &fluidState,
+                           int phaseIdx)
+    {
+        if (phaseIdx == wPhaseIdx) {
+            return H2O::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+        }
+        else if (phaseIdx == nPhaseIdx) {
+            return HeavyOil::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+        }
+        else if (phaseIdx == gPhaseIdx) {  // gas phase enthalpy depends strongly on composition
+            Scalar hgc = HeavyOil::gasEnthalpy(fluidState.temperature(phaseIdx),
+                                           fluidState.pressure(phaseIdx));
+            Scalar hgw = H2O::gasEnthalpy(fluidState.temperature(phaseIdx),
+                                          fluidState.pressure(phaseIdx));
+
+            Scalar result = 0;
+            result += hgw * fluidState.massFraction(gPhaseIdx, H2OIdx);
+            result += hgc * fluidState.massFraction(gPhaseIdx, NAPLIdx);
+
+            return result;
+        }
+        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+    }
+
+    using Base::heatCapacity;
+    template <class FluidState>
+    static Scalar heatCapacity(const FluidState &fluidState,
+                               int phaseIdx)
+    {
+        DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2ONAPL::heatCapacity()");
+    }
+
+    using Base::thermalConductivity;
+    template <class FluidState>
+    static Scalar thermalConductivity(const FluidState &fluidState,
+                                      int phaseIdx)
+    {
+        const Scalar temperature  = fluidState.temperature(phaseIdx) ;
+        const Scalar pressure = fluidState.pressure(phaseIdx);
+        if (phaseIdx == wPhaseIdx)
+        {
+            return H2O::liquidThermalConductivity(temperature, pressure);
+        }
+        else if (phaseIdx == nPhaseIdx)
+        {
+            return HeavyOil::liquidThermalConductivity(temperature, pressure);
+        }
+        else if (phaseIdx == gPhaseIdx)
+        {
+            return H2O::gasThermalConductivity(temperature, pressure);
+        }
+        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+    }
+
+private:
+
+};
+} // end namespace FluidSystems
+
+#ifdef DUMUX_PROPERTIES_HH
+// forward definitions of the property tags
+namespace Properties {
+    NEW_PROP_TAG(Scalar);
+    NEW_PROP_TAG(Components);
+}
+
+/*!
+ * \brief A threephase fluid system with water, heavyoil as components.
+ *
+ * This is an adapter to use Dumux::H2OheavyoilFluidSystem<TypeTag>, as is
+ * done with most other classes in Dumux.
+ *  This fluidsystem is applied by default with the tabulated version of
+ *  water of the IAPWS-formulation.
+ *
+ *  To change the component formulation (ie to use nontabulated or
+ *  incompressible water), or to switch on verbosity of tabulation,
+ *  use the property system and the property "Components":
+ *
+ *        // Select desired version of the component
+ *        SET_PROP(myApplicationProperty, Components) : public GET_PROP(TypeTag, DefaultComponents)
+ *        {
+ *            typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ *
+ *        // Do not use the defaults !
+ *        //    typedef Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar> > H2O;
+ *
+ *        // Apply e.g. untabulated water:
+ *        typedef Dumux::H2O<Scalar> H2O;
+ *        };
+ *
+ *   Also remember to initialize tabulated components (FluidSystem::init()), while this
+ *   is not necessary for non-tabularized ones.
+ */
+template<class TypeTag>
+class H2OHeavyOilFluidSystem
+: public FluidSystems::H2OHeavyOil<TypeTag, typename GET_PROP_TYPE(TypeTag, Scalar),
+                                        typename GET_PROP(TypeTag, Components)::H2O>
+{};
+#endif
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/material/fluidsystems/h2oheavyoilfluidsystem.hh b/dumux/material/fluidsystems/h2oheavyoilfluidsystem.hh
index e08409edfdb5b2e5912bd6c94fe9b064af0c7946..3a2a90a5f2fb1843b9f7d5f6dfb303cf77cec51f 100644
--- a/dumux/material/fluidsystems/h2oheavyoilfluidsystem.hh
+++ b/dumux/material/fluidsystems/h2oheavyoilfluidsystem.hh
@@ -1,489 +1,7 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program is free software: you can redistribute it and/or modify    *
- *   it under the terms of the GNU General Public License as published by    *
- *   the Free Software Foundation, either version 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program is distributed in the hope that it will be useful,         *
- *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
- *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
- *   GNU General Public License for more details.                            *
- *                                                                           *
- *   You should have received a copy of the GNU General Public License       *
- *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief A fluid system with water, gas and NAPL as phases and
- *        \f$H_2O\f$ and \f$NAPL (contaminant)\f$ as components.
- * It uses heavyoil Properties, but allows for a scaling of density
- * thus enabling an artifical DNAPL if desired
- */
-#ifndef DUMUX_H2O_HEAVYOIL_FLUID_SYSTEM_HH
-#define DUMUX_H2O_HEAVYOIL_FLUID_SYSTEM_HH
+#ifndef DUMUX_H2O_HEAVYOIL_FLUID_SYSTEM_OLD_HH
+#define DUMUX_H2O_HEAVYOIL_FLUID_SYSTEM_OLD_HH
 
-#include <dumux/material/idealgas.hh>
-#include <dumux/material/components/h2o.hh>
-#include <dumux/material/components/tabulatedcomponent.hh>
-#include <dumux/material/components/simpleh2o.hh>
-#include <dumux/material/components/heavyoil.hh>
-
-#include <dumux/material/binarycoefficients/h2o_heavyoil.hh>
-
-#include <dumux/material/fluidsystems/base.hh>
-
-namespace Dumux
-{
-namespace FluidSystems
-{
-
-/*!
- * \brief A compositional fluid with water and heavy oil
- *        components in both, the liquid and the gas phase.
- */
-template <class TypeTag, class Scalar,
-          class H2OType = Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar> > >
-class H2OHeavyOil
-    : public BaseFluidSystem<Scalar, H2OHeavyOil<TypeTag, Scalar, H2OType> >
-{
-    typedef H2OHeavyOil<TypeTag, Scalar, H2OType> ThisType;
-    typedef BaseFluidSystem<Scalar, ThisType> Base;
-
-public:
-    typedef Dumux::HeavyOil<Scalar> HeavyOil;
-    typedef H2OType H2O;
-
-
-    static const int numPhases = 3;
-    static const int numComponents = 2;
-
-    static const int wPhaseIdx = 0; // index of the water phase
-    static const int nPhaseIdx = 1; // index of the NAPL phase
-    static const int gPhaseIdx = 2; // index of the gas phase
-
-    static const int H2OIdx = 0;
-    static const int NAPLIdx = 1;
-
-    // export component indices to indicate the main component
-    // of the corresponding phase at atmospheric pressure 1 bar
-    // and room temperature 20°C:
-    static const int wCompIdx = H2OIdx;
-    static const int nCompIdx = NAPLIdx;
-
-    /*!
-     * \brief Initialize the fluid system's static parameters generically
-     *
-     * If a tabulated H2O component is used, we do our best to create
-     * tables that always work.
-     */
-    static void init()
-    {
-        init(/*tempMin=*/273.15,
-             /*tempMax=*/623.15,
-             /*numTemp=*/100,
-             /*pMin=*/0.0,
-             /*pMax=*/20e6,
-             /*numP=*/200);
-    }
-
-    /*!
-     * \brief Initialize the fluid system's static parameters using
-     *        problem specific temperature and pressure ranges
-     *
-     * \param tempMin The minimum temperature used for tabulation of water [K]
-     * \param tempMax The maximum temperature used for tabulation of water [K]
-     * \param nTemp The number of ticks on the temperature axis of the  table of water
-     * \param pressMin The minimum pressure used for tabulation of water [Pa]
-     * \param pressMax The maximum pressure used for tabulation of water [Pa]
-     * \param nPress The number of ticks on the pressure axis of the  table of water
-     */
-    static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
-                     Scalar pressMin, Scalar pressMax, unsigned nPress)
-    {
-        if (H2O::isTabulated) {
-            std::cout << "Initializing tables for the H2O fluid properties ("
-                      << nTemp*nPress
-                      << " entries).\n";
-
-            H2O::init(tempMin, tempMax, nTemp,
-                      pressMin, pressMax, nPress);
-        }
-    }
-
-
-    /*!
-     * \brief Return whether a phase is liquid
-     *
-     * \param phaseIdx The index of the fluid phase to consider
-     */
-    static bool isLiquid(int phaseIdx)
-    {
-        assert(0 <= phaseIdx && phaseIdx < numPhases);
-        return phaseIdx != gPhaseIdx;
-    }
-
-    static bool isIdealGas(int phaseIdx)
-    { return phaseIdx == gPhaseIdx && H2O::gasIsIdeal() && HeavyOil::gasIsIdeal(); }
-
-    /*!
-     * \brief Returns true if and only if a fluid phase is assumed to
-     *        be an ideal mixture.
-     *
-     * We define an ideal mixture as a fluid phase where the fugacity
-     * coefficients of all components times the pressure of the phase
-     * are indepent on the fluid composition. This assumtion is true
-     * if Henry's law and Raoult's law apply. If you are unsure what
-     * this function should return, it is safe to return false. The
-     * only damage done will be (slightly) increased computation times
-     * in some cases.
-     *
-     * \param phaseIdx The index of the fluid phase to consider
-     */
-    static bool isIdealMixture(int phaseIdx)
-    {
-        assert(0 <= phaseIdx && phaseIdx < numPhases);
-        // we assume Henry's and Raoult's laws for the water phase and
-        // and no interaction between gas molecules of different
-        // components, so all phases are ideal mixtures!
-        return true;
-    }
-
-    /*!
-     * \brief Returns true if and only if a fluid phase is assumed to
-     *        be compressible.
-     *
-     * Compressible means that the partial derivative of the density
-     * to the fluid pressure is always larger than zero.
-     *
-     * \param phaseIdx The index of the fluid phase to consider
-     */
-    static bool isCompressible(int phaseIdx)
-    {
-        assert(0 <= phaseIdx && phaseIdx < numPhases);
-        // gases are always compressible
-        if (phaseIdx == gPhaseIdx)
-            return true;
-        else if (phaseIdx == wPhaseIdx)
-            // the water component decides for the water phase...
-            return H2O::liquidIsCompressible();
-
-        // the NAPL component decides for the napl phase...
-        return HeavyOil::liquidIsCompressible();
-    }
-
-    /*!
-     * \brief Return the human readable name of a phase (used in indices)
-     */
-    static std::string phaseName(int phaseIdx)
-    {
-        switch (phaseIdx) {
-        case wPhaseIdx: return "w";
-        case nPhaseIdx: return "n";
-        case gPhaseIdx: return "g";;
-        };
-        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
-    }
-
-    /*!
-     * \brief Return the human readable name of a component (used in indices)
-     */
-    static std::string componentName(int compIdx)
-    {
-        switch (compIdx) {
-        case H2OIdx: return H2O::name();
-        case NAPLIdx: return HeavyOil::name();
-        };
-        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
-    }
-
-    /*!
-     * \brief Return the molar mass of a component in [kg/mol].
-     */
-    static Scalar molarMass(int compIdx)
-    {
-        switch (compIdx) {
-        case H2OIdx: return H2O::molarMass();
-        case NAPLIdx: return HeavyOil::molarMass();
-        };
-        DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
-    }
-
-    /*!
-     * \brief Given all mole fractions in a phase, return the phase
-     *        density [kg/m^3].
-     */
-    using Base::density;
-    template <class FluidState>
-    static Scalar density(const FluidState &fluidState, int phaseIdx)
-    {
-        if (phaseIdx == wPhaseIdx) {
-            // See: doctoral thesis of Steffen Ochs 2007
-            // Steam injection into saturated porous media : process analysis including experimental and numerical investigations
-            // http://elib.uni-stuttgart.de/bitstream/11682/271/1/Diss_Ochs_OPUS.pdf
-            Scalar rholH2O = H2O::liquidDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
-            Scalar clH2O = rholH2O/H2O::molarMass();
-
-            // this assumes each dissolved molecule displaces exactly one
-            // water molecule in the liquid
-            return
-                clH2O*(H2O::molarMass()*fluidState.moleFraction(wPhaseIdx, H2OIdx)
-                       +
-                       HeavyOil::molarMass()*fluidState.moleFraction(wPhaseIdx, NAPLIdx));
-        }
-        else if (phaseIdx == nPhaseIdx) {
-            // assume pure NAPL for the NAPL phase
-            Scalar pressure = HeavyOil::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
-            return HeavyOil::liquidDensity(fluidState.temperature(phaseIdx), pressure);
-        }
-
-        assert (phaseIdx == gPhaseIdx);
-        Scalar pH2O =
-            fluidState.moleFraction(gPhaseIdx, H2OIdx)  *
-            fluidState.pressure(gPhaseIdx);
-        Scalar pNAPL =
-            fluidState.moleFraction(gPhaseIdx, NAPLIdx)  *
-            fluidState.pressure(gPhaseIdx);
-        return
-            H2O::gasDensity(fluidState.temperature(phaseIdx), pH2O) +
-            HeavyOil::gasDensity(fluidState.temperature(phaseIdx), pNAPL);
-    }
-
-    /*!
-     * \brief Return the viscosity of a phase.
-     */
-    using Base::viscosity;
-    template <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
-                            int phaseIdx)
-    {
-        if (phaseIdx == wPhaseIdx) {
-            // assume pure water viscosity
-            return H2O::liquidViscosity(fluidState.temperature(phaseIdx),
-                                        fluidState.pressure(phaseIdx));
-        }
-        else if (phaseIdx == nPhaseIdx) {
-            // assume pure NAPL viscosity
-            return HeavyOil::liquidViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
-        }
-
-        assert (phaseIdx == gPhaseIdx);
-
-        /* Wilke method. See:
-         *
-         * See: R. Reid, et al.: The Properties of Gases and Liquids,
-         * 4th edition, McGraw-Hill, 1987, 407-410
-         * 5th edition, McGraw-Hill, 20001, p. 9.21/22
-         *
-         * in this case, we use a simplified version in order to avoid
-         * computationally costly evaluation of sqrt and pow functions and
-         * divisions
-         * -- compare e.g. with Promo Class p. 32/33
-         */
-        const Scalar mu[numComponents] = {
-            H2O::gasViscosity(fluidState.temperature(phaseIdx), H2O::vaporPressure(fluidState.temperature(phaseIdx))),
-            HeavyOil::gasViscosity(fluidState.temperature(phaseIdx), HeavyOil::vaporPressure(fluidState.temperature(phaseIdx)))
-        };
-
-        return mu[H2OIdx]*fluidState.moleFraction(gPhaseIdx, H2OIdx)
-                + mu[NAPLIdx]*fluidState.moleFraction(gPhaseIdx, NAPLIdx);
-    }
-
-
-    using Base::binaryDiffusionCoefficient;
-    template <class FluidState>
-    static Scalar diffusionCoefficient(const FluidState &fluidState, int phaseIdx)
-    {
-        const Scalar T = fluidState.temperature(phaseIdx);
-        const Scalar p = fluidState.pressure(phaseIdx);
-
-        // liquid phase
-        if (phaseIdx == wPhaseIdx)
-            return BinaryCoeff::H2O_HeavyOil::liquidDiffCoeff(T, p);
-
-        // gas phase
-        else if (phaseIdx == gPhaseIdx)
-            return BinaryCoeff::H2O_HeavyOil::gasDiffCoeff(T, p);
-
-        else
-            DUNE_THROW(Dune::InvalidStateException, "non-existent diffusion coefficient for phase index " << phaseIdx);
-    }
-
-     /* Henry coefficients
-     */
-    template <class FluidState>
-    static Scalar henryCoefficient(const FluidState &fluidState,
-                                   int phaseIdx,
-                                   int compIdx)
-    {
-        assert(0 <= phaseIdx  && phaseIdx < numPhases);
-        assert(0 <= compIdx  && compIdx < numComponents);
-
-        const Scalar T = fluidState.temperature(phaseIdx);
-        const Scalar p = fluidState.pressure(phaseIdx);
-
-        if (compIdx == NAPLIdx && phaseIdx == wPhaseIdx)
-            return Dumux::BinaryCoeff::H2O_HeavyOil::henryOilInWater(T)/p;
-
-        else if (phaseIdx == nPhaseIdx && compIdx == H2OIdx)
-            return Dumux::BinaryCoeff::H2O_HeavyOil::henryWaterInOil(T)/p;
-
-        else
-            DUNE_THROW(Dune::InvalidStateException, "non-existent henry coefficient for phase index " << phaseIdx
-                                                     << " and component index " << compIdx);
-    }
-
-     /*  partial pressures in the gas phase, taken from saturation vapor pressures
-     */
-    template <class FluidState>
-    static Scalar partialPressureGas(const FluidState &fluidState, int phaseIdx,
-                                      int compIdx)
-    {
-        assert(0 <= compIdx  && compIdx < numComponents);
-
-        const Scalar T = fluidState.temperature(phaseIdx);
-        if (compIdx == NAPLIdx)
-            return HeavyOil::vaporPressure(T);
-        else if (compIdx == H2OIdx)
-            return H2O::vaporPressure(T);
-        else
-            DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
-    }
-
-     /*  inverse vapor pressures, taken from inverse saturation vapor pressures
-     */
-    template <class FluidState>
-    static Scalar inverseVaporPressureCurve(const FluidState &fluidState,
-                                            int phaseIdx,
-                                            int compIdx)
-    {
-        assert(0 <= compIdx  && compIdx < numComponents);
-
-        const Scalar pressure = fluidState.pressure(phaseIdx);
-        if (compIdx == NAPLIdx)
-            return HeavyOil::vaporTemperature(pressure);
-        else if (compIdx == H2OIdx)
-            return H2O::vaporTemperature(pressure);
-        else
-            DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
-    }
-
-
-
-    /*!
-     * \brief Given all mole fractions in a phase, return the specific
-     *        phase enthalpy [J/kg].
-     */
-    /*!
-     *  \todo This system neglects the contribution of gas-molecules in the liquid phase.
-     *        This contribution is probably not big. Somebody would have to find out the enthalpy of solution for this system. ...
-     */
-    using Base::enthalpy;
-    template <class FluidState>
-    static Scalar enthalpy(const FluidState &fluidState,
-                           int phaseIdx)
-    {
-        if (phaseIdx == wPhaseIdx) {
-            return H2O::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
-        }
-        else if (phaseIdx == nPhaseIdx) {
-            return HeavyOil::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
-        }
-        else if (phaseIdx == gPhaseIdx) {  // gas phase enthalpy depends strongly on composition
-            Scalar hgc = HeavyOil::gasEnthalpy(fluidState.temperature(phaseIdx),
-                                           fluidState.pressure(phaseIdx));
-            Scalar hgw = H2O::gasEnthalpy(fluidState.temperature(phaseIdx),
-                                          fluidState.pressure(phaseIdx));
-
-            Scalar result = 0;
-            result += hgw * fluidState.massFraction(gPhaseIdx, H2OIdx);
-            result += hgc * fluidState.massFraction(gPhaseIdx, NAPLIdx);
-
-            return result;
-        }
-        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
-    }
-
-    using Base::heatCapacity;
-    template <class FluidState>
-    static Scalar heatCapacity(const FluidState &fluidState,
-                               int phaseIdx)
-    {
-        DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2ONAPL::heatCapacity()");
-    }
-
-    using Base::thermalConductivity;
-    template <class FluidState>
-    static Scalar thermalConductivity(const FluidState &fluidState,
-                                      int phaseIdx)
-    {
-        const Scalar temperature  = fluidState.temperature(phaseIdx) ;
-        const Scalar pressure = fluidState.pressure(phaseIdx);
-        if (phaseIdx == wPhaseIdx)
-        {
-            return H2O::liquidThermalConductivity(temperature, pressure);
-        }
-        else if (phaseIdx == nPhaseIdx)
-        {
-            return HeavyOil::liquidThermalConductivity(temperature, pressure);
-        }
-        else if (phaseIdx == gPhaseIdx)
-        {
-            return H2O::gasThermalConductivity(temperature, pressure);
-        }
-        DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
-    }
-
-private:
-
-};
-} // end namespace FluidSystems
-
-#ifdef DUMUX_PROPERTIES_HH
-// forward definitions of the property tags
-namespace Properties {
-    NEW_PROP_TAG(Scalar);
-    NEW_PROP_TAG(Components);
-}
-
-/*!
- * \brief A threephase fluid system with water, heavyoil as components.
- *
- * This is an adapter to use Dumux::H2OheavyoilFluidSystem<TypeTag>, as is
- * done with most other classes in Dumux.
- *  This fluidsystem is applied by default with the tabulated version of
- *  water of the IAPWS-formulation.
- *
- *  To change the component formulation (ie to use nontabulated or
- *  incompressible water), or to switch on verbosity of tabulation,
- *  use the property system and the property "Components":
- *
- *        // Select desired version of the component
- *        SET_PROP(myApplicationProperty, Components) : public GET_PROP(TypeTag, DefaultComponents)
- *        {
- *            typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- *
- *        // Do not use the defaults !
- *        //    typedef Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar> > H2O;
- *
- *        // Apply e.g. untabulated water:
- *        typedef Dumux::H2O<Scalar> H2O;
- *        };
- *
- *   Also remember to initialize tabulated components (FluidSystem::init()), while this
- *   is not necessary for non-tabularized ones.
- */
-template<class TypeTag>
-class H2OHeavyOilFluidSystem
-: public FluidSystems::H2OHeavyOil<TypeTag, typename GET_PROP_TYPE(TypeTag, Scalar),
-                                        typename GET_PROP(TypeTag, Components)::H2O>
-{};
-#endif
-
-} // end namespace Dumux
+#warning The file h2oheavyoilfluidsystem.hh is deprecated, use h2oheavyoil.hh instead.
+#include "h2oheavyoil.hh"
 
 #endif
diff --git a/dumux/material/fluidsystems/h2on2o2.hh b/dumux/material/fluidsystems/h2on2o2.hh
index 24f4bb2c694b3704429e3257ddb7c461874259a3..99d29746df090a2610406a5b3d0b388995eb3ec8 100644
--- a/dumux/material/fluidsystems/h2on2o2.hh
+++ b/dumux/material/fluidsystems/h2on2o2.hh
@@ -552,7 +552,7 @@ public:
      * where \f$p_\alpha\f$ is the pressure of the fluid phase.
      *
      * For liquids with very low miscibility this boils down to the
-     * inverse Henry constant for the solutes and the saturated vapor pressure
+     * Henry constant for the solutes and the saturated vapor pressure
      * both divided by phase pressure.
      */
     template <class FluidState>
diff --git a/dumux/material/fluidsystems/purewatersimple.hh b/dumux/material/fluidsystems/purewatersimple.hh
index 338bce8fd5104d2fddba639d3628413dcb60c0b2..8c51ebbb2b57169ba721cabc9317aab2d1276706 100644
--- a/dumux/material/fluidsystems/purewatersimple.hh
+++ b/dumux/material/fluidsystems/purewatersimple.hh
@@ -388,7 +388,7 @@ public:
 
     using Base::fugacityCoefficient;
     /*!
-     * \brief Calculate the fugacity coefficient \f$\mathrm{[Pa]}\f$ of an individual
+     * \brief Calculate the fugacity coefficient \f$\mathrm{[-]}\f$ of an individual
      *        component in a fluid phase
      *
      * The fugacity coefficient \f$\phi^\kappa_\alpha\f$ of
diff --git a/dumux/material/fluidsystems/spe5.hh b/dumux/material/fluidsystems/spe5.hh
index a8a679f0eb591f37db132a71190246044539bb37..57aacf47c0e4360fdecf687a6a8ba9c40d444399 100644
--- a/dumux/material/fluidsystems/spe5.hh
+++ b/dumux/material/fluidsystems/spe5.hh
@@ -415,7 +415,7 @@ public:
     }
 
     /*!
-     * \brief Calculate the fugacity coefficient \f$\mathrm{[Pa]}\f$ of an individual
+     * \brief Calculate the fugacity coefficient \f$\mathrm{[-]}\f$ of an individual
      *        component in a fluid phase
      *
      * The fugacity coefficient \f$\mathrm{\phi^\kappa_\alpha}\f$ is connected
diff --git a/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh b/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh
index e7094f2b2fe082e0b74a274847f8eb063976ab6d..d82af8a1ce69611619596e54ef54e7358fb49bba 100644
--- a/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh
+++ b/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh
@@ -28,7 +28,7 @@
 #include <dumux/porousmediumflow/implicit/problem.hh>
 
 #include <dumux/porousmediumflow/3pwateroil/model.hh>
-#include <dumux/material/fluidsystems/h2oheavyoilfluidsystem.hh>
+#include <dumux/material/fluidsystems/h2oheavyoil.hh>
 #include "3pwateroilsagdspatialparams.hh"
 
 namespace Dumux