diff --git a/dumux/boxmodels/1p/1pproperties.hh b/dumux/boxmodels/1p/1pproperties.hh
index 6040cf40f53dbf63d40262cb7db4de99f34ea166..ada38886cf652d5b2a60bb9313a42ba26eb08cce 100644
--- a/dumux/boxmodels/1p/1pproperties.hh
+++ b/dumux/boxmodels/1p/1pproperties.hh
@@ -53,13 +53,11 @@ NEW_TYPE_TAG(BoxOneP, INHERITS_FROM(BoxModel));
 NEW_PROP_TAG(NumPhases);   //!< Number of fluid phases in the system
 NEW_PROP_TAG(OnePIndices); //!< Enumerations for the 1p models
 NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters object
-NEW_PROP_TAG(Fluid); //!< The fluid for the single-phase problems
+NEW_PROP_TAG(FluidSystem); //!< The type of the fluid system to use
+NEW_PROP_TAG(Fluid); //!< The fluid used for the default fluid system
 NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
 NEW_PROP_TAG(UpwindWeight); //!< Returns weight of the upwind cell when calculating fluxes
 
-// use the CG solver preconditioned by ILU-0
-SET_TYPE_PROP(BoxOneP, LinearSolver, Dumux::BoxCGILU0Solver<TypeTag> );
-
 // \}
 };
 
diff --git a/dumux/boxmodels/1p/1ppropertydefaults.hh b/dumux/boxmodels/1p/1ppropertydefaults.hh
index 106cc14691ec88e41ed4cd34b782034dc6fa5047..2781aa9d9c0a88c017af466c1cc29bb831b1dfe6 100644
--- a/dumux/boxmodels/1p/1ppropertydefaults.hh
+++ b/dumux/boxmodels/1p/1ppropertydefaults.hh
@@ -37,6 +37,12 @@
 #include "1pfluxvariables.hh"
 #include "1pindices.hh"
 
+#include <dumux/material/fluidsystems/gasphase.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/components/nullcomponent.hh>
+
+#include <dumux/material/MpNcfluidsystems/1pfluidsystem.hh>
+
 namespace Dumux
 {
 // \{
@@ -69,6 +75,23 @@ SET_TYPE_PROP(BoxOneP, OnePIndices, OnePIndices);
 //! fluxes. Use central differences by default.
 SET_SCALAR_PROP(BoxOneP, UpwindWeight, 0.5);
 
+//! The fluid system to use by default
+SET_PROP(BoxOneP, FluidSystem)
+{ private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Fluid)) Fluid;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+    
+public:
+    typedef OnePFluidSystem<Scalar, Fluid> type;
+};
+
+SET_PROP(BoxOneP, Fluid)
+{ private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+    typedef Dumux::LiquidPhase<Scalar, Dumux::NullComponent<Scalar> > type;
+};
+
 // \}
 };
 
diff --git a/dumux/boxmodels/1p/1pvolumevariables.hh b/dumux/boxmodels/1p/1pvolumevariables.hh
index 8cb85d474be617493e17acab49b6185c2b2d4bdf..356d29d5dc9c6b5fa0c8c7fe80b8534c724ef5b2 100644
--- a/dumux/boxmodels/1p/1pvolumevariables.hh
+++ b/dumux/boxmodels/1p/1pvolumevariables.hh
@@ -29,6 +29,8 @@
 #include "1pproperties.hh"
 #include <dumux/boxmodels/common/boxvolumevariables.hh>
 
+#include <dumux/material/MpNcfluidstates/immisciblefluidstate.hh>
+
 namespace Dumux
 {
 
@@ -42,13 +44,13 @@ template <class TypeTag>
 class OnePVolumeVariables : public BoxVolumeVariables<TypeTag>
 {
     typedef BoxVolumeVariables<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) Implementation;
 
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
 
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
-    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) Implementation;
 
     typedef typename GridView::template Codim<0>::Entity Element;
 
@@ -59,12 +61,15 @@ class OnePVolumeVariables : public BoxVolumeVariables<TypeTag>
 
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePIndices)) Indices;
-    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Fluid)) Fluid;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
 
     typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
     typedef Dune::FieldVector<Scalar, dim> LocalPosition;
 
 public:
+    //! Type of the fluid state
+    typedef Dumux::ImmiscibleFluidState<Scalar, FluidSystem> FluidState;
+
     /*!
      * \brief Update all quantities for a given control volume.
      *
@@ -84,20 +89,35 @@ public:
     {
         ParentType::update(priVars, problem, element, elemGeom, scvIdx, isOldSol);
 
-        this->asImp_().updateTemperature_(priVars,
-                                          element,
-                                          elemGeom,
-                                          scvIdx,
-                                          problem);
+        // set the phase temperature and pressure
+        asImp_().updateTemperature_(priVars,
+                                    element,
+                                    elemGeom,
+                                    scvIdx,
+                                    problem);
+        fluidState_.setPressure(/*phaseIdx=*/0, priVars[Indices::pressureIdx]);
 
-        pressure_ = priVars[Indices::pressureIdx];
-        density_ = Fluid::density(temperature(), pressure());
-        viscosity_ = Fluid::viscosity(temperature(), pressure());
+        // saturation in a single phase is always 1 and thus redundant
+        // to set. But since we use the fluid state shared by the
+        // immiscible multi-phase models, so we have to set it here...
+        fluidState_.setSaturation(/*phaseIdx=*/0, 1.0);
+        
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updatePhase(fluidState_, /*phaseIdx=*/0);
+
+        Scalar value = FluidSystem::density(fluidState_, paramCache,  /*phaseIdx=*/0);
+        fluidState_.setDensity(/*phaseIdx=*/0, value);
+
+        value = FluidSystem::viscosity(fluidState_, paramCache,  /*phaseIdx=*/0);
+        fluidState_.setViscosity(/*phaseIdx=*/0, value);
 
         // porosity
         porosity_ = problem.spatialParameters().porosity(element,
                                                          elemGeom,
                                                          scvIdx);
+
+        // energy related quantities
+        asImp_().updateEnergy_(paramCache, priVars, problem, element, elemGeom, scvIdx, isOldSol);
     };
 
     /*!
@@ -108,28 +128,28 @@ public:
      * identical.
      */
     Scalar temperature() const
-    { return temperature_; }
+    { return fluidState_.temperature(); }
 
     /*!
      * \brief Returns the effective pressure of a given phase within
      *        the control volume.
      */
     Scalar pressure() const
-    { return pressure_; }
+    { return fluidState_.pressure(/*phaseIdx=*/0); }
 
     /*!
      * \brief Returns the mass density of a given phase within the
      *        control volume.
      */
     Scalar density() const
-    { return density_; }
+    { return fluidState_.density(/*phaseIdx=*/0); }
 
     /*!
      * \brief Returns the dynamic viscosity of the fluid within the
      *        control volume.
      */
     Scalar viscosity() const
-    { return viscosity_; }
+    { return fluidState_.viscosity(/*phaseIdx=*/0); }
 
     /*!
      * \brief Returns the average porosity within the control volume.
@@ -137,19 +157,42 @@ public:
     Scalar porosity() const
     { return porosity_; }
 
+    /*!
+     * \brief Returns the fluid state of the control volume.
+     */
+    const FluidState &fluidState() const
+    { return fluidState_; }
+    
 protected:
     void updateTemperature_(const PrimaryVariables &priVars,
                             const Element &element,
                             const FVElementGeometry &elemGeom,
                             int scvIdx,
                             const Problem &problem)
-    { temperature_ = problem.boxTemperature(element, elemGeom, scvIdx); }
+    { fluidState_.setTemperature(problem.boxTemperature(element, elemGeom, scvIdx)); }
 
-    Scalar temperature_;
-    Scalar pressure_;
-    Scalar density_;
-    Scalar viscosity_;
+    /*!
+     * \brief Called by update() to compute the energy related quantities
+     */
+    template <class ParameterCache>
+    void updateEnergy_(ParameterCache &paramCache,
+                       const PrimaryVariables &sol,
+                       const Problem &problem,
+                       const Element &element,
+                       const FVElementGeometry &elemGeom,
+                       int vertIdx,
+                       bool isOldSol)
+    { }
+
+    FluidState fluidState_;
     Scalar porosity_;
+
+private:
+    Implementation &asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    const Implementation &asImp_() const
+    { return *static_cast<const Implementation*>(this); }
 };
 
 }
diff --git a/dumux/boxmodels/2p/2pvolumevariables.hh b/dumux/boxmodels/2p/2pvolumevariables.hh
index 97f9216239e1ca8bfb09d7b3d718ab0ff27e6fa4..50ec3933830f4414ebaf0124b70fdd4760fde41e 100644
--- a/dumux/boxmodels/2p/2pvolumevariables.hh
+++ b/dumux/boxmodels/2p/2pvolumevariables.hh
@@ -102,12 +102,12 @@ public:
                            scvIdx,
                            isOldSol);
 
-        asImp().updateTemperature_(priVars,
-                                   element,
-                                   elemGeom,
-                                   scvIdx,
-                                   problem);
-
+        asImp_().updateTemperature_(priVars,
+                                    element,
+                                    elemGeom,
+                                    scvIdx,
+                                    problem);
+        
         // material law parameters
         const MaterialLawParams &materialParams =
             problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
@@ -126,6 +126,16 @@ public:
         porosity_ = problem.spatialParameters().porosity(element,
                                                          elemGeom,
                                                          scvIdx);
+
+
+#warning "TODO: it would be nice to avoid creating the parameter cache a second time here," \
+    " but that requires to abandon/modify model.completeFluidState()."  \
+    " (alternatively a completeFluidState() method can be introduced in the 2pni model.)"
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updateAll(fluidState_);
+
+        // energy related quantities
+        asImp_().updateEnergy_(paramCache, priVars, problem, element, elemGeom, scvIdx, isOldSol);
     }
 
     /*!
@@ -202,15 +212,28 @@ protected:
         fluidState_.setTemperature(problem.boxTemperature(element, elemGeom, scvIdx));
     }
 
+    /*!
+     * \brief Called by update() to compute the energy related quantities
+     */
+    template <class ParameterCache>
+    void updateEnergy_(ParameterCache &paramCache,
+                       const PrimaryVariables &sol,
+                       const Problem &problem,
+                       const Element &element,
+                       const FVElementGeometry &elemGeom,
+                       int vertIdx,
+                       bool isOldSol)
+    { }
+
     FluidState fluidState_;
     Scalar porosity_;
     Scalar mobility_[numPhases];
 
 private:
-    Implementation &asImp()
+    Implementation &asImp_()
     { return *static_cast<Implementation*>(this); }
 
-    const Implementation &asImp() const
+    const Implementation &asImp_() const
     { return *static_cast<const Implementation*>(this); }
 };
 
diff --git a/dumux/boxmodels/2p2c/2p2cvolumevariables.hh b/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
index 01beda45f64de865abdcc995badb2eaef097ecf6..6ee0c2461c1d7c4f5a8a3ef4558ba956cf5c75d7 100644
--- a/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cvolumevariables.hh
@@ -266,8 +266,7 @@ protected:
                        const FVElementGeometry &elemGeom,
                        int vertIdx,
                        bool isOldSol)
-    {
-    }
+    { }
 
     Scalar porosity_;        //!< Effective porosity within the control volume
     Scalar relativePermeability_[numPhases];  //!< Relative permeability within the control volume
diff --git a/dumux/boxmodels/2pni/2pnivolumevariables.hh b/dumux/boxmodels/2pni/2pnivolumevariables.hh
index f5d636b67f32fd7269b9b87ec7f73c555139cc8f..eb04350fab0f6a96b5feb9543d89ba6383080d1d 100644
--- a/dumux/boxmodels/2pni/2pnivolumevariables.hh
+++ b/dumux/boxmodels/2pni/2pnivolumevariables.hh
@@ -64,76 +64,6 @@ class TwoPNIVolumeVariables : public TwoPVolumeVariables<TypeTag>
     typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
 
 public:
-    /*!
-     * \brief Update all quantities for a given control volume.
-     *
-     * \param priVars The local primary variable vector
-     * \param problem The problem object
-     * \param element The current element
-     * \param elemGeom The finite-volume geometry in the box scheme
-     * \param scvIdx The local index of the SCV (sub-control volume)
-     * \param isOldSol Evaluate function with solution of current or previous time step
-     *
-     */
-    void update(const PrimaryVariables &priVars,
-                const Problem &problem,
-                const Element &element,
-                const FVElementGeometry &elemGeom,
-                int scvIdx,
-                bool isOldSol)
-    {
-        typedef Indices I;
-
-        // vertex update data for the mass balance
-        ParentType::update(priVars,
-                           problem,
-                           element,
-                           elemGeom,
-                           scvIdx,
-                           isOldSol);
-
-        typename FluidSystem::ParameterCache paramCache;
-        // TODO: this calculates the cached parameters a second time
-        // and is thus inefficient!
-        paramCache.updateAll(this->fluidState());
-
-        // the internal energies
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            internalEnergy_[phaseIdx] =
-                FluidSystem::internalEnergy(this->fluidState(),
-                                            paramCache,
-                                            phaseIdx);
-        }
-        Valgrind::CheckDefined(internalEnergy_);
-    }
-
-    /*!
-        * \brief Update the temperature for a given control volume.
-        *
-        * \param priVars The local primary variable vector
-        * \param element The current element
-        * \param elemGeom The finite-volume geometry in the box scheme
-        * \param scvIdx The local index of the SCV (sub-control volume)
-        * \param problem The problem object
-        *
-        */
-
-    // this method gets called by the parent class
-    void updateTemperature_(const PrimaryVariables &priVars,
-                            const Element &element,
-                            const FVElementGeometry &elemGeom,
-                            int scvIdx,
-                            const Problem &problem)
-    {
-        // retrieve temperature from primary variables
-        this->fluidState_.setTemperature(priVars[temperatureIdx]);
-
-        heatCapacity_ =
-            problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
-
-        Valgrind::CheckDefined(heatCapacity_);
-    }
-
     /*!
      * \brief Returns the total internal energy of a phase in the
      *        sub-control volume.
@@ -142,7 +72,7 @@ public:
      *
      */
     Scalar internalEnergy(int phaseIdx) const
-    { return internalEnergy_[phaseIdx]; };
+    { return this->fluidState_.internalEnergy(phaseIdx); };
 
     /*!
      * \brief Returns the total enthalpy of a phase in the sub-control
@@ -151,10 +81,7 @@ public:
      *  \param phaseIdx The phase index
      */
     Scalar enthalpy(int phaseIdx) const
-    {
-        return internalEnergy_[phaseIdx]
-            + this->fluidState().pressure(phaseIdx)
-            / this->fluidState().density(phaseIdx); };
+    { return this->fluidState_.enthalpy(phaseIdx); };
 
     /*!
      * \brief Returns the total heat capacity \f$\mathrm{[J/K*m^3]}\f$ of the rock matrix in
@@ -164,7 +91,54 @@ public:
     { return heatCapacity_; };
 
 protected:
-    Scalar internalEnergy_[numPhases];
+    // this method gets called by the parent class. since this method
+    // is protected, we are friends with our parent..
+    friend class TwoPVolumeVariables<TypeTag>;
+
+    /*!
+     * \brief Update the temperature for a given control volume.
+     *
+     * \param priVars The local primary variable vector
+     * \param element The current element
+     * \param elemGeom The finite-volume geometry in the box scheme
+     * \param scvIdx The local index of the SCV (sub-control volume)
+     * \param problem The problem object
+     *
+     */
+    void updateTemperature_(const PrimaryVariables &priVars,
+                            const Element &element,
+                            const FVElementGeometry &elemGeom,
+                            int scvIdx,
+                            const Problem &problem)
+    {
+        // retrieve temperature from primary variables
+        this->fluidState_.setTemperature(priVars[temperatureIdx]);
+    }
+    
+    /*!
+     * \brief Called by update() to compute the energy related quantities
+     */
+    template <class ParameterCache>
+    void updateEnergy_(ParameterCache &paramCache,
+                       const PrimaryVariables &sol,
+                       const Problem &problem,
+                       const Element &element,
+                       const FVElementGeometry &elemGeom,
+                       int vertIdx,
+                       bool isOldSol)
+    {
+        // copmute and set the internal energies of the fluid phases
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            Scalar u = FluidSystem::internalEnergy(this->fluidState_, paramCache, phaseIdx);
+
+            this->fluidState_.setInternalEnergy(phaseIdx, u);
+        }
+
+        // copmute and set the heat capacity of the solid phase
+        heatCapacity_ = problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
+        Valgrind::CheckDefined(heatCapacity_);
+    }
+
     Scalar heatCapacity_;
 };
 
diff --git a/dumux/material/MpNcfluidsystems/1pfluidsystem.hh b/dumux/material/MpNcfluidsystems/1pfluidsystem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..bac653903f6bb7ba21c57855234ff290b2642c44
--- /dev/null
+++ b/dumux/material/MpNcfluidsystems/1pfluidsystem.hh
@@ -0,0 +1,371 @@
+/*****************************************************************************
+ *   Copyright (C) 2011 by Andreas Lauser                                    *
+ *   Institute of Hydraulic Engineering                                      *
+ *   University of Stuttgart, Germany                                        *
+ *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
+ *                                                                           *
+ *   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 for single phase models.
+ */
+#ifndef DUMUX_1P_FLUIDSYSTEM_HH
+#define DUMUX_1P_FLUIDSYSTEM_HH
+
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/fluidsystems/gasphase.hh>
+
+#include <dune/common/exceptions.hh>
+
+#include "nullparametercache.hh"
+
+namespace Dumux {
+/*!
+ * \ingroup Fluidsystems
+ *
+ * \brief A fluid system for single phase models.
+ *
+ * The fluid is defined as a template parameter. For existing
+ * components the Dumux::LiquidPhase<Component> and
+ * Dumux::GasPhase<Component> may be used.
+ */
+template <class Scalar, class Fluid>
+class OnePFluidSystem
+{
+    typedef OnePFluidSystem<Scalar, Fluid> ThisType;
+
+public:
+    //! The type of parameter cache objects
+    typedef Dumux::NullParameterCache ParameterCache;
+
+    /****************************************
+     * Fluid phase related static parameters
+     ****************************************/
+
+    //! Number of phases in the fluid system
+    static constexpr int numPhases = 1;
+
+    /*!
+     * \brief Return the human readable name of a fluid phase
+     *
+     * \param phaseIdx The index of the fluid phase to consider
+     */
+    static const char *phaseName(int phaseIdx)
+    {
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+        return Fluid::name();
+    }
+
+    /*!
+     * \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 Fluid::isLiquid();
+    }
+
+    /*!
+     * \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 only a single component is involved. 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 immisibility
+        return true;
+    }
+
+    /****************************************
+     * Component related static parameters
+     ****************************************/
+
+    //! Number of components in the fluid system
+    static constexpr int numComponents = 1;
+
+    /*!
+     * \brief Return the human readable name of a component
+     *
+     * \param compIdx The index of the component to consider
+     */
+    static const char *componentName(int compIdx)
+    {
+        assert(0 <= compIdx && compIdx < numComponents);
+
+        return Fluid::name();
+    }
+
+    /*!
+     * \brief Return the molar mass of a component in [kg/mol].
+     *
+     * \param compIdx index of the component
+     */
+    static Scalar molarMass(int compIdx)
+    {
+        assert(0 <= compIdx && compIdx < numComponents);
+
+        return Fluid::molarMass();
+    }
+
+    /*!
+     * \brief Critical temperature of a component [K].
+     *
+     * \param compIdx The index of the component to consider
+     */
+    static Scalar criticalTemperature(int compIdx)
+    {
+        assert(0 <= compIdx && compIdx < numComponents);
+
+        return Fluid::criticalTemperature();
+    };
+
+    /*!
+     * \brief Critical pressure of a component [Pa].
+     *
+     * \param compIdx The index of the component to consider
+     */
+    static Scalar criticalPressure(int compIdx)
+    {
+        assert(0 <= compIdx && compIdx < numComponents);
+
+        return Fluid::criticalPressure();
+    };
+
+    /*!
+     * \brief The acentric factor of a component [].
+     *
+     * \param compIdx The index of the component to consider
+     */
+    static Scalar acentricFactor(int compIdx)
+    {
+        assert(0 <= compIdx && compIdx < numComponents);
+
+        return Fluid::acentricFactor();
+    };
+
+    /****************************************
+     * thermodynamic relations
+     ****************************************/
+
+    /*!
+     * \brief Initialize the fluid system's static parameters
+     */
+    static void init()
+    { }
+
+    /*!
+     * \brief Calculate the density [kg/m^3] of a fluid phase
+     *
+     * \param fluidState An abitrary fluid state
+     * \param paramCache The fluid system's parameter cache
+     * \param phaseIdx The index of the fluid phase to consider
+     */
+    template <class FluidState>
+    static Scalar density(const FluidState &fluidState,
+                          const ParameterCache &paramCache,
+                          int phaseIdx)
+    {
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+        Scalar temperature = fluidState.temperature(phaseIdx);
+        Scalar pressure = fluidState.pressure(phaseIdx);
+        return Fluid::density(temperature, pressure);
+    }
+
+    /*!
+     * \brief Calculate the dynamic viscosity of a fluid phase [Pa*s]
+     *
+     * \param fluidState An abitrary fluid state
+     * \param paramCache The fluid system's parameter cache
+     * \param phaseIdx The index of the fluid phase to consider
+     */
+    template <class FluidState>
+    static Scalar viscosity(const FluidState &fluidState,
+                            const ParameterCache &paramCache,
+                            int phaseIdx)
+    {
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+        Scalar temperature = fluidState.temperature(phaseIdx);
+        Scalar pressure = fluidState.pressure(phaseIdx);
+        return Fluid::viscosity(temperature, pressure);
+    }
+
+    /*!
+     * \brief Calculate the fugacity coefficient [Pa] of an individual
+     *        component in a fluid phase
+     *
+     * The fugacity coefficient \f$\phi_\kappa\f$ is connected to the
+     * fugacity \f$f_\kappa\f$ and the component's molarity
+     * \f$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,
+                                      int phaseIdx,
+                                      int compIdx)
+    {
+        assert(0 <= phaseIdx  && phaseIdx < numPhases);
+        assert(0 <= compIdx  && compIdx < numComponents);
+
+        if (phaseIdx == compIdx)
+            // TODO (?): calculate the real fugacity coefficient of
+            // the component in the fluid. Probably that's not worth
+            // the effort, since the fugacity coefficient of the other
+            // component is infinite anyway...
+            return 1.0;
+        return std::numeric_limits<Scalar>::infinity();
+    }
+
+    /*!
+     * \brief Calculate the molecular diffusion coefficient for a
+     *        component in a fluid phase [mol^2 * s / (kg*m^3)]
+     *
+     * Molecular diffusion of a compoent \f$\kappa\f$ is caused by a
+     * gradient of the chemical potential and follows the law
+     *
+     * \f[ J = - D \grad mu_\kappa \f]
+     *
+     * where \f$\mu_\kappa\f$ is the component's chemical potential,
+     * \f$D\f$ is the diffusion coefficient and \f$J\f$ is the
+     * diffusive flux. \f$mu_\kappa\f$ is connected to the component's
+     * fugacity \f$f_\kappa\f$ by the relation
+     *
+     * \f[ \mu_\kappa = R T_\alpha \mathrm{ln} \frac{f_\kappa}{p_\alpha} \f]
+     *
+     * where \f$p_\alpha\f$ and \f$T_\alpha\f$ are the fluid phase'
+     * pressure and temperature.
+     *
+     * \param fluidState An abitrary fluid state
+     * \param paramCache The fluid system's parameter cache
+     * \param phaseIdx The index of the fluid phase to consider
+     * \param compIdx The index of the component to consider
+     */
+    template <class FluidState>
+    static Scalar diffusionCoefficient(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx,
+                                       int compIdx)
+    {
+        DUNE_THROW(Dune::InvalidStateException, "Not applicable: Diffusion coefficients");
+    };
+    
+    /*!
+     * \brief Given a phase's composition, temperature and pressure,
+     *        return the binary diffusion coefficient for components
+     *        \f$i\f$ and \f$j\f$ in this phase.
+     *
+     * \param fluidState An abitrary fluid state
+     * \param paramCache The fluid system's parameter cache
+     * \param phaseIdx The index of the fluid phase to consider
+     * \param compIIdx The index of the first component to consider
+     * \param compJIdx The index of the second component to consider
+     */
+    template <class FluidState>
+    static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
+                                             const ParameterCache &paramCache,
+                                             int phaseIdx,
+                                             int compIIdx,
+                                             int compJIdx)
+
+    {
+        DUNE_THROW(Dune::InvalidStateException, "Not applicable: Binary diffusion coefficients");
+    };
+ 
+    /*!
+     * \brief Given a phase's composition, temperature, pressure and
+     *        density, calculate its specific internal energy [J/kg].
+     *
+     * \param fluidState An abitrary fluid state
+     * \param paramCache The fluid system's parameter cache
+     * \param phaseIdx The index of the fluid phase to consider
+     */
+    template <class FluidState>
+    static Scalar internalEnergy(const FluidState &fluidState,
+                                 const ParameterCache &paramCache,
+                                 int phaseIdx)
+    {
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+        Scalar temperature = fluidState.temperature(phaseIdx);
+        Scalar pressure = fluidState.pressure(phaseIdx);
+        return Fluid::internalEnergy(temperature, pressure);
+    }
+
+    /*!
+     * \brief Thermal conductivity of a fluid phase [W/(m^2 K/m)].
+     *
+     * Use the conductivity of air and water as a first approximation.
+     * Source:
+     * http://en.wikipedia.org/wiki/List_of_thermal_conductivities
+     *
+     * \param fluidState An abitrary fluid state
+     * \param paramCache The fluid system's parameter cache
+     * \param phaseIdx The index of the fluid phase to consider
+     */
+    template <class FluidState>
+    static Scalar thermalConductivity(const FluidState &fluidState,
+                                      const ParameterCache &paramCache,
+                                      int phaseIdx)
+    {
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+        Scalar temperature = fluidState.temperature(phaseIdx);
+        Scalar pressure = fluidState.pressure(phaseIdx);
+        return Fluid::thermalConductivity(temperature, pressure);
+    }
+
+    /*!
+     * \brief Specific isobaric heat capacity of a fluid phase.
+     *        \f$\mathrm{[J/kg]}\f$.
+     *
+     * \param fluidState An abitrary fluid state
+     * \param paramCache The fluid system's parameter cache
+     * \param phaseIdx The index of the fluid phase to consider
+     */
+    template <class FluidState>
+    static Scalar heatCapacity(const FluidState &fluidState,
+                               const ParameterCache &paramCache,
+                               int phaseIdx)
+    {
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+        Scalar temperature = fluidState.temperature(phaseIdx);
+        Scalar pressure = fluidState.pressure(phaseIdx);
+        return Fluid::heatCapacity(temperature, pressure);
+    }
+    
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh b/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh
index 85b5ddfc0850f7eff8495471d4714e3d4cee9cc7..dc36667a279578a7a320ea98791de16b798437c0 100644
--- a/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh
+++ b/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh
@@ -234,6 +234,24 @@ public:
         return NonWettingPhase::density(temperature, pressure);
     }
 
+    /*!
+     * \brief Return the viscosity of a phase [Pa*s].
+     *
+     */
+    using Base::viscosity;
+    template <class FluidState>
+    static Scalar viscosity(const FluidState &fluidState,
+                            int phaseIdx)
+    {
+        assert(0 <= phaseIdx  && phaseIdx < numPhases);
+
+        Scalar temperature = fluidState.temperature(phaseIdx);
+        Scalar pressure = fluidState.pressure(phaseIdx);
+        if (phaseIdx == wPhaseIdx)
+            return WettingPhase::viscosity(temperature, pressure);
+        return NonWettingPhase::viscosity(temperature, pressure);
+    }
+
     /*!
      * \brief Calculate the fugacity coefficient [Pa] of an individual
      *        component in a fluid phase
@@ -262,30 +280,6 @@ public:
         return std::numeric_limits<Scalar>::infinity();
     }
 
-    /*!
-     * \brief Return the viscosity of a phase [Pa*s].
-     *
-     * \param phaseIdx index of the phase
-     * \param temperature phase temperature in [K]
-     * \param pressure phase pressure in [Pa]
-     * \param fluidState The fluid state of the two-phase model
-     * \tparam FluidState the fluid state class of the two-phase model
-     * \return returns the viscosity of the phase [Pa*s]
-     */
-    using Base::viscosity;
-    template <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
-                            int phaseIdx)
-    {
-        assert(0 <= phaseIdx  && phaseIdx < numPhases);
-
-        Scalar temperature = fluidState.temperature(phaseIdx);
-        Scalar pressure = fluidState.pressure(phaseIdx);
-        if (phaseIdx == wPhaseIdx)
-            return WettingPhase::viscosity(temperature, pressure);
-        return NonWettingPhase::viscosity(temperature, pressure);
-    }
-
     /*!
      * \brief Calculate the binary molecular diffusion coefficient for
      *        a component in a fluid phase [mol^2 * s / (kg*m^3)]
diff --git a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
index fcece55f1e2d2407aff3c4b9895689f2b2156067..f643546ed46699cf8783c25d43e97ccb4c7febc0 100644
--- a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
+++ b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
@@ -279,7 +279,7 @@ public:
     }
 
     /*!
-     * \brief Calculate the molar volume [m^3/mol] of a fluid phase
+     * \brief Calculate the density [kg/m^3] of a fluid phase
      *
      * \param fluidState An abitrary fluid state
      * \param paramCache The fluid system's parameter cache
@@ -306,6 +306,31 @@ public:
             * fluidState.averageMolarMass(gPhaseIdx);
     };
 
+    /*!
+     * \brief Calculate the dynamic viscosity of a fluid phase [Pa*s]
+     *
+     * \param fluidState An abitrary fluid state
+     * \param paramCache The fluid system's parameter cache
+     * \param phaseIdx The index of the fluid phase to consider
+     */
+    template <class FluidState>
+    static Scalar viscosity(const FluidState &fluidState,
+                            const ParameterCache &paramCache,
+                            int phaseIdx)
+    {
+        assert(0 <= phaseIdx  && phaseIdx < numPhases);
+
+        Scalar T = fluidState.temperature(phaseIdx);
+        Scalar p = fluidState.pressure(phaseIdx);
+        if (phaseIdx == lPhaseIdx) {
+            // assume pure water for the liquid phase
+            return H2O::liquidViscosity(T, p);
+        }
+        
+        // assume pure nitrogen for the gas phase
+        return N2::gasViscosity(T, p);
+    };
+
     /*!
      * \brief Calculate the fugacity coefficient [Pa] of an individual
      *        component in a fluid phase
@@ -342,30 +367,6 @@ public:
         return 1.0; // ideal gas
     }
 
-    /*!
-     * \brief Calculate the dynamic viscosity of a fluid phase [Pa*s]
-     *
-     * \param fluidState An abitrary fluid state
-     * \param paramCache The fluid system's parameter cache
-     * \param phaseIdx The index of the fluid phase to consider
-     */
-    template <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
-                            const ParameterCache &paramCache,
-                            int phaseIdx)
-    {
-        assert(0 <= phaseIdx  && phaseIdx < numPhases);
-
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
-        if (phaseIdx == lPhaseIdx) {
-            // assume pure water for the liquid phase
-            return H2O::liquidViscosity(T, p);
-        }
-        
-        // assume pure nitrogen for the gas phase
-        return N2::gasViscosity(T, p);
-    };
 
     /*!
      * \brief Calculate the molecular diffusion coefficient for a
@@ -455,7 +456,7 @@ public:
 
     /*!
      * \brief Given a phase's composition, temperature, pressure and
-     *        density, calculate its specific enthalpy [J/kg].
+     *        density, calculate its specific internal energy [J/kg].
      *
      *  \todo This fluid system neglects the contribution of
      *        gas-molecules in the liquid phase. This contribution is
diff --git a/test/boxmodels/1p/1ptestproblem.hh b/test/boxmodels/1p/1ptestproblem.hh
index 9e67cf863c51550e2c8b182bc7f276b2e59522c0..2a5aad85cd8b873e81b53b8e9d0177391823e206 100644
--- a/test/boxmodels/1p/1ptestproblem.hh
+++ b/test/boxmodels/1p/1ptestproblem.hh
@@ -78,7 +78,8 @@ SET_PROP(OnePTestProblem, SpatialParameters)
 
 // Linear solver settings
 SET_TYPE_PROP(OnePTestProblem, LinearSolver, Dumux::BoxCGILU0Solver<TypeTag> );
-SET_INT_PROP(OnePTestProblem, LinearSolverVerbosity, 1);
+SET_INT_PROP(OnePTestProblem, LinearSolverVerbosity, 0);
+SET_SCALAR_PROP(OnePTestProblem, LinearSolverResidualReduction, 1e-12);
 SET_INT_PROP(OnePTestProblem, PreconditionerIterations, 1);
 SET_SCALAR_PROP(OnePTestProblem, PreconditionerRelaxation, 1.0);
 
diff --git a/test/boxmodels/1p/grids/test_1p_2d.dgf b/test/boxmodels/1p/grids/test_1p_2d.dgf
index 536c2946ca6804f87ac0be91f3e29a27ee531966..74b5ecc050539884116616c14edcd28407f5ec29 100644
--- a/test/boxmodels/1p/grids/test_1p_2d.dgf
+++ b/test/boxmodels/1p/grids/test_1p_2d.dgf
@@ -2,7 +2,7 @@ DGF
 Interval
 0 0   % first corner 
 1 1   % second corner
-10 10 % cells in x and y direction 
+100 100 % cells in x and y direction 
 #
 
 BOUNDARYDOMAIN