diff --git a/dumux/boxmodels/MpNc/MpNcproperties.hh b/dumux/boxmodels/MpNc/MpNcproperties.hh
index 7cc4bf98cb7c05fef48375443aac362a3c32e419..7050f9c65aff5e5cb67e53aca839bfcf9bb12755 100644
--- a/dumux/boxmodels/MpNc/MpNcproperties.hh
+++ b/dumux/boxmodels/MpNc/MpNcproperties.hh
@@ -67,7 +67,7 @@ NEW_PROP_TAG(MPNCVtkAddVarPressures);
 NEW_PROP_TAG(MPNCVtkAddVelocities);
 NEW_PROP_TAG(MPNCVtkAddDensities);
 NEW_PROP_TAG(MPNCVtkAddMobilities);
-NEW_PROP_TAG(MPNCVtkAddMeanMolarMass);
+NEW_PROP_TAG(MPNCVtkAddAverageMolarMass);
 NEW_PROP_TAG(MPNCVtkAddMassFractions);
 NEW_PROP_TAG(MPNCVtkAddMoleFractions);
 NEW_PROP_TAG(MPNCVtkAddMolarities);
diff --git a/dumux/boxmodels/MpNc/MpNcpropertydefaults.hh b/dumux/boxmodels/MpNc/MpNcpropertydefaults.hh
index b4006423ecd4af474ce250d2b8f0946c11288a85..e67ab9defb1dd595153049c1de154ea8d5d6d8e1 100644
--- a/dumux/boxmodels/MpNc/MpNcpropertydefaults.hh
+++ b/dumux/boxmodels/MpNc/MpNcpropertydefaults.hh
@@ -214,7 +214,7 @@ SET_BOOL_PROP(BoxMPNC, MPNCVtkAddVarPressures, false);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddVelocities, false);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddDensities, true);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMobilities, true);
-SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMeanMolarMass, false);
+SET_BOOL_PROP(BoxMPNC, MPNCVtkAddAverageMolarMass, false);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMassFractions, false);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMoleFractions, true);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMolarities, false);
diff --git a/dumux/boxmodels/MpNc/MpNcvolumevariables.hh b/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
index 38253ce229df2f787a16d8e3d00698d02c3f1197..55e846d944c5985823be671c72448c10f74919d0 100644
--- a/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
+++ b/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
@@ -90,7 +90,6 @@ class MPNCVolumeVariables
 
     typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
 
-
 public:
     //! The return type of the fluidState() method
     typedef typename MassVolumeVariables::FluidState FluidState;
@@ -124,8 +123,9 @@ public:
                            elemGeom,
                            scvIdx,
                            isOldSol);
+        ParentType::checkDefined();
 
-        typename FluidSystem::MutableParameters mutParams;
+        typename FluidSystem::ParameterCache paramCache;
 
         /////////////
         // set the phase saturations
@@ -133,24 +133,21 @@ public:
         Scalar sumSat = 0;
         for (int i = 0; i < numPhases - 1; ++i) {
             sumSat += priVars[S0Idx + i];
-            mutParams.setSaturation(i, priVars[S0Idx + i]);
+            fluidState_.setSaturation(i, priVars[S0Idx + i]);
         }
         Valgrind::CheckDefined(sumSat);
-        mutParams.setSaturation(numPhases - 1, 1.0 - sumSat);
+        fluidState_.setSaturation(numPhases - 1, 1.0 - sumSat);
 
         /////////////
         // set the fluid phase temperatures
         /////////////
-        // update the temperature part of the energy module
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx){
-            Scalar T = EnergyVolumeVariables::getTemperature(priVars,
-                                                             element,
-                                                             elemGeom,
-                                                             scvIdx,
-                                                             problem,
-                                                             phaseIdx);
-            mutParams.setTemperature(phaseIdx, T);
-        }
+        EnergyVolumeVariables::updateTemperatures(fluidState_,
+                                                  paramCache,
+                                                  priVars,
+                                                  element,
+                                                  elemGeom,
+                                                  scvIdx,
+                                                  problem);
 
         /////////////
         // set the phase pressures
@@ -161,16 +158,17 @@ public:
             problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
         // capillary pressures
         Scalar capPress[numPhases];
-        MaterialLaw::capillaryPressures(capPress, materialParams, mutParams);
+        MaterialLaw::capillaryPressures(capPress, materialParams, fluidState_);
         // add to the pressure of the first fluid phase
         Scalar p0 = priVars[p0Idx];
         for (int i = 0; i < numPhases; ++ i)
-            mutParams.setPressure(i, p0 - capPress[0] + capPress[i]);
+            fluidState_.setPressure(i, p0 - capPress[0] + capPress[i]);
 
         /////////////
         // set the fluid compositions
         /////////////
-        MassVolumeVariables::update(mutParams,
+        MassVolumeVariables::update(fluidState_,
+                                    paramCache,
                                     priVars,
                                     hint_,
                                     problem,
@@ -178,7 +176,6 @@ public:
                                     elemGeom,
                                     scvIdx);
         MassVolumeVariables::checkDefined();
-        ParentType::checkDefined();
 
         /////////////
         // Porosity
@@ -189,7 +186,6 @@ public:
                                                          elemGeom,
                                                          scvIdx);
         Valgrind::CheckDefined(porosity_);
-        ParentType::checkDefined();
 
         /////////////
         // Phase mobilities
@@ -198,31 +194,30 @@ public:
         // relative permeabilities
         MaterialLaw::relativePermeabilities(relativePermeability_,
                                             materialParams, 
-                                            mutParams);
+                                            fluidState_);
 
         // dynamic viscosities
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             // viscosities
-            Scalar mu = FluidSystem::computeViscosity(mutParams, phaseIdx);
-            mutParams.setViscosity(phaseIdx, mu);
+            Scalar mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
+            fluidState_.setViscosity(phaseIdx, mu);
         }
-        ParentType::checkDefined();
 
         /////////////
         // diffusion
         /////////////
 
         // update the diffusion part of the volume data
-        DiffusionVolumeVariables::update(mutParams, *this, problem);
+        DiffusionVolumeVariables::update(fluidState_, paramCache, *this, problem);
         DiffusionVolumeVariables::checkDefined();
-        ParentType::checkDefined();
 
         /////////////
         // energy
         /////////////
 
         // update the remaining parts of the energy module
-        EnergyVolumeVariables::update(mutParams,
+        EnergyVolumeVariables::update(fluidState_,
+                                      paramCache,
                                       priVars,
                                       element,
                                       elemGeom,
@@ -230,18 +225,16 @@ public:
                                       problem);
 
         EnergyVolumeVariables::checkDefined();
-        ParentType::checkDefined();
-
 
-        // assign the fluid state to be stored (needed in the next update)
-        fluidState_.assign(mutParams);
+        // make sure the quantities in the fluid state are well-defined
         fluidState_.checkDefined();
 
         // specific interfacial area,
         // well also all the dimensionless numbers :-)
         // well, also the mass transfer rate
         IAVolumeVariables::update(*this,
-                                  mutParams,
+                                  fluidState_,
+                                  paramCache,
                                   priVars,
                                   problem,
                                   element,
@@ -249,11 +242,7 @@ public:
                                   scvIdx);
         IAVolumeVariables::checkDefined();
 
-        ParentType::checkDefined();
-
         checkDefined();
-
-        ParentType::checkDefined();
     }
 
     /*!
@@ -330,7 +319,7 @@ public:
         // difference of sum of mole fractions in the phase from 100%
         Scalar a = 1;
         for (int i = 0; i < numComponents; ++i)
-            a -= fluidState.moleFrac(phaseIdx, i);
+            a -= fluidState.moleFraction(phaseIdx, i);
         return a;
     }
 
diff --git a/dumux/boxmodels/MpNc/MpNcvolumevariablesia.hh b/dumux/boxmodels/MpNc/MpNcvolumevariablesia.hh
index 23764d1a10bdda3b77795a63e1621a24bf30ec59..ffc2b3e2106b77bc401f0a0057c944381c368719 100644
--- a/dumux/boxmodels/MpNc/MpNcvolumevariablesia.hh
+++ b/dumux/boxmodels/MpNc/MpNcvolumevariablesia.hh
@@ -61,9 +61,10 @@ public:
     /*!
      * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
      */
-    template <class MutableParams>
+    template <class FluidState, class ParameterCache>
     void update(const VolumeVariables & volVars,
-                const MutableParams & mutParams,
+                const FluidState &fluidState,
+                const ParameterCache &paramCache,
                 const PrimaryVariables &priVars,
                 const Problem &problem,
                 const Element & element,
diff --git a/dumux/boxmodels/MpNc/MpNcvtkwritercommon.hh b/dumux/boxmodels/MpNc/MpNcvtkwritercommon.hh
index a6b53eedae6fe34ea7ce93d8a18ba29f99ea79a8..9d751c5a67d1b4635e8c8635cee8791d066da378 100644
--- a/dumux/boxmodels/MpNc/MpNcvtkwritercommon.hh
+++ b/dumux/boxmodels/MpNc/MpNcvtkwritercommon.hh
@@ -85,7 +85,7 @@ public:
         velocityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddVelocities);
         densityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddDensities);
         mobilityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMobilities);
-        meanMolarMassOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMeanMolarMass);
+        averageMolarMassOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddAverageMolarMass);
         massFracOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMassFractions);
         moleFracOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMoleFractions);
         molarityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMolarities);
@@ -116,7 +116,7 @@ public:
         if (pressureOutput_) this->resizePhaseBuffer_(pressure_);
         if (densityOutput_) this->resizePhaseBuffer_(density_);
         if (mobilityOutput_) this->resizePhaseBuffer_(mobility_);
-        if (meanMolarMassOutput_) this->resizePhaseBuffer_(meanMolarMass_);
+        if (averageMolarMassOutput_) this->resizePhaseBuffer_(averageMolarMass_);
 
         if (moleFracOutput_) this->resizePhaseComponentBuffer_(moleFrac_);
         if (massFracOutput_) this->resizePhaseComponentBuffer_(massFrac_);
@@ -154,10 +154,10 @@ public:
                 if (pressureOutput_) pressure_[phaseIdx][I] = volVars.fluidState().pressure(phaseIdx);
                 if (densityOutput_) density_[phaseIdx][I] = volVars.fluidState().density(phaseIdx);
                 if (mobilityOutput_) mobility_[phaseIdx][I] = volVars.mobility(phaseIdx);
-                if (meanMolarMassOutput_) meanMolarMass_[phaseIdx][I] = volVars.fluidState().meanMolarMass(phaseIdx);
+                if (averageMolarMassOutput_) averageMolarMass_[phaseIdx][I] = volVars.fluidState().averageMolarMass(phaseIdx);
                 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-                    if (moleFracOutput_) moleFrac_[phaseIdx][compIdx][I] = volVars.fluidState().moleFrac(phaseIdx, compIdx);
-                    if (massFracOutput_) massFrac_[phaseIdx][compIdx][I] = volVars.fluidState().massFrac(phaseIdx, compIdx);
+                    if (moleFracOutput_) moleFrac_[phaseIdx][compIdx][I] = volVars.fluidState().moleFraction(phaseIdx, compIdx);
+                    if (massFracOutput_) massFrac_[phaseIdx][compIdx][I] = volVars.fluidState().massFraction(phaseIdx, compIdx);
                     if (molarityOutput_) molarity_[phaseIdx][compIdx][I] = volVars.fluidState().molarity(phaseIdx, compIdx);
                 }
             }
@@ -213,8 +213,8 @@ public:
         if (densityOutput_)
             this->commitPhaseBuffer_(writer, "rho_%s", density_);
 
-        if (meanMolarMassOutput_)
-            this->commitPhaseBuffer_(writer, "M_%s", meanMolarMass_);
+        if (averageMolarMassOutput_)
+            this->commitPhaseBuffer_(writer, "M_%s", averageMolarMass_);
 
         if (mobilityOutput_)
             this->commitPhaseBuffer_(writer, "lambda_%s", mobility_);
@@ -263,13 +263,13 @@ private:
     bool massFracOutput_;
     bool moleFracOutput_;
     bool molarityOutput_;
-    bool meanMolarMassOutput_;
+    bool averageMolarMassOutput_;
 
     PhaseBuffer saturation_;
     PhaseBuffer pressure_;
     PhaseBuffer density_;
     PhaseBuffer mobility_;
-    PhaseBuffer meanMolarMass_;
+    PhaseBuffer averageMolarMass_;
 
     PhaseVelocityField velocity_;
     ScalarBuffer boxSurface_;
diff --git a/dumux/boxmodels/MpNc/diffusion/volumevariables.hh b/dumux/boxmodels/MpNc/diffusion/volumevariables.hh
index 096724bc9b29d6eb56be001d4da01ae28388c907..b8be0213826152da406a96492f3d41b8f821aa52 100644
--- a/dumux/boxmodels/MpNc/diffusion/volumevariables.hh
+++ b/dumux/boxmodels/MpNc/diffusion/volumevariables.hh
@@ -129,18 +129,17 @@ class MPNCVolumeVariablesDiffusion<TypeTag, false>
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
-    typedef typename FluidSystem::MutableParameters MutableParameters;
 
 public:
     MPNCVolumeVariablesDiffusion()
     {}
 
-    void update(MutableParameters &mutParams,
-        const VolumeVariables &volVars,
+    template <class FluidState, class ParameterCache>
+    void update(FluidState &fluidState,
+                ParameterCache &paramCache,
+                const VolumeVariables &volVars,
                 const Problem &problem)
-    {
-    };
+    { };
 
     Scalar diffCoeffL(int compIdx) const
     { return 0; };
diff --git a/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh b/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh
index 8d0bba1ecf544607e5f4a92a826761a8c8038614..5da3e6ab9d6dfa0c696a3daa983333967d72be58 100644
--- a/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh
+++ b/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh
@@ -27,6 +27,8 @@
 #ifndef DUMUX_MPNC_ENERGY_VOLUME_VARIABLES_HH
 #define DUMUX_MPNC_ENERGY_VOLUME_VARIABLES_HH
 
+#include <dumux/material/MpNcfluidstates/equilibriumfluidstate.hh>
+
 namespace Dumux
 {
 /*!
@@ -57,19 +59,24 @@ class MPNCVolumeVariablesEnergy
     enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
     //typedef typename GET_PROP_TYPE(TypeTag, PTAG(MPNCEnergyIndices)) EnergyIndices;
 
-
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+    typedef typename FluidSystem::ParameterCache ParameterCache; 
+    typedef Dumux::EquilibriumFluidState<Scalar, FluidSystem> FluidState;
+    
 public:
     /*!
      * \brief Update the temperature of the sub-control volume.
      */
-    Scalar getTemperature(const PrimaryVariables &sol,
-                          const Element &element,
-                          const FVElementGeometry &elemGeom,
-                          int scvIdx,
-                          const Problem &problem,
-                          const int temperatureIdx) const
+    void updateTemperatures(FluidState &fs,
+                            ParameterCache &paramCache,
+                            const PrimaryVariables &sol,
+                            const Element &element,
+                            const FVElementGeometry &elemGeom,
+                            int scvIdx,
+                            const Problem &problem) const
     {
-        return problem.boxTemperature(element, elemGeom, scvIdx);
+        Scalar T = problem.boxTemperature(element, elemGeom, scvIdx);
+        fs.setTemperature(T);
     }
 
 
@@ -79,39 +86,23 @@ public:
      *
      * Since we are isothermal, we don't need to do anything!
      */
-    template <class MutableParams>
-    void update(MutableParams &mutParams,
+    void update(FluidState &fs,
+                ParameterCache &paramCache,
                 const PrimaryVariables &sol,
                 const Element &element,
                 const FVElementGeometry &elemGeom,
                 int scvIdx,
                 const Problem &problem)
     {
-        temperature_ = problem.boxTemperature(element, elemGeom, scvIdx);
-
-        // set the fluid temperatures
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-            mutParams.setTemperature(phaseIdx, temperature_);
     }
 
-    /*!
-     * \brief Return the temperature of any given phase [K]
-     */
-    Scalar temperature(int phaseIdx = 0) const
-    { return temperature_; }
-
     /*!
      * \brief If running under valgrind this produces an error message
      *        if some of the object's attributes is undefined.
      */
     void checkDefined() const
     {
-        Valgrind::CheckDefined(temperature_);
     }
-
-protected:
-    Scalar temperature_;
-
 };
 
 /*!
@@ -130,7 +121,6 @@ class MPNCVolumeVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*kineticEnergyT
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
 
-    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(MPNCIndices)) Indices;
 
     enum { numPhases        = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
@@ -139,28 +129,33 @@ class MPNCVolumeVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*kineticEnergyT
     enum { numEnergyEqs     = Indices::NumPrimaryEnergyVars};
     enum { temperature0Idx = Indices::temperatureIdx };
 
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+    typedef typename FluidSystem::ParameterCache ParameterCache; 
+    typedef Dumux::EquilibriumFluidState<Scalar, FluidSystem> FluidState;
+    
 public:
     /*!
      * \brief Update the temperature of the sub-control volume.
      */
-    Scalar getTemperature(const PrimaryVariables &sol,
-                          const Element &element,
-                          const FVElementGeometry &elemGeom,
-                          int scvIdx,
-                          const Problem &problem,
-                          const int dummy) const
+    void updateTemperatures(FluidState &fs,
+                            ParameterCache &paramCache,
+                            const PrimaryVariables &sol,
+                            const Element &element,
+                            const FVElementGeometry &elemGeom,
+                            int scvIdx,
+                            const Problem &problem) const
     {
         // retrieve temperature from solution vector
-        return sol[temperatureIdx];
+        Scalar T = sol[temperatureIdx];
+        fs.setTemperature(T);
     }
 
     /*!
      * \brief Update the enthalpy and the internal energy for a given
      *        control volume.
      */
-    template <class MutableParams>
-    void update(MutableParams &mutParams,
-                const PrimaryVariables &sol,
+    void update(FluidState &fs,
+                ParameterCache &paramCache,
                 const Element &element,
                 const FVElementGeometry &elemGeom,
                 int scvIdx,
@@ -168,24 +163,20 @@ public:
     {
         Valgrind::SetUndefined(*this);
 
-        // set the fluid temperatures
-        temperature_ = sol[temperature0Idx];
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-            mutParams.setTemperature(phaseIdx, temperature_);
-
+        // heat capacities of the fluids plus the porous medium
         heatCapacity_ =
             problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
         Valgrind::CheckDefined(heatCapacity_);
 
         soilDensity_ =
-                problem.spatialParameters().soilDensity(element, elemGeom, scvIdx);
+            problem.spatialParameters().soilDensity(element, elemGeom, scvIdx);
         Valgrind::CheckDefined(soilDensity_);
 
         // set the enthalpies
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            Scalar h =
-                FluidSystem::computeEnthalpy(mutParams, phaseIdx);
-            mutParams.setEnthalpy(phaseIdx, h);
+            Scalar u =
+                FluidSystem::internalEnergy(fs, paramCache, phaseIdx);
+            fs.setInternalEnergy(phaseIdx, u);
         }
     }
 
@@ -196,13 +187,6 @@ public:
     Scalar heatCapacity() const
     { return heatCapacity_; };
 
-    /*!
-     * \brief Returns the temperature in fluid / solid phase(s)
-     *        the sub-control volume.
-     */
-    Scalar temperature(int phaseIdx = 0) const
-    { return temperature_; }
-
     /*!
      * \brief Returns the total density of the given soil [kg / m^3] in
      *        the sub-control volume.
@@ -218,13 +202,11 @@ public:
     {
         Valgrind::CheckDefined(heatCapacity_);
         Valgrind::CheckDefined(soilDensity_);
-        Valgrind::CheckDefined(temperature_);
     };
 
 protected:
     Scalar heatCapacity_;
     Scalar soilDensity_;
-    Scalar temperature_;
 };
 
 } // end namepace
diff --git a/dumux/boxmodels/MpNc/mass/MpNcvolumevariablesmass.hh b/dumux/boxmodels/MpNc/mass/MpNcvolumevariablesmass.hh
index 3887239dbfa4507219992b06f783cf57ae743822..d67f6574e1aa6e11ccb60c10163e21e1b678104e 100644
--- a/dumux/boxmodels/MpNc/mass/MpNcvolumevariablesmass.hh
+++ b/dumux/boxmodels/MpNc/mass/MpNcvolumevariablesmass.hh
@@ -60,7 +60,8 @@ class MPNCVolumeVariablesMass
     enum { fug0Idx = Indices::fug0Idx };
 
     typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
-
+    
+    typedef typename FluidSystem::ParameterCache ParameterCache;
 public:
     /*!
      * \brief The fluid state which is used by the volume variables to
@@ -75,11 +76,10 @@ public:
      * \brief Update composition of all phases in the mutable
      *        parameters from the primary variables.
      */
-    template <class MutableParams>
-    void update(MutableParams &mutParams,
+    void update(FluidState &fs,
+                ParameterCache &paramCache,
                 const PrimaryVariables &priVars,
                 const VolumeVariables *hint,
-
                 const Problem &problem,
                 const Element &element,
                 const FVElementGeometry &elemGeom,
@@ -97,10 +97,10 @@ public:
                 Scalar x_ij = 1.0/numComponents;
                 if (hint)
                     // use the hint for the initial mole fraction!
-                    x_ij = hint->fluidState().moleFrac(phaseIdx, compIdx);
+                    x_ij = hint->fluidState().moleFraction(phaseIdx, compIdx);
 
                 // set initial guess of the component's mole fraction
-                mutParams.setMoleFrac(phaseIdx,
+                fs.setMoleFraction(phaseIdx,
                                       compIdx,
                                       x_ij);
 
@@ -111,18 +111,18 @@ public:
             if (!hint)
                 // if no hint was given, we ask the constraint solver
                 // to make an initial guess
-                CompositionFromFugacitiesSolver::guessInitial(mutParams, phaseIdx, fug);
-            CompositionFromFugacitiesSolver::solve(mutParams, phaseIdx, fug);
+                CompositionFromFugacitiesSolver::guessInitial(fs, paramCache, phaseIdx, fug);
+            CompositionFromFugacitiesSolver::solve(fs, paramCache, phaseIdx, fug);
 
             /*
               std::cout << "final composition: " << FluidSystem::phaseName(phaseIdx) << "="
-              << mutParams.moleFrac(phaseIdx, 0) << " "
-              << mutParams.moleFrac(phaseIdx, 1) << " "
-              << mutParams.moleFrac(phaseIdx, 2) << " "
-              << mutParams.moleFrac(phaseIdx, 3) << " "
-              << mutParams.moleFrac(phaseIdx, 4) << " "
-              << mutParams.moleFrac(phaseIdx, 5) << " "
-              << mutParams.moleFrac(phaseIdx, 6) << "\n";
+              << fs.moleFrac(phaseIdx, 0) << " "
+              << fs.moleFrac(phaseIdx, 1) << " "
+              << fs.moleFrac(phaseIdx, 2) << " "
+              << fs.moleFrac(phaseIdx, 3) << " "
+              << fs.moleFrac(phaseIdx, 4) << " "
+              << fs.moleFrac(phaseIdx, 5) << " "
+              << fs.moleFrac(phaseIdx, 6) << "\n";
             */
 
         }
diff --git a/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh b/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
index a3fde6d5c87a8c5098a788c00fc8b563a3dd9245..1b40b3548454fcc7ccb0a95806d9185075c109db 100644
--- a/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
+++ b/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
@@ -38,18 +38,22 @@ namespace Dumux {
 template <class Scalar, class FluidSystem>
 class CompositionFromFugacities
 {
-    typedef typename FluidSystem::MutableParameters MutableParameters;
-
     enum { numPhases = FluidSystem::numPhases };
     enum { numComponents = FluidSystem::numComponents };
-    
+
+    typedef typename FluidSystem::ParameterCache ParameterCache;
+
 public:
     typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
 
     /*!
      * \brief Guess an initial value for the composition of the phase.
      */
-    static void guessInitial(MutableParameters &mutParams, int phaseIdx, const ComponentVector &fugVec)
+    template <class FluidState>
+    static void guessInitial(FluidState &fluidState,
+                             ParameterCache &paramCache,
+                             int phaseIdx,
+                             const ComponentVector &fugVec)
     {
         if (FluidSystem::isIdealMixture(phaseIdx)) 
             return;
@@ -57,9 +61,9 @@ public:
         // Pure component fugacities
         for (int i = 0; i < numComponents; ++ i) {
             //std::cout << f << " -> " << mutParams.fugacity(phaseIdx, i)/f << "\n";
-            mutParams.setMoleFrac(phaseIdx,
-                                  i, 
-                                  1.0/numComponents);
+            fluidState.setMoleFraction(phaseIdx,
+                                   i, 
+                                   1.0/numComponents);
         }
     };
 
@@ -69,14 +73,16 @@ public:
      *
      * The phase's fugacities must already be set.
      */
-    static void solve(MutableParameters &mutParams, 
+    template <class FluidState>
+    static void solve(FluidState &fluidState,
+                      ParameterCache &paramCache,
                       int phaseIdx, 
                       const ComponentVector &targetFug)
     {
         // use a much more efficient method in case the phase is an
         // ideal mixture
         if (FluidSystem::isIdealMixture(phaseIdx)) {
-            solveIdealMix_(mutParams, phaseIdx, targetFug);
+            solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
             return;
         }
 
@@ -85,7 +91,7 @@ public:
         // save initial composition in case something goes wrong
         Dune::FieldVector<Scalar, numComponents> xInit;
         for (int i = 0; i < numComponents; ++i) {
-            xInit[i] = mutParams.moleFrac(phaseIdx, i);
+            xInit[i] = fluidState.moleFraction(phaseIdx, i);
         };
         
         /////////////////////////
@@ -99,22 +105,25 @@ public:
         // right hand side
         Dune::FieldVector<Scalar, numComponents> b;
         
+        fluidState.updateAverageMolarMass(phaseIdx);
+        paramCache.updatePhase(fluidState, phaseIdx);
+        
         // maximum number of iterations
         const int nMax = 25;
         for (int nIdx = 0; nIdx < nMax; ++nIdx) {
             // calculate Jacobian matrix and right hand side
-            linearize_(J, b, mutParams, phaseIdx, targetFug); 
+            linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug); 
             Valgrind::CheckDefined(J);
             Valgrind::CheckDefined(b);
 
             /*
             std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
             for (int i = 0; i < FluidSystem::numComponents; ++i)
-                std::cout << mutParams.moleFrac(phaseIdx, i) << " ";
+                std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
             std::cout << "\n";
             std::cout << FluidSystem::phaseName(phaseIdx) << "Phase phi: ";
             for (int i = 0; i < FluidSystem::numComponents; ++i)
-                std::cout << mutParams.fugacityCoeff(phaseIdx, i) << " ";
+                std::cout << fluidState.fugacityCoeff(phaseIdx, i) << " ";
             std::cout << "\n";
             */
 
@@ -131,10 +140,10 @@ public:
             /*
             std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
             for (int i = 0; i < FluidSystem::numComponents; ++i)
-                std::cout << mutParams.moleFrac(phaseIdx, i) << " ";
+                std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
             std::cout << "\n";
             std::cout << "J: " << J << "\n";
-            std::cout << "rho: " << mutParams.density(phaseIdx) << "\n";
+            std::cout << "rho: " << fluidState.density(phaseIdx) << "\n";
             std::cout << "delta: " << x << "\n";
             std::cout << "defect: " << b << "\n";
 
@@ -145,14 +154,12 @@ public:
 
             // update the fluid composition. b is also used to store
             // the defect for the next iteration.
-            Scalar relError = update_(mutParams, x, b, phaseIdx, targetFug);
+            Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
             //std::cout << "relError: " << relError << "\n";
             
             if (relError < 1e-9) {
-                mutParams.updateMeanMolarMass(phaseIdx);
-                
-                Scalar Vm = FluidSystem::computeMolarVolume(mutParams, phaseIdx);
-                mutParams.setMolarVolume(phaseIdx, Vm);
+                Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
+                fluidState.setDensity(phaseIdx, rho);
 
                 //std::cout << "num iterations: " << nIdx << "\n";
                 return;
@@ -163,8 +170,8 @@ public:
                    "Calculating the " << FluidSystem::phaseName(phaseIdx)
                    << "Phase composition failed. Initial {x} = {"
                    << xInit 
-                   << "}, {fug_t} = {" << targetFug << "}, p = " << mutParams.pressure(phaseIdx)
-                   << ", T = " << mutParams.temperature(phaseIdx));
+                   << "}, {fug_t} = {" << targetFug << "}, p = " << fluidState.pressure(phaseIdx)
+                   << ", T = " << fluidState.temperature(phaseIdx));
     };
 
 
@@ -172,26 +179,35 @@ protected:
     // update the phase composition in case the phase is an ideal
     // mixture, i.e. the component's fugacity coefficients are
     // independent of the phase's composition.
-    static void solveIdealMix_(MutableParameters &mutParams, 
+    template <class FluidState>
+    static void solveIdealMix_(FluidState &fluidState, 
+                               ParameterCache &paramCache,
                                int phaseIdx,
                                const ComponentVector &fugacities)
     {        
         for (int i = 0; i < numComponents; ++ i) {
-            Scalar phi = FluidSystem::computeFugacityCoeff(mutParams, phaseIdx, i);
-            Scalar gamma = phi * mutParams.pressure(phaseIdx);
-            mutParams.setFugacityCoeff(phaseIdx, i, phi);
-            mutParams.setMoleFrac(phaseIdx, i, fugacities[i]/gamma);
+            Scalar phi = FluidSystem::fugacityCoefficient(fluidState,
+                                                          paramCache,
+                                                          phaseIdx,
+                                                          i);
+            Scalar gamma = phi * fluidState.pressure(phaseIdx);
+            fluidState.setFugacityCoefficient(phaseIdx, i, phi);
+            fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
         };
-        mutParams.updateMeanMolarMass(phaseIdx);
         
-        Scalar Vm = FluidSystem::computeMolarVolume(mutParams, phaseIdx);
-        mutParams.setMolarVolume(phaseIdx, Vm);
+        fluidState.updateAverageMolarMass(phaseIdx);
+        paramCache.updatePhase(fluidState, phaseIdx);
+        
+        Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
+        fluidState.setDensity(phaseIdx, rho);
         return;
     };
 
+    template <class FluidState>
     static void linearize_(Dune::FieldMatrix<Scalar, numComponents, numComponents> &J,
                            Dune::FieldVector<Scalar, numComponents> &defect,
-                           MutableParameters &mutParams, 
+                           FluidState &fluidState, 
+                           ParameterCache &paramCache,
                            int phaseIdx,
                            const ComponentVector &targetFug)
     {       
@@ -201,11 +217,12 @@ protected:
         // calculate the defect (deviation of the current fugacities
         // from the target fugacities)
         for (int i = 0; i < numComponents; ++ i) {
-            Scalar phi = FluidSystem::computeFugacityCoeff(mutParams,
-                                                           phaseIdx,
-                                                           i);
-            Scalar f = phi*mutParams.pressure(phaseIdx)*mutParams.moleFrac(phaseIdx, i);
-            mutParams.setFugacityCoeff(phaseIdx, i, phi);
+            Scalar phi = FluidSystem::fugacityCoefficient(fluidState,
+                                                          paramCache,
+                                                          phaseIdx,
+                                                          i);
+            Scalar f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
+            fluidState.setFugacityCoefficient(phaseIdx, i, phi);
                 
             defect[i] = targetFug[i] - f;
         }
@@ -221,21 +238,24 @@ protected:
             // forward differences
 
             // deviate the mole fraction of the i-th component
-            Scalar x_i = mutParams.moleFrac(phaseIdx, i);
-            mutParams.setMoleFrac(phaseIdx, i, x_i + eps);
+            Scalar x_i = fluidState.moleFraction(phaseIdx, i);
+            fluidState.setMoleFraction(phaseIdx, i, x_i + eps);
+            fluidState.updateAverageMolarMass(phaseIdx);
+            paramCache.updatePhaseSingleMoleFraction(fluidState, phaseIdx, i);
 
             // compute new defect and derivative for all component
             // fugacities
             for (int j = 0; j < numComponents; ++j) {
                 // compute the j-th component's fugacity coefficient ...
-                Scalar phi = FluidSystem::computeFugacityCoeff(mutParams,
-                                                               phaseIdx,
-                                                               j);
+                Scalar phi = FluidSystem::fugacityCoefficient(fluidState,
+                                                              paramCache,
+                                                              phaseIdx,
+                                                              j);
                 // ... and its fugacity ...
                 Scalar f = 
                     phi *
-                    mutParams.pressure(phaseIdx) *
-                    mutParams.moleFrac(phaseIdx, j);
+                    fluidState.pressure(phaseIdx) *
+                    fluidState.moleFraction(phaseIdx, j);
                 // as well as the defect for this fugacity
                 Scalar defJPlusEps = targetFug[j] - f;
 
@@ -245,14 +265,18 @@ protected:
             }
             
             // reset composition to original value
-            mutParams.setMoleFrac(phaseIdx, i, x_i);
-            
+            fluidState.setMoleFraction(phaseIdx, i, x_i);
+            fluidState.updateAverageMolarMass(phaseIdx);
+            paramCache.updatePhaseSingleMoleFraction(fluidState, phaseIdx, i);
+
             // end forward differences
             ////////
         }
     }
 
-    static Scalar update_(MutableParameters &mutParams,
+    template <class FluidState>
+    static Scalar update_(FluidState &fluidState,
+                          ParameterCache &paramCache,
                           Dune::FieldVector<Scalar, numComponents> &x,
                           Dune::FieldVector<Scalar, numComponents> &b,
                           int phaseIdx,
@@ -264,10 +288,10 @@ protected:
         Scalar sumDelta = 0;
         Scalar sumx = 0;
         for (int i = 0; i < numComponents; ++i) {
-            origComp[i] = mutParams.moleFrac(phaseIdx, i);
+            origComp[i] = fluidState.moleFraction(phaseIdx, i);
             relError = std::max(relError, std::abs(x[i]));
             
-            sumx += std::abs(mutParams.moleFrac(phaseIdx, i));
+            sumx += std::abs(fluidState.moleFraction(phaseIdx, i));
             sumDelta += std::abs(x[i]);
         };
 
@@ -278,7 +302,7 @@ protected:
             x /= (sumDelta/maxDelta);
 #endif
 
-        //Scalar curDefect = calculateDefect_(mutParams, phaseIdx, targetFug);
+        //Scalar curDefect = calculateDefect_(fluidState, phaseIdx, targetFug);
         //Scalar nextDefect;
         //Scalar sumMoleFrac = 0.0;
         //ComponentVector newB(1e100);
@@ -297,16 +321,19 @@ protected:
                 else
                     newx = 0;
 #endif
-                mutParams.setMoleFrac(phaseIdx, i, newx);
+                fluidState.setMoleFraction(phaseIdx, i, newx);
                 //sumMoleFrac += std::abs(newx);
             }
 
+            fluidState.updateAverageMolarMass(phaseIdx);
+            paramCache.updatePhaseComposition(fluidState, phaseIdx);
+
             /*
             // if the sum of the mole fractions gets 0, we take the
             // original composition divided by 100
             if (sumMoleFrac < 1e-10) {
                 for (int i = 0; i < numComponents; ++i) {
-                    mutParams.setMoleFrac(phaseIdx, i, origComp[i]/100);
+                    fluidState.setMoleFraction(phaseIdx, i, origComp[i]/100);
                 }
                 return relError;
             }
@@ -315,13 +342,13 @@ protected:
             /*
             // calculate new residual
             for (int i = 0; i < numComponents; ++i) {
-                Scalar phi = FluidSystem::computeFugacityCoeff(mutParams,
+                Scalar phi = FluidSystem::computeFugacityCoeff(fluidState,
                                                                phaseIdx,
                                                                i);
-                mutParams.setFugacityCoeff(phaseIdx, i, phi);
+                fluidState.setFugacityCoeff(phaseIdx, i, phi);
             }
 
-            nextDefect = calculateDefect_(mutParams, phaseIdx, targetFug);
+            nextDefect = calculateDefect_(fluidState, phaseIdx, targetFug);
             //std::cout << "try delta: " << x << "\n";
             //std::cout << "defect: old=" << curDefect << " new=" << nextDefect << "\n"; 
             if (nextDefect <= curDefect)
@@ -335,7 +362,8 @@ protected:
         return relError;
     }
 
-    static Scalar calculateDefect_(const MutableParameters &params, 
+    template <class FluidState>
+    static Scalar calculateDefect_(const FluidState &params, 
                                    int phaseIdx, 
                                    const ComponentVector &targetFug)
     {
diff --git a/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh b/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh
index 5f2e2d611c661cf89049dedfe524f225cfbd4ae2..d5f5eff514c247f44e2ff7d8088bc9fe1806c3ac 100644
--- a/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh
+++ b/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh
@@ -36,12 +36,12 @@ namespace Dumux
  *        multi-phase, multi-component fluid system assuming
  *        thermodynamic equilibrium.
  */
-template <class Scalar, class StaticParameters>
+template <class Scalar, class FluidSystem>
 class EquilibriumFluidState
 {
 public:
-    enum { numComponents = StaticParameters::numComponents };
-    enum { numPhases = StaticParameters::numPhases };
+    enum { numPhases = FluidSystem::numPhases };
+    enum { numComponents = FluidSystem::numComponents };
 
     EquilibriumFluidState()
     { Valgrind::SetUndefined(*this); }
@@ -63,18 +63,18 @@ public:
     /*!
      * \brief The mole fraction of a component in a phase []
      */
-    Scalar moleFrac(int phaseIdx, int compIdx) const
-    { return moleFrac_[phaseIdx][compIdx]; }
+    Scalar moleFraction(int phaseIdx, int compIdx) const
+    { return moleFraction_[phaseIdx][compIdx]; }
 
     /*!
      * \brief The mass fraction of a component in a phase []
      */
-    Scalar massFrac(int phaseIdx, int compIdx) const
+    Scalar massFraction(int phaseIdx, int compIdx) const
     { 
         return
-            sumMassFrac(phaseIdx)*
+            sumMassFractions(phaseIdx)*
             molarity(phaseIdx, compIdx)*
-            StaticParameters::molarMass(compIdx)
+            FluidSystem::molarMass(compIdx)
             /
             density(phaseIdx);
     }
@@ -84,22 +84,22 @@ public:
      *
      * We define this to be the same as the sum of all mass fractions.
      */
-    Scalar sumMoleFrac(int phaseIdx) const
-    { return sumMoleFrac_[phaseIdx]; }
+    Scalar sumMoleFractions(int phaseIdx) const
+    { return sumMoleFractions_[phaseIdx]; }
 
     /*!
      * \brief The sum of all component mass fractions in a phase []
      *
      * We define this to be the same as the sum of all mole fractions.
      */
-    Scalar sumMassFrac(int phaseIdx) const
-    { return sumMoleFrac_[phaseIdx]; }
+    Scalar sumMassFractions(int phaseIdx) const
+    { return sumMoleFractions_[phaseIdx]; }
 
     /*!
-     * \brief The mean molar mass of a fluid phase [kg/mol]
+     * \brief The average molar mass of a fluid phase [kg/mol]
      */
-    Scalar meanMolarMass(int phaseIdx) const
-    { return meanMolarMass_[phaseIdx]; }
+    Scalar averageMolarMass(int phaseIdx) const
+    { return averageMolarMass_[phaseIdx]; }
 
     /*!
      * \brief The concentration of a component in a phase [mol/m^3]
@@ -111,37 +111,37 @@ public:
      * http://en.wikipedia.org/wiki/Concentration
      */
     Scalar molarity(int phaseIdx, int compIdx) const
-    { return molarDensity(phaseIdx)*moleFrac(phaseIdx, compIdx); }
+    { return molarDensity(phaseIdx)*moleFraction(phaseIdx, compIdx); }
 
     /*!
      * \brief The fugacity of a component in a phase [Pa]
      */
     Scalar fugacity(int phaseIdx, int compIdx) const
-    { return fugacityCoeff(phaseIdx, compIdx)*moleFrac(phaseIdx, compIdx)*pressure(phaseIdx); }
+    { return fugacityCoefficient(phaseIdx, compIdx)*moleFraction(phaseIdx, compIdx)*pressure(phaseIdx); }
 
     /*!
      * \brief The fugacity coefficient of a component in a phase [Pa]
      */
-    Scalar fugacityCoeff(int phaseIdx, int compIdx) const
-    { return fugacityCoeff_[phaseIdx][compIdx]; }
+    Scalar fugacityCoefficient(int phaseIdx, int compIdx) const
+    { return fugacityCoefficient_[phaseIdx][compIdx]; }
 
     /*!
      * \brief The molar volume of a fluid phase [m^3/mol]
      */
     Scalar molarVolume(int phaseIdx) const
-    { return molarVolume_[phaseIdx]; }
+    { return 1/molarDensity(phaseIdx); }
 
     /*!
      * \brief The mass density of a fluid phase [kg/m^3]
      */
     Scalar density(int phaseIdx) const
-    { return molarDensity(phaseIdx)*meanMolarMass(phaseIdx); }
+    { return density_[phaseIdx]; }
 
     /*!
      * \brief The molar density of a fluid phase [mol/m^3]
      */
     Scalar molarDensity(int phaseIdx) const
-    { return 1/molarVolume(phaseIdx); }
+    { return density_[phaseIdx]/averageMolarMass(phaseIdx); }
 
     /*!
      * \brief The temperature of a fluid phase [K]
@@ -159,13 +159,13 @@ public:
      * \brief The specific enthalpy of a fluid phase [J/kg]
      */
     Scalar enthalpy(int phaseIdx) const
-    { return enthalpy_[phaseIdx]; }
+    { return internalEnergy_[phaseIdx] + pressure(phaseIdx)/(density(phaseIdx)); }
 
     /*!
      * \brief The specific internal energy of a fluid phase [J/kg]
      */
     Scalar internalEnergy(int phaseIdx) const
-    { return enthalpy(phaseIdx) - pressure(phaseIdx)/(density(phaseIdx)); }
+    { return internalEnergy_[phaseIdx]; }
 
     /*!
      * \brief The dynamic viscosity of a fluid phase [Pa s]
@@ -223,14 +223,14 @@ public:
     {
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-                moleFrac_[phaseIdx][compIdx] = fs.moleFrac(phaseIdx, compIdx);
-                fugacityCoeff_[phaseIdx][compIdx] = fs.fugacityCoeff(phaseIdx, compIdx);
+                moleFraction_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
+                fugacityCoefficient_[phaseIdx][compIdx] = fs.fugacityCoefficient(phaseIdx, compIdx);
             }
-            meanMolarMass_[phaseIdx] = fs.meanMolarMass(phaseIdx);
+            averageMolarMass_[phaseIdx] = fs.averageMolarMass(phaseIdx);
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
-            molarVolume_[phaseIdx] = fs.molarVolume(phaseIdx);
-            enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
+            density_[phaseIdx] = fs.density(phaseIdx);
+            internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
         }
         temperature_ = fs.temperature(0);
@@ -257,26 +257,26 @@ public:
     /*!
      * \brief Set the mole fraction of a component in a phase []
      */
-    void setMoleFrac(int phaseIdx, int compIdx, Scalar value)
-    { moleFrac_[phaseIdx][compIdx] = value; }   
+    void setMoleFraction(int phaseIdx, int compIdx, Scalar value)
+    { moleFraction_[phaseIdx][compIdx] = value; }   
 
     /*!
      * \brief Set the fugacity of a component in a phase []
      */
-    void setFugacityCoeff(int phaseIdx, int compIdx, Scalar value)
-    { fugacityCoeff_[phaseIdx][compIdx] = value; }   
+    void setFugacityCoefficient(int phaseIdx, int compIdx, Scalar value)
+    { fugacityCoefficient_[phaseIdx][compIdx] = value; }   
 
     /*!
-     * \brief Set the molar volume of a phase [m^3/mol]
+     * \brief Set the density of a phase [kg / m^3]
      */
-    void setMolarVolume(int phaseIdx, Scalar value)
-    { molarVolume_[phaseIdx] = value; }   
+    void setDensity(int phaseIdx, Scalar value)
+    { density_[phaseIdx] = value; }   
 
     /*!
      * \brief Set the specific enthalpy of a phase [J/m^3]
      */
-    void setEnthalpy(int phaseIdx, Scalar value)
-    { enthalpy_[phaseIdx] = value; }   
+    void setInternalEnergy(int phaseIdx, Scalar value)
+    { internalEnergy_[phaseIdx] = value; }   
 
     /*!
      * \brief Set the dynamic viscosity of a phase [Pa s]
@@ -288,18 +288,18 @@ public:
      * \brief Calculatate the mean molar mass of a phase given that
      *        all mole fractions have been set
      */
-    void updateMeanMolarMass(int phaseIdx)
+    void updateAverageMolarMass(int phaseIdx)
     {
-        meanMolarMass_[phaseIdx] = 0;
-        sumMoleFrac_[phaseIdx] = 0;
+        averageMolarMass_[phaseIdx] = 0;
+        sumMoleFractions_[phaseIdx] = 0;
      
         for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-            sumMoleFrac_[phaseIdx] += moleFrac_[phaseIdx][compIdx];
-            meanMolarMass_[phaseIdx] += moleFrac_[phaseIdx][compIdx]*StaticParameters::molarMass(compIdx);
+            sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx];
+            averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx);
         }
-        sumMoleFrac_[phaseIdx] = std::max(1e-15, sumMoleFrac_[phaseIdx]);
-        Valgrind::CheckDefined(meanMolarMass_[phaseIdx]);
-        Valgrind::CheckDefined(sumMoleFrac_[phaseIdx]);
+        sumMoleFractions_[phaseIdx] = std::max(1e-15, sumMoleFractions_[phaseIdx]);
+        Valgrind::CheckDefined(averageMolarMass_[phaseIdx]);
+        Valgrind::CheckDefined(sumMoleFractions_[phaseIdx]);
     }
     
     /*!
@@ -313,17 +313,17 @@ public:
     void checkDefined() const
     {
 #if HAVE_VALGRIND && ! defined NDEBUG
-        for (int i = 0; i < numPhases; ++i) {
-            for (int j = 0; j < numComponents; ++j) {
-                Valgrind::CheckDefined(moleFrac_[i][j]);
-                Valgrind::CheckDefined(fugacityCoeff_[i][j]);
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+                Valgrind::CheckDefined(moleFraction_[phaseIdx][compIdx]);
+                Valgrind::CheckDefined(fugacityCoefficient_[phaseIdx][compIdx]);
             }
-            Valgrind::CheckDefined(meanMolarMass_[i]);
-            Valgrind::CheckDefined(pressure_[i]);
-            Valgrind::CheckDefined(saturation_[i]);
-            Valgrind::CheckDefined(molarVolume_[i]);
-            //Valgrind::CheckDefined(enthalpy_[i]);
-            Valgrind::CheckDefined(viscosity_[i]);
+            Valgrind::CheckDefined(averageMolarMass_[phaseIdx]);
+            Valgrind::CheckDefined(pressure_[phaseIdx]);
+            Valgrind::CheckDefined(saturation_[phaseIdx]);
+            Valgrind::CheckDefined(density_[phaseIdx]);
+            //Valgrind::CheckDefined(enthalpy_[phaseIdx]);
+            Valgrind::CheckDefined(viscosity_[phaseIdx]);
         }
 
         Valgrind::CheckDefined(temperature_);
@@ -331,15 +331,15 @@ public:
     }
 
 protected:
-    Scalar moleFrac_[numPhases][numComponents];
-    Scalar fugacityCoeff_[numPhases][numComponents];
+    Scalar moleFraction_[numPhases][numComponents];
+    Scalar fugacityCoefficient_[numPhases][numComponents];
 
-    Scalar meanMolarMass_[numPhases];
-    Scalar sumMoleFrac_[numPhases];
+    Scalar averageMolarMass_[numPhases];
+    Scalar sumMoleFractions_[numPhases];
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
-    Scalar molarVolume_[numPhases];
-    Scalar enthalpy_[numPhases];
+    Scalar density_[numPhases];
+    Scalar internalEnergy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_;
 };
diff --git a/dumux/material/MpNcfluidstates/genericfluidstate.hh b/dumux/material/MpNcfluidstates/genericfluidstate.hh
index 16a00e2038635ebc10d8d2033e465db78e0974e0..d19bfce8dfb7e5ba11219d59f3a0e6720659f531 100644
--- a/dumux/material/MpNcfluidstates/genericfluidstate.hh
+++ b/dumux/material/MpNcfluidstates/genericfluidstate.hh
@@ -39,12 +39,12 @@ namespace Dumux
  *        multi-phase, multi-component fluid system without using
  *        any assumptions.
  */
-template <class Scalar, class StaticParameters>
+template <class Scalar, class FluidSystem>
 class GenericFluidState
 {
 public:
-    enum { numComponents    = StaticParameters::numComponents };
-    enum { numPhases        = StaticParameters::numPhases };
+    enum { numPhases = FluidSystem::numPhases };
+    enum { numComponents = FluidSystem::numComponents };
 
     GenericFluidState()
     { Valgrind::SetUndefined(*this); }
@@ -62,20 +62,20 @@ public:
     /*!
      * \brief The mole fraction of a component in a phase []
      */
-    Scalar moleFrac(int phaseIdx, int compIdx) const
-    { return moleFrac_[phaseIdx][compIdx]; }
+    Scalar moleFraction(int phaseIdx, int compIdx) const
+    { return moleFraction_[phaseIdx][compIdx]; }
 
 
     /*!
      * \brief The mass fraction of a component in a phase []
      */
-    Scalar massFrac(int phaseIdx, int compIdx) const
+    Scalar massFraction(int phaseIdx, int compIdx) const
     { 
         return
-            sumMassFrac(phaseIdx) *
-            moleFrac_[phaseIdx][compIdx] *
-            StaticParameters::molarMass(compIdx) /
-            meanMolarMass_[phaseIdx];
+            sumMassFraction(phaseIdx)
+            * moleFraction_[phaseIdx][compIdx]
+            * FluidSystem::molarMass(compIdx)
+            / averageMolarMass_[phaseIdx];
     }
 
     /*!
@@ -83,22 +83,22 @@ public:
      *
      * We define this to be the same as the sum of all mass fractions.
      */
-    Scalar sumMoleFrac(int phaseIdx) const
-    { return sumMoleFrac_[phaseIdx]; }
+    Scalar sumMoleFractions(int phaseIdx) const
+    { return sumMoleFractions_[phaseIdx]; }
 
     /*!
      * \brief The sum of all component mass fractions in a phase []
      *
      * We define this to be the same as the sum of all mole fractions.
      */
-    Scalar sumMassFrac(int phaseIdx) const
-    { return sumMoleFrac_[phaseIdx]; }
+    Scalar sumMassFraction(int phaseIdx) const
+    { return sumMoleFractions_[phaseIdx]; }
 
     /*!
      * \brief The mean molar mass of a fluid phase [kg/mol]
      */
-    Scalar meanMolarMass(int phaseIdx) const
-    { return meanMolarMass_[phaseIdx]; }
+    Scalar averageMolarMass(int phaseIdx) const
+    { return averageMolarMass_[phaseIdx]; }
 
     /*!
      * \brief The molar concentration of a component in a phase [mol/m^3]
@@ -109,40 +109,38 @@ public:
      * http://en.wikipedia.org/wiki/Concentration
      */
     Scalar molarity(int phaseIdx, int compIdx) const
-    { return molarDensity(phaseIdx)*moleFrac(phaseIdx, compIdx); }
+    { return molarDensity(phaseIdx)*moleFraction(phaseIdx, compIdx); }
 
     /*!
      * \brief The fugacity coefficient of a component in a phase []
      */
-    Scalar fugacityCoeff(int phaseIdx, int compIdx) const
-    { return fugacityCoeff_[phaseIdx][compIdx]; }
+    Scalar fugacityCoefficient(int phaseIdx, int compIdx) const
+    { return fugacityCoefficient_[phaseIdx][compIdx]; }
 
     /*!
      * \brief The fugacity of a component in a phase [Pa]
      */
     Scalar fugacity(int phaseIdx, int compIdx) const
-    { return pressure_[phaseIdx]*fugacityCoeff_[phaseIdx][compIdx]*moleFrac_[phaseIdx][compIdx]; }
+    { return pressure_[phaseIdx]*fugacityCoefficient_[phaseIdx][compIdx]*moleFraction_[phaseIdx][compIdx]; }
 
     /*!
      * \brief The molar volume of a fluid phase [m^3/mol]
      */
     Scalar molarVolume(int phaseIdx) const
-    { return molarVolume_[phaseIdx]; }
+    { return 1/molarDensity(phaseIdx); }
 
     /*!
      * \brief The mass density of a fluid phase [kg/m^3]
      */
     Scalar density(int phaseIdx) const
-    {
-        return 
-            meanMolarMass_[phaseIdx]/molarVolume_[phaseIdx];
-    }
+    { return density_[phaseIdx]; }
 
     /*!
      * \brief The molar density of a fluid phase [mol/m^3]
      */
     Scalar molarDensity(int phaseIdx) const
-    { return 1/molarVolume(phaseIdx); }
+    { return density_[phaseIdx]/averageMolarMass(phaseIdx); }
+
 
     /*!
      * \brief The temperature of a fluid phase [K]
@@ -157,16 +155,16 @@ public:
     { return pressure_[phaseIdx]; };
 
     /*!
-     * \brief The specific enthalpy of a fluid phase [J/kg]
+     * \brief The specific internal energy of a fluid phase [J/kg]
      */
-    Scalar enthalpy(int phaseIdx) const
-    { return enthalpy_[phaseIdx]; };
+    Scalar internalEnergy(int phaseIdx) const
+    { return internalEnergy_[phaseIdx]; };
 
     /*!
-     * \brief The specific internal energy of a fluid phase [J/kg]
+     * \brief The specific enthalpy of a fluid phase [J/kg]
      */
-    Scalar internalEnergy(int phaseIdx) const
-    { return enthalpy(phaseIdx) - pressure(phaseIdx)/(density(phaseIdx)); };
+    Scalar enthalpy(int phaseIdx) const
+    { return internalEnergy(phaseIdx) + pressure(phaseIdx)/density(phaseIdx); };
 
     /*!
      * \brief The dynamic viscosity of a fluid phase [Pa s]
@@ -200,28 +198,26 @@ public:
     /*!
      * \brief Set the mole fraction of a component in a phase []
      */
-    void setMoleFrac(int phaseIdx, int compIdx, Scalar value)
-    { moleFrac_[phaseIdx][compIdx] = value; }   
+    void setMoleFraction(int phaseIdx, int compIdx, Scalar value)
+    { moleFraction_[phaseIdx][compIdx] = value; }   
 
     /*!
      * \brief Set the fugacity of a component in a phase []
      */
-    void setFugacityCoeff(int phaseIdx, int compIdx, Scalar value)
-    { fugacityCoeff_[phaseIdx][compIdx] = value; }   
+    void setFugacityCoefficient(int phaseIdx, int compIdx, Scalar value)
+    { fugacityCoefficient_[phaseIdx][compIdx] = value; }   
 
     /*!
-     * \brief Set the molar volume of a phase [m^3/mol]
+     * \brief Set the density of a phase [kg / m^3]
      */
-    void setMolarVolume(int phaseIdx, Scalar value)
-    {
-        molarVolume_[phaseIdx] = value;
-    }
+    void setDensity(int phaseIdx, Scalar value)
+    { density_[phaseIdx] = value; }
 
     /*!
-     * \brief Set the specific enthalpy of a phase [J/m^3]
+     * \brief Set the specific internal energy of a phase [J/m^3]
      */
-    void setEnthalpy(int phaseIdx, Scalar value)
-    { enthalpy_[phaseIdx] = value; }
+    void setInternalEnergy(int phaseIdx, Scalar value)
+    { internalEnergy_[phaseIdx] = value; }
 
     /*!
      * \brief Set the dynamic viscosity of a phase [Pa s]
@@ -233,18 +229,18 @@ public:
      * \brief Calculatate the mean molar mass of a phase given that
      *        all mole fractions have been set
      */
-    void updateMeanMolarMass(int phaseIdx)
+    void updateAverageMolarMass(int phaseIdx)
     {
-        meanMolarMass_[phaseIdx] = 0;
-        sumMoleFrac_[phaseIdx] = 0;
+        averageMolarMass_[phaseIdx] = 0;
+        sumMoleFractions_[phaseIdx] = 0;
      
         for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-            sumMoleFrac_[phaseIdx] += moleFrac_[phaseIdx][compIdx];
-            meanMolarMass_[phaseIdx] += moleFrac_[phaseIdx][compIdx]*StaticParameters::molarMass(compIdx);
+            sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx];
+            averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx);
         }
-        sumMoleFrac_[phaseIdx] = std::max(1e-15, sumMoleFrac_[phaseIdx]);
-        Valgrind::CheckDefined(meanMolarMass_[phaseIdx]);
-        Valgrind::CheckDefined(sumMoleFrac_[phaseIdx]);
+        sumMoleFractions_[phaseIdx] = std::max(1e-15, sumMoleFractions_[phaseIdx]);
+        Valgrind::CheckDefined(averageMolarMass_[phaseIdx]);
+        Valgrind::CheckDefined(sumMoleFractions_[phaseIdx]);
     }
     
     /*!
@@ -256,15 +252,15 @@ public:
     {
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-                moleFrac_[phaseIdx][compIdx] = fs.moleFrac(phaseIdx, compIdx);
-                fugacityCoeff_[phaseIdx][compIdx] = fs.fugacityCoeff(phaseIdx, compIdx);
+                moleFraction_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
+                fugacityCoefficient_[phaseIdx][compIdx] = fs.fugacityCoefficient(phaseIdx, compIdx);
             }
-            updateMeanMolarMass(phaseIdx);
+            updateAverageMolarMass(phaseIdx);
 
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
-            molarVolume_[phaseIdx] = fs.molarVolume(phaseIdx);
-            enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
+            density_[phaseIdx] = fs.density(phaseIdx);
+            internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
             temperature_[phaseIdx] = fs.temperature(phaseIdx);
         }
@@ -283,30 +279,30 @@ public:
 #if HAVE_VALGRIND && ! defined NDEBUG
         for (int i = 0; i < numPhases; ++i) {
             for (int j = 0; j < numComponents; ++j) {
-                Valgrind::CheckDefined(fugacityCoeff_[i][j]);
-                Valgrind::CheckDefined(moleFrac_[i][j]);
+                Valgrind::CheckDefined(fugacityCoefficient_[i][j]);
+                Valgrind::CheckDefined(moleFraction_[i][j]);
             }
-            Valgrind::CheckDefined(meanMolarMass_[i]);
+            Valgrind::CheckDefined(averageMolarMass_[i]);
             Valgrind::CheckDefined(pressure_[i]);
             Valgrind::CheckDefined(saturation_[i]);
-            Valgrind::CheckDefined(molarVolume_[i]);
+            Valgrind::CheckDefined(density_[i]);
             Valgrind::CheckDefined(temperature_[i]);
-            //Valgrind::CheckDefined(enthalpy_[i]);
+            //Valgrind::CheckDefined(internalEnergy_[i]);
             Valgrind::CheckDefined(viscosity_[i]);
         }
 #endif // HAVE_VALGRIND
     }
 
 protected:
-    Scalar moleFrac_[numPhases][numComponents];
-    Scalar fugacityCoeff_[numPhases][numComponents];
+    Scalar moleFraction_[numPhases][numComponents];
+    Scalar fugacityCoefficient_[numPhases][numComponents];
 
-    Scalar meanMolarMass_[numPhases];
-    Scalar sumMoleFrac_[numPhases];
+    Scalar averageMolarMass_[numPhases];
+    Scalar sumMoleFractions_[numPhases];
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
-    Scalar molarVolume_[numPhases];
-    Scalar enthalpy_[numPhases];
+    Scalar density_[numPhases];
+    Scalar internalEnergy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_[numPhases];
 };
diff --git a/dumux/material/MpNcfluidstates/immisciblefluidstate.hh b/dumux/material/MpNcfluidstates/immisciblefluidstate.hh
index 586128df24c72aff68874bbd0854b9a73a9177fe..b3b548fcd35882a95ce6037bfe198bb4d9e3ba40 100644
--- a/dumux/material/MpNcfluidstates/immisciblefluidstate.hh
+++ b/dumux/material/MpNcfluidstates/immisciblefluidstate.hh
@@ -65,13 +65,13 @@ public:
     /*!
      * \brief The mole fraction of a component in a phase []
      */
-    Scalar moleFrac(int phaseIdx, int compIdx) const
+    Scalar moleFraction(int phaseIdx, int compIdx) const
     { return (phaseIdx == compIdx)?1.0:0.0; }
 
     /*!
      * \brief The mass fraction of a component in a phase []
      */
-    Scalar massFrac(int phaseIdx, int compIdx) const
+    Scalar massFraction(int phaseIdx, int compIdx) const
     { return (phaseIdx == compIdx)?1.0:0.0; }
 
     /*!
@@ -79,7 +79,7 @@ public:
      *
      * We define this to be the same as the sum of all mass fractions.
      */
-    Scalar sumMoleFrac(int phaseIdx) const
+    Scalar sumMoleFractions(int phaseIdx) const
     { return 1.0; }
 
     /*!
@@ -87,13 +87,13 @@ public:
      *
      * We define this to be the same as the sum of all mole fractions.
      */
-    Scalar sumMassFrac(int phaseIdx) const
+    Scalar sumMassFractions(int phaseIdx) const
     { return 1.0; }
 
     /*!
-     * \brief The mean molar mass of a fluid phase [kg/mol]
+     * \brief The average molar mass of a fluid phase [kg/mol]
      */
-    Scalar meanMolarMass(int phaseIdx) const
+    Scalar averageMolarMass(int phaseIdx) const
     { return StaticParameters::molarMass(/*compIdx=*/phaseIdx); }
 
     /*!
@@ -106,7 +106,7 @@ public:
      * http://en.wikipedia.org/wiki/Concentration
      */
     Scalar molarity(int phaseIdx, int compIdx) const
-    { return molarDensity(phaseIdx)*moleFrac(phaseIdx, compIdx); }
+    { return molarDensity(phaseIdx)*moleFraction(phaseIdx, compIdx); }
 
     /*!
      * \brief The fugacity of a component in a phase [Pa]
@@ -120,9 +120,9 @@ public:
     Scalar fugacity(int phaseIdx, int compIdx) const
     { 
         if (phaseIdx == compIdx)
-            return fugacityCoeff(phaseIdx, compIdx)*moleFrac(phaseIdx, compIdx)*pressure(phaseIdx);
+            return pressure(phaseIdx);
         else
-            return 0.0;
+            return 0;
     };
     
     /*!
@@ -133,10 +133,10 @@ public:
      * infinite. Beware that this will very likely break your code if
      * you don't keep that in mind.
      */
-    Scalar fugacityCoeff(int phaseIdx, int compIdx) const
+    Scalar fugacityCoefficient(int phaseIdx, int compIdx) const
     {
         if (phaseIdx == compIdx)
-            return fugacityCoeff_[phaseIdx];
+            return 1.0;
         else
             return std::numeric_limits<Scalar>::infinity();
     }
@@ -145,19 +145,19 @@ public:
      * \brief The molar volume of a fluid phase [m^3/mol]
      */
     Scalar molarVolume(int phaseIdx) const
-    { return molarVolume_[phaseIdx]; }
+    { return 1/molarDensity(phaseIdx); }
 
     /*!
      * \brief The mass density of a fluid phase [kg/m^3]
      */
     Scalar density(int phaseIdx) const
-    { return molarDensity(phaseIdx)*meanMolarMass(phaseIdx); }
+    { return density_[phaseIdx]; }
 
     /*!
      * \brief The molar density of a fluid phase [mol/m^3]
      */
     Scalar molarDensity(int phaseIdx) const
-    { return 1/molarVolume(phaseIdx); }
+    { return density_[phaseIdx]/averageMolarMass(phaseIdx); }
 
     /*!
      * \brief The temperature of a fluid phase [K]
@@ -175,13 +175,13 @@ public:
      * \brief The specific enthalpy of a fluid phase [J/kg]
      */
     Scalar enthalpy(int phaseIdx) const
-    { return enthalpy_[phaseIdx]; }
+    { return internalEnergy_[phaseIdx] + pressure(phaseIdx)/(density(phaseIdx)); }
 
     /*!
      * \brief The specific internal energy of a fluid phase [J/kg]
      */
     Scalar internalEnergy(int phaseIdx) const
-    { return enthalpy(phaseIdx) - pressure(phaseIdx)/density(phaseIdx); }
+    { return internalEnergy_[phaseIdx]; }
 
     /*!
      * \brief The dynamic viscosity of a fluid phase [Pa s]
@@ -228,11 +228,10 @@ public:
     void assign(const FluidState &fs)
     {
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            fugacityCoeff_[phaseIdx] = fs.fugacityCoeff(phaseIdx, phaseIdx);
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
-            molarVolume_[phaseIdx] = fs.molarVolume(phaseIdx);
-            enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
+            density_[phaseIdx] = fs.density(phaseIdx);
+            internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
         }
         temperature_ = fs.temperature(0);
@@ -257,25 +256,16 @@ public:
     { saturation_[phaseIdx] = value; }   
 
     /*!
-     * \brief Set the fugacity coefficient of a component in a phase []
-     *
-     * The phase is always index the same as the component, since we
-     * assume immiscibility.
+     * \brief Set the density of a phase [kg / m^3]
      */
-    void setFugacityCoeff(int compIdx, Scalar value)
-    { fugacityCoeff_[compIdx] = value; }   
+    void setDensity(int phaseIdx, Scalar value)
+    { density_[phaseIdx] = value; }   
 
     /*!
-     * \brief Set the molar volume of a phase [m^3/mol]
+     * \brief Set the specific internal energy of a phase [J/m^3]
      */
-    void setMolarVolume(int phaseIdx, Scalar value)
-    { molarVolume_[phaseIdx] = value; }   
-
-    /*!
-     * \brief Set the specific enthalpy of a phase [J/m^3]
-     */
-    void setEnthalpy(int phaseIdx, Scalar value)
-    { enthalpy_[phaseIdx] = value; }   
+    void setInternalEnergy(int phaseIdx, Scalar value)
+    { internalEnergy_[phaseIdx] = value; }   
 
     /*!
      * \brief Set the dynamic viscosity of a phase [Pa s]
@@ -296,12 +286,12 @@ public:
 #if HAVE_VALGRIND && ! defined NDEBUG
         for (int i = 0; i < numPhases; ++i) {
             //for (int j = 0; j < numComponents; ++j) {
-            //    Valgrind::CheckDefined(fugacityCoeff_[i][j]);
+            //    Valgrind::CheckDefined(fugacityCoefficient_[i][j]);
             //}
             Valgrind::CheckDefined(pressure_[i]);
             Valgrind::CheckDefined(saturation_[i]);
-            Valgrind::CheckDefined(molarVolume_[i]);
-            //Valgrind::CheckDefined(enthalpy_[i]);
+            Valgrind::CheckDefined(density_[i]);
+            //Valgrind::CheckDefined(internalEnergy_[i]);
             Valgrind::CheckDefined(viscosity_[i]);
         }
 
@@ -310,12 +300,10 @@ public:
     }
 
 protected:
-    Scalar fugacityCoeff_[numComponents];
-
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
-    Scalar molarVolume_[numPhases];
-    Scalar enthalpy_[numPhases];
+    Scalar density_[numPhases];
+    Scalar internalEnergy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_;
 };
diff --git a/dumux/material/MpNcfluidstates/trackinggenericfluidstate.hh b/dumux/material/MpNcfluidstates/trackinggenericfluidstate.hh
deleted file mode 100644
index 8960a49dc4386f66d9159310c1101181c5dfac23..0000000000000000000000000000000000000000
--- a/dumux/material/MpNcfluidstates/trackinggenericfluidstate.hh
+++ /dev/null
@@ -1,398 +0,0 @@
-/*****************************************************************************
- *   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 GenericFluidState which keeps track on which variables changed.
- */
-#ifndef DUMUX_TRACKING_GENERIC_FLUID_STATE_HH
-#define DUMUX_TRACKING_GENERIC_FLUID_STATE_HH
-
-#include "genericfluidstate.hh"
-
-#include <dumux/common/valgrind.hh>
-
-#include <cmath>
-#include <algorithm>
-
-namespace Dumux
-{
-/*!
- * \brief A GenericFluidState which keeps track on which variables changed.
- */
-template <class Scalar, class StaticParameters>
-class TrackingGenericFluidState
-{
-    typedef Dumux::GenericFluidState<Scalar, StaticParameters> GenericFluidState;
-public:
-    enum { numComponents = StaticParameters::numComponents };
-    enum { numPhases = StaticParameters::numPhases };
-
-    TrackingGenericFluidState()
-    {
-        for (int i = 0; i < numPhases; ++i)
-            changed_[i].raw = ~((unsigned) 0x00);
-    }
-
-    /*****************************************************
-     * Generic access to fluid properties (No assumptions 
-     * on thermodynamic equilibrium required)
-     *****************************************************/
-    /*!
-     * \brief Returns the saturation of a phase []
-     */
-    Scalar saturation(int phaseIdx) const
-    { return fs_.saturation(phaseIdx); }
-
-    /*!
-     * \brief The mole fraction of a component in a phase []
-     */
-    Scalar moleFrac(int phaseIdx, int compIdx) const
-    { return fs_.moleFrac(phaseIdx, compIdx); }
-
-    /*!
-     * \brief The mass fraction of a component in a phase []
-     */
-    Scalar massFrac(int phaseIdx, int compIdx) const
-    { return fs_.massFrac(phaseIdx, compIdx); }
-
-    /*!
-     * \brief The sum of all component mole fractions in a phase []
-     *
-     * We define this to be the same as the sum of all mass fractions.
-     */
-    Scalar sumMoleFrac(int phaseIdx) const
-    { return fs_.sumMoleFrac(phaseIdx); }
-
-    /*!
-     * \brief The sum of all component mass fractions in a phase []
-     *
-     * We define this to be the same as the sum of all mole fractions.
-     */
-    Scalar sumMassFrac(int phaseIdx) const
-    { return fs_.sumMassFrac(phaseIdx); }
-
-    /*!
-     * \brief The mean molar mass of a fluid phase [kg/mol]
-     */
-    Scalar meanMolarMass(int phaseIdx) const
-    { return fs_.meanMolarMass(phaseIdx); }
-
-    /*!
-     * \brief The molar concentration of a component in a phase [mol/m^3]
-     *
-     * This is often just called "concentration", but there are many
-     * other (though less common) measures for concentration.
-     *
-     * http://en.wikipedia.org/wiki/Concentration
-     */
-    Scalar molarity(int phaseIdx, int compIdx) const
-    { return fs_.molarity(phaseIdx, compIdx); }
-
-    /*!
-     * \brief The fugacity of a component in a phase [Pa]
-     */
-    Scalar fugacity(int phaseIdx, int compIdx) const
-    { return fs_.fugacity(phaseIdx, compIdx); }
-
-    /*!
-     * \brief The fugacity coefficient of a component in a phase []
-     */
-    Scalar fugacityCoeff(int phaseIdx, int compIdx) const
-    { return fs_.fugacityCoeff(phaseIdx, compIdx); }
-
-    /*!
-     * \brief The molar volume of a fluid phase [m^3/mol]
-     */
-    Scalar molarVolume(int phaseIdx) const
-    { return fs_.molarVolume(phaseIdx); }
-
-    /*!
-     * \brief The mass density of a fluid phase [kg/m^3]
-     */
-    Scalar density(int phaseIdx) const
-    { return fs_.density(phaseIdx); }
-
-    /*!
-     * \brief The molar density of a fluid phase [mol/m^3]
-     */
-    Scalar molarDensity(int phaseIdx) const
-    { return fs_.molarDensity(phaseIdx); }
-
-    /*!
-     * \brief The temperature of a fluid phase [K]
-     */
-    Scalar temperature(int phaseIdx) const
-    { return fs_.temperature(phaseIdx); }
-    
-    /*!
-     * \brief The pressure of a fluid phase [Pa]
-     */
-    Scalar pressure(int phaseIdx) const
-    { return fs_.pressure(phaseIdx); }
-
-    /*!
-     * \brief The specific enthalpy of a fluid phase [J/kg]
-     */
-    Scalar enthalpy(int phaseIdx) const
-    { return fs_.enthalpy(phaseIdx); }
-
-    /*!
-     * \brief The specific internal energy of a fluid phase [J/kg]
-     */
-    Scalar internalEnergy(int phaseIdx) const
-    { return fs_.internalEnergy(phaseIdx); }
-
-    /*!
-     * \brief The dynamic viscosity of a fluid phase [Pa s]
-     */
-    Scalar viscosity(int phaseIdx) const
-    { return fs_.viscosity(phaseIdx); };
-
-    /*****************************************************
-     * Setter methods. Note that these are not part of the 
-     * generic FluidState interface but specific for each
-     * implementation...
-     *****************************************************/
-    /*!
-     * \brief Set the temperature [K] of a fluid phase
-     */
-    void setTemperature(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.temperature = 1;
-        fs_.setTemperature(phaseIdx, value);
-    };
-
-    /*!
-     * \brief Set the fluid pressure of a phase [Pa]
-     */
-    void setPressure(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.pressure = 1;
-        fs_.setPressure(phaseIdx, value);
-    };
-
-    /*!
-     * \brief Set the saturation of a phase []
-     */
-    void setSaturation(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.saturation = 1;
-        fs_.setSaturation(phaseIdx, value);
-    };
-
-    /*!
-     * \brief Set the mole fraction of a component in a phase []
-     */
-    void setMoleFrac(int phaseIdx, int compIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.moleFrac |= 1 << compIdx;
-        fs_.setMoleFrac(phaseIdx, compIdx, value);
-    };
-
-
-    /*!
-     * \brief Set the fugacity of a component in a phase []
-     */
-    void setFugacityCoeff(int phaseIdx, int compIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.fugacityCoeff |= 1 << compIdx;
-        fs_.setFugacityCoeff(phaseIdx, compIdx, value);
-    }
-
-    /*!
-     * \brief Set the molar volume of a phase [m^3/mol]
-     */
-    void setMolarVolume(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.molarVolume = 1;
-        fs_.setMolarVolume(phaseIdx, value);
-    }
-
-    /*!
-     * \brief Set the specific enthalpy of a phase [J/m^3]
-     */
-    void setEnthalpy(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.enthalpy = 1;
-        fs_.setEnthalpy(phaseIdx, value);
-    }   
-
-    /*!
-     * \brief Set the dynamic viscosity of a phase [Pa s]
-     */
-    void setViscosity(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.viscosity = 1;
-        fs_.setViscosity(phaseIdx, value);
-    }   
-
-    /*!
-     * \brief Calculatate the mean molar mass of a phase given that
-     *        all mole fractions have been set
-     */
-    void updateMeanMolarMass(int phaseIdx)
-    { fs_.updateMeanMolarMass(phaseIdx); }
-    
-
-    /*!
-     * \brief Reset the dirty tracking for a phase.
-     */
-    void setPhaseClean(int phaseIdx)
-    {
-        // set all tracking bits for the phase to zero
-        changed_[phaseIdx].raw = 0;
-    }
-
-    /*!
-     * \brief Reset the dirty tracking of all phases.
-     */
-    void setClean()
-    {
-        for (int i = 0; i < numPhases; ++i)
-            changed_[i].raw = 0;
-    }
-
-    /*!
-     * \brief Return if all phases are clean.
-     */
-    bool isClean() const
-    { 
-        for (int i = 0; i < numPhases; ++i)
-            if (changed_[i].raw != 0)
-                return false;
-        return true;
-    };
-
-    /*!
-     * \brief Find out whether the values for a phase been since the
-     *        last call to cleanPhase().
-     */
-    bool phaseClean(int phaseIdx) const
-    { return changed_[phaseIdx].raw == 0; }
-
-
-    /*!
-     * \brief Find out whether a specific phase's temperature has been
-     *        changed.
-     */
-    bool temperatureClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.temperature == 0; }
-
-    /*!
-     * \brief Find out whether a specific phase's saturation has been
-     *        changed.
-     */
-    bool saturationClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.saturation == 0; }
-
-    /*!
-     * \brief Find out whether a specific phase's pressure has been
-     *        changed.
-     */
-    bool pressureClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.pressure == 0; }
-
-    /*!
-     * \brief Find out whether a specific phase's enthalpy has been
-     *        changed.
-     */
-    bool enthalpyClean(int phaseIdx) const
-    { return changed_[phaseIdx].enthalpy == 0; }
-
-    /*!
-     * \brief Find out whether a specific phase's temperature has been
-     *        changed.
-     */
-    bool molarVolumeClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.molarVolume == 0; }
-
-    /*!
-     * \brief Find out whether some of the mole fractions of a phase
-     *        have been changed since the last call to
-     *        setCleanPhase().
-     */
-    bool moleFracsClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.moleFrac == 0; }
-
-    /*!
-     * \brief Find out whether a specific mole fraction of a component
-     *        in a phase has been changed since the last call to
-     *        setCleanPhase().
-     */
-    bool moleFracsClean(int phaseIdx, int compIdx) const
-    { return (changed_[phaseIdx].fields.moleFrac  & (1 << compIdx)) == 0; }
-
-    /*!
-     * \brief Find out whether at least one fugacity of a component in
-     *        a phase has been changed since the last call to
-     *        setCleanPhase().
-     */
-    bool fugacityCoeffsClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.fugacitCoeffs == 0; }
-
-    /*!
-     * \brief Find out whether a specific fugacity of a component in a
-     *        phase has been changed since the last call to
-     *        setCleanPhase().
-     */
-    bool fugacityCoeffClean(int phaseIdx, int compIdx) const
-    { return (changed_[phaseIdx].fields.fugacityCoeffs  & (1 << compIdx)) == 0; }
-    
-
-    /*!
-     * \brief Make sure that all attributes are defined.
-     *
-     * This method does not do anything if the program is not run
-     * under valgrind. If it is, then valgrind will print an error
-     * message if some attributes of the object have not been properly
-     * defined.
-     */
-    void checkDefined() const
-    {
-#if HAVE_VALGRIND && ! defined NDEBUG
-        Valgrind::CheckDefined(changed_);
-        fs_.checkDefined();
-#endif // HAVE_VALGRIND
-    }
-
-protected:
-    static_assert(sizeof(unsigned)*8 >= 5 + 2*numComponents,
-                  "Too many components to use the bitfield trick");
-    union {
-        unsigned raw;
-        struct {
-            unsigned char moleFrac : numComponents;
-            unsigned char fugacityCoeff : numComponents;
-
-            unsigned char molarVolume : 1;
-            unsigned char temperature : 1;
-            unsigned char saturation : 1;
-            unsigned char pressure : 1;
-            unsigned char enthalpy : 1;
-            unsigned char viscosity : 1;
-        } __attribute__((__packed__)) fields;
-    } changed_[numPhases];
-    GenericFluidState fs_;
-};
-
-} // end namepace Dumux
-
-#endif
diff --git a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
index c03aca5fe40414bf220cb4e61c935b95cc98e191..283e9385a3557d5193c153bd29aa2acbe9eede98 100644
--- a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
+++ b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
@@ -1,5 +1,5 @@
 /*****************************************************************************
- *   Copyright (C) 2009-2010 by Andreas Lauser                               *
+ *   Copyright (C) 2009-2011 by Andreas Lauser                               *
  *   Institute of Hydraulic Engineering                                      *
  *   University of Stuttgart, Germany                                        *
  *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
@@ -38,14 +38,32 @@
 
 #include <dumux/common/exceptions.hh>
 
+#include "nullparametercache.hh"
+
 namespace Dumux
 {
+/*!
+ * \brief A twophase fluid system with water and nitrogen as components.
+ */
+template <class Scalar>
+class H2ON2FluidSystem
+{
+    // convenience typedefs
+    typedef Dumux::IdealGas<Scalar> IdealGas;
 
+    typedef Dumux::H2O<Scalar> IapwsH2O;
+    typedef Dumux::SimpleH2O<Scalar> SimpleH2O;
+
+    typedef Dumux::TabulatedComponent<Scalar, IapwsH2O > TabulatedH2O;
+
+    typedef Dumux::N2<Scalar> SimpleN2;
+
+public:
+    //! The type of parameter cache objects
+    typedef Dumux::NullParameterCache ParameterCache;
 
-template <class Scalar>
-struct H2ON2StaticParameters {
     /****************************************
-     * Fluid phase parameters
+     * Fluid phase related static parameters
      ****************************************/
 
     //! Number of phases in the fluid system
@@ -55,38 +73,15 @@ struct H2ON2StaticParameters {
     static constexpr int lPhaseIdx = 0;
     //! Index of the gas phase
     static constexpr int gPhaseIdx = 1;
-    //! Index of the solid phase
-    static constexpr int sPhaseIdx = 2;
     
-    //! The component for pure water 
-    typedef Dumux::SimpleH2O<Scalar> SimpleH2O;
-    typedef Dumux::H2O<Scalar> IapwsH2O;
-    typedef Dumux::TabulatedComponent<Scalar, IapwsH2O> TabulatedH2O;
-
-//    typedef IapwsH2O H2O;
+    //! The components for pure water
     typedef TabulatedH2O H2O;
-//    typedef SimpleH2O H2O;
-    
-    //! The component for pure nitrogen
-    typedef Dumux::N2<Scalar> N2;
-
-    /*!
-     * \brief Initialize the static parameters.
-     */
-    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";
-            
-            TabulatedH2O::init(tempMin, tempMax, nTemp,
-                               pressMin, pressMax, nPress);
-        }
+    //typedef SimpleH2O H2O;
+    //typedef IapwsH2O H2O;
 
-    }
-   
+    //! The components for pure nitrogen
+    typedef SimpleN2 N2;
+ 
     /*!
      * \brief Return the human readable name of a fluid phase
      */
@@ -94,11 +89,10 @@ struct H2ON2StaticParameters {
     {
         static const char *name[] = {
             "l",
-            "g",
-            "s"
+            "g"          
         };
 
-        assert(0 <= phaseIdx && phaseIdx < numPhases + 1);
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
         return name[phaseIdx];
     }
 
@@ -133,7 +127,7 @@ struct H2ON2StaticParameters {
     }
 
     /****************************************
-     * Component related parameters
+     * Component related static parameters
      ****************************************/
 
     //! Number of components in the fluid system
@@ -220,81 +214,68 @@ struct H2ON2StaticParameters {
         assert(0 <= compIdx && compIdx < numComponents);
         return accFac[compIdx];
     };
-};
-
-/*!
- * \brief A twophase fluid system with water and nitrogen as components.
- */
-template <class Scalar>
-class H2ON2FluidSystem : public H2ON2StaticParameters<Scalar>
-{
-public:
-    typedef Dumux::H2ON2StaticParameters<Scalar> StaticParameters;
-    typedef Dumux::GenericFluidState<Scalar, StaticParameters> MutableParameters;
-
-private:
-    // convenience typedefs
-    typedef StaticParameters SP;
-    typedef Dumux::IdealGas<Scalar> IdealGas;
 
-public:
+    /*!
+     * \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=*/-10,
+             /*pMax=*/20e6,
+             /*numP=*/200);
+    }
 
     /*!
-     * \brief Initialize the fluid system's static parameters
+     * \brief Initialize the fluid system's static parameters using 
+     *        problem specific temperature and pressure ranges
      */
     static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
-    				 Scalar pressMin, Scalar pressMax, unsigned nPress)
+                     Scalar pressMin, Scalar pressMax, unsigned nPress)
     {
-        SP::init(tempMin, tempMax, nTemp,
-                 pressMin, pressMax, nPress);
+        if (H2O::isTabulated) {
+            std::cout << "Initializing tables for the H2O fluid properties ("
+                      << nTemp*nPress
+                      << " entries).\n";
+            
+            TabulatedH2O::init(tempMin, tempMax, nTemp,
+                               pressMin, pressMax, nPress);
+        }
     }
 
     /*!
      * \brief Calculate the molar volume [m^3/mol] of a fluid phase
      */
-    static Scalar computeMolarVolume(MutableParameters &params,
-                                     int phaseIdx)
+    template <class FluidState>
+    static Scalar density(const FluidState &fluidState,
+                          const ParameterCache &paramCache,
+                          int phaseIdx)
     {
-        assert(0 <= phaseIdx  && phaseIdx <= SP::numPhases);
+        assert(0 <= phaseIdx  && phaseIdx <= numPhases);
 
-        params.updateMeanMolarMass(phaseIdx);
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);
+        Scalar T = fluidState.temperature(phaseIdx);
+        Scalar p = fluidState.pressure(phaseIdx);
 
         switch (phaseIdx) {
-        case SP::lPhaseIdx:
+        case lPhaseIdx:
             // assume pure water where one water molecule gets
             // replaced by one nitrogen molecule
-            return SP::H2O::molarMass()/SP::H2O::liquidDensity(T, p);
-        case SP::gPhaseIdx:
+            return H2O::liquidDensity(T, p);
+        case gPhaseIdx:
             // assume ideal gas
-            return 1.0 / IdealGas::concentration(T, p);
+            return
+                IdealGas::molarDensity(T, p) 
+                * fluidState.averageMolarMass(gPhaseIdx);
         }
         
         DUNE_THROW(Dune::InvalidStateException, "Unhandled phase index " << phaseIdx);
     };
 
-    /*!
-     * \brief Calculate the fugacity of a component in a fluid phase
-     *        [Pa]
-     *
-     * The components chemical \f$mu_\kappa\f$ potential 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.
-     */
-    static Scalar computeFugacity(const MutableParameters &params,
-                                  int phaseIdx,
-                                  int compIdx)
-    {
-        Scalar x = params.moleFrac(phaseIdx,compIdx);
-        Scalar p = params.pressure(phaseIdx);
-        return x*p*computeFugacityCoeff(params, phaseIdx, compIdx);
-    };
-
     /*!
      * \brief Calculate the fugacity coefficient [Pa] of an individual
      *        component in a fluid phase
@@ -305,23 +286,24 @@ public:
      *
      * \f[ f_\kappa = \phi_\kappa * x_{\kappa} \f]
      */
-    static Scalar computeFugacityCoeff(MutableParameters &params,
-                                       int phaseIdx,
-                                       int compIdx)
+    template <class FluidState>
+    static Scalar fugacityCoefficient(const FluidState &fluidState,
+                                      const ParameterCache &paramCache,
+                                      int phaseIdx,
+                                      int compIdx)
     {
-        assert(0 <= phaseIdx  && phaseIdx <= SP::numPhases);
-        assert(0 <= compIdx  && compIdx <= SP::numComponents);
+        assert(0 <= phaseIdx  && phaseIdx <= numPhases);
+        assert(0 <= compIdx  && compIdx <= numComponents);
 
-        params.updateMeanMolarMass(phaseIdx);
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);
+        Scalar T = fluidState.temperature(phaseIdx);
+        Scalar p = fluidState.pressure(phaseIdx);
         switch (phaseIdx) {
-        case SP::lPhaseIdx: 
+        case lPhaseIdx: 
             switch (compIdx) {
-            case SP::H2OIdx: return SP::H2O::vaporPressure(T)/p;
-            case SP::N2Idx: return BinaryCoeff::H2O_N2::henry(T)/p;
+            case H2OIdx: return H2O::vaporPressure(T)/p;
+            case N2Idx: return BinaryCoeff::H2O_N2::henry(T)/p;
             };
-        case SP::gPhaseIdx:
+        case gPhaseIdx:
             return 1.0; // ideal gas
         };
 
@@ -331,21 +313,22 @@ public:
     /*!
      * \brief Calculate the dynamic viscosity of a fluid phase [Pa*s]
      */
-    static Scalar computeViscosity(MutableParameters &params,
-                                   int phaseIdx)
+    template <class FluidState>
+    static Scalar viscosity(const FluidState &fluidState,
+                            const ParameterCache &paramCache,
+                            int phaseIdx)
     {
-        assert(0 <= phaseIdx  && phaseIdx <= SP::numPhases);
+        assert(0 <= phaseIdx  && phaseIdx <= numPhases);
 
-        params.updateMeanMolarMass(phaseIdx);
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);
+        Scalar T = fluidState.temperature(phaseIdx);
+        Scalar p = fluidState.pressure(phaseIdx);
         switch (phaseIdx) {
-        case SP::lPhaseIdx:
+        case lPhaseIdx:
             // assume pure water for the liquid phase
-            return SP::H2O::liquidViscosity(T, p);
-        case SP::gPhaseIdx:
+            return H2O::liquidViscosity(T, p);
+        case gPhaseIdx:
             // assume pure water for the gas phase
-            return SP::N2::gasViscosity(T, p);
+            return N2::gasViscosity(T, p);
         }
         
         DUNE_THROW(Dune::InvalidStateException, "Unhandled phase index " << phaseIdx);
@@ -370,9 +353,11 @@ public:
      * where \f$p_\alpha\f$ and \f$T_\alpha\f$ are the fluid phase'
      * pressure and temperature.
      */
-    static Scalar computeDiffusionCoeff(const MutableParameters &params,
-                                        int phaseIdx,
-                                        int compIdx)
+    template <class FluidState>
+    static Scalar diffusionCoefficient(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx,
+                                       int compIdx)
     {
         // TODO!
         DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients");
@@ -383,20 +368,21 @@ public:
      *        return the binary diffusion coefficient for components
      *        \f$i\f$ and \f$j\f$ in this phase.
      */
-    static Scalar computeBinaryDiffCoeff(MutableParameters &params, 
-                                         int phaseIdx,
-                                         int compIIdx,
-                                         int compJIdx)
+    template <class FluidState>
+    static Scalar binaryDiffusionCoefficient(const FluidState &fluidState, 
+                                             const ParameterCache &paramCache,
+                                             int phaseIdx,
+                                             int compIIdx,
+                                             int compJIdx)
                                   
     {
-        params.updateMeanMolarMass(phaseIdx);
         if (compIIdx > compJIdx)
             std::swap(compIIdx, compJIdx);
 
 #ifndef NDEBUG
         if (compIIdx == compJIdx ||
-            phaseIdx > SP::numPhases - 1 ||
-            compJIdx > SP::numComponents - 1)
+            phaseIdx > numPhases - 1 ||
+            compJIdx > numComponents - 1)
         {
             DUNE_THROW(Dune::InvalidStateException,
                        "Binary diffusion coefficient of components "
@@ -405,26 +391,26 @@ public:
         }
 #endif
 
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);       
+        Scalar T = fluidState.temperature(phaseIdx);
+        Scalar p = fluidState.pressure(phaseIdx);       
 
         switch (phaseIdx) {
-        case SP::lPhaseIdx:
+        case lPhaseIdx:
             switch (compIIdx) {
-            case SP::H2OIdx:
+            case H2OIdx:
                 switch (compJIdx) {
-                case SP::N2Idx: return BinaryCoeff::H2O_N2::liquidDiffCoeff(T, p);
+                case N2Idx: return BinaryCoeff::H2O_N2::liquidDiffCoeff(T, p);
                 }
             default:
                 DUNE_THROW(Dune::InvalidStateException,
                            "Binary diffusion coefficients of trace "
                            "substances in liquid phase is undefined!\n");
             }
-        case SP::gPhaseIdx:
+        case gPhaseIdx:
             switch (compIIdx) {
-            case SP::H2OIdx:
+            case H2OIdx:
                 switch (compJIdx) {
-                case SP::N2Idx: return BinaryCoeff::H2O_N2::gasDiffCoeff(T, p);
+                case N2Idx: return BinaryCoeff::H2O_N2::gasDiffCoeff(T, p);
                 }
             }
         }
@@ -438,98 +424,87 @@ public:
     /*!
      * \brief Given a phase's composition, temperature, pressure and
      *        density, calculate its specific enthalpy [J/kg].
+     *
+     *  \todo This fluid 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. ...
      */
-
-    /*!
-     *  \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. ...
-     */
-    static Scalar computeEnthalpy(MutableParameters &params, 
-                                  int phaseIdx)
+    template <class FluidState>
+    static Scalar internalEnergy(const FluidState &fluidState, 
+                                 const ParameterCache &paramCache,
+                                 int phaseIdx)
     {
-        params.updateMeanMolarMass(phaseIdx);
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);
+        Scalar T = fluidState.temperature(phaseIdx);
+        Scalar p = fluidState.pressure(phaseIdx);
         Valgrind::CheckDefined(T);
         Valgrind::CheckDefined(p);
-        if (phaseIdx == SP::lPhaseIdx) {
+        if (phaseIdx == lPhaseIdx) {
             // TODO: correct way to deal with the solutes???
             return 
-                SP::H2O::liquidEnthalpy(T, p) ;
+                H2O::liquidInternalEnergy(T, p) ;
         }
         else {
             // assume ideal gas
-            Scalar XH2O = params.massFrac(SP::gPhaseIdx, SP::H2OIdx);
-            Scalar XN2 = params.massFrac(SP::gPhaseIdx, SP::N2Idx);          
+            Scalar XH2O = fluidState.massFrac(gPhaseIdx, H2OIdx);
+            Scalar XN2 = fluidState.massFrac(gPhaseIdx, N2Idx);          
             Scalar result = 0;
-            result += XH2O*SP::H2O::gasEnthalpy(T, p);
-            result += XN2*SP::N2::gasEnthalpy(T, p);
+            result += XH2O*H2O::gasInternalEnergy(T, p);
+            result += XN2*N2::gasInternalEnergy(T, p);
             return result;
         }
     }
 
     /*!
-     * \brief Given a phase's composition, temperature, pressure and
-     *        density, calculate its specific internal energy [J/kg].
-     */
-    static Scalar computeInternalEnergy(MutableParameters &params, 
-                                        int phaseIdx)
-    {
-        params.updateMeanMolarMass(phaseIdx);
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);
-        Scalar rho = params.density(phaseIdx);
-        return
-            computeEnthalpy(params, phaseIdx) -
-            p/rho;
-    }
-
-    /*!
-     * \brief Thermal conductivity of phases.
+     * \brief Thermal conductivity of a fluid phase.
      *
      * Use the conductivity of air and water as a first approximation.
      * Source:
      * http://en.wikipedia.org/wiki/List_of_thermal_conductivities
      */
-    static Scalar computeThermalConductivity(const MutableParameters &params,
-                                             int phaseIdx)
+    template <class FluidState>
+    static Scalar thermalConductivity(const FluidState &fluidState,
+                                      const ParameterCache &paramCache,
+                                      int phaseIdx)
     {
 //    	TODO thermal conductivity is a function of:
-//        Scalar p = params.pressure(phaseIdx);
-//        Scalar T = params.temperature(phaseIdx);
-//        Scalar x = params.moleFrac(phaseIdx,compIdx);
+//        Scalar p = fluidState.pressure(phaseIdx);
+//        Scalar T = fluidState.temperature(phaseIdx);
+//        Scalar x = fluidState.moleFrac(phaseIdx,compIdx);
 #warning: so far rough estimates from wikipedia
         switch (phaseIdx) {
-        case SP::lPhaseIdx: // use conductivity of pure water
+        case lPhaseIdx: // use conductivity of pure water
             return  0.6;   // conductivity of water[W / (m K ) ]
-        case SP::gPhaseIdx:// use conductivity of pure air
+        case gPhaseIdx:// use conductivity of pure air
             return 0.025; // conductivity of air [W / (m K ) ]
         }
         DUNE_THROW(Dune::InvalidStateException, "Unhandled phase index " << phaseIdx);
     }
     
     /*!
-     * \brief Specific isobaric heat capacity of liquid water / air
+     * \brief Specific isobaric heat capacity of a fluid phase.
      *        \f$\mathrm{[J/kg]}\f$.
      *
      * \param params    mutable parameters
      * \param phaseIdx  for which phase to give back the heat capacity
      */
-    static Scalar computeHeatCapacity(const MutableParameters &params,
-                                      int phaseIdx)
+    template <class FluidState>
+    static Scalar heatCapacity(const FluidState &fluidState,
+                               const ParameterCache &paramCache,
+                               int phaseIdx)
     {
 //        http://en.wikipedia.org/wiki/Heat_capacity
 #warning: so far rough estimates from wikipedia
 //      TODO heatCapacity is a function of composition.
-//        Scalar p = params.pressure(phaseIdx);
-//        Scalar T = params.temperature(phaseIdx);
-//        Scalar x = params.moleFrac(phaseIdx,compIdx);
+//        Scalar p = fluidState.pressure(phaseIdx);
+//        Scalar T = fluidState.temperature(phaseIdx);
+//        Scalar x = fluidState.moleFrac(phaseIdx,compIdx);
         switch (phaseIdx) {
-        case SP::lPhaseIdx: // use heat capacity of pure liquid water
+        case lPhaseIdx: // use heat capacity of pure liquid water
             return  4181.3;  // @(25°C) !!!
             /* [J/(kg K)]*/ /* not working because ddgamma_ddtau is not defined*/ /* Dumux::H2O<Scalar>::liquidHeatCap_p(T,
                                              p); */
-        case SP::gPhaseIdx:
+        case gPhaseIdx:
             return  1003.5 ; // @ (0°C) !!!
             /* [J/(kg K)]*/ /* not working because ddgamma_ddtau is not defined*/ /*Dumux::H2O<Scalar>::gasHeatCap_p(T,
                                           p) ;*/
diff --git a/dumux/material/MpNcfluidsystems/nullparametercache.hh b/dumux/material/MpNcfluidsystems/nullparametercache.hh
new file mode 100644
index 0000000000000000000000000000000000000000..25395197802d42dbd9226e6b5e4052d3720e53e3
--- /dev/null
+++ b/dumux/material/MpNcfluidsystems/nullparametercache.hh
@@ -0,0 +1,56 @@
+/*****************************************************************************
+ *   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 The a parameter cache which does nothing
+ */
+#ifndef DUMUX_NULL_PARAMETER_CACHE_HH
+#define DUMUX_NULL_PARAMETER_CACHE_HH
+
+#include "parametercachebase.hh"
+
+namespace Dumux
+{
+/*!
+ * \brief The a parameter cache which does nothing
+ */
+class NullParameterCache : public ParameterCacheBase<NullParameterCache>
+{
+public:
+    NullParameterCache() 
+    {};
+    
+    template <class FluidState>
+    void updateAll(const FluidState &fs)
+    {
+    };
+    
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     */
+    template <class FluidState>
+    void updatePhase(const FluidState &fs, int phaseIdx)
+    {};
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/material/MpNcfluidsystems/parametercachebase.hh b/dumux/material/MpNcfluidsystems/parametercachebase.hh
new file mode 100644
index 0000000000000000000000000000000000000000..71b2984f6042287b3c53e00ce41b049feaba07d5
--- /dev/null
+++ b/dumux/material/MpNcfluidsystems/parametercachebase.hh
@@ -0,0 +1,120 @@
+/*****************************************************************************
+ *   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 The base class of the parameter cache classes for fluid systems
+ */
+#ifndef DUMUX_PARAMETER_CACHE_BASE_HH
+#define DUMUX_PARAMETER_CACHE_BASE_HH
+
+namespace Dumux
+{
+/*!
+ * \brief The base class of the parameter cache classes for fluid systems
+ */
+template <class Implementation>
+class ParameterCacheBase
+{
+public:
+    ParameterCacheBase() 
+    {};
+    
+    template <class FluidState>
+    void updateAll(const FluidState &fs)
+    {
+        for (int phaseIdx = 0; phaseIdx < FluidState::numPhases; ++phaseIdx)
+            updatePhase(fs, phaseIdx);
+    };
+    
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     */
+    template <class FluidState>
+    void updatePhase(const FluidState &fs, int phaseIdx)
+    {};
+
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     *        which depend on temperature
+     *
+     * *Only* use this method if only the temperature of a phase
+     * changed between two update*() calls. If more changed, call
+     * updatePhase()!
+     */
+    template <class FluidState> void
+    updatePhaseTemperature(const FluidState &fs, int phaseIdx)
+    {
+        asImp_().updatePhase(fs, phaseIdx);
+    };
+
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     *        which depend on pressure
+     *
+     * *Only* use this method if only the pressure of a phase changed
+     * between two update*() calls. If more changed, call
+     * updatePhase()!
+     */
+    template <class FluidState>
+    void updatePhasePressure(const FluidState &fs, int phaseIdx)
+    {
+        asImp_().updatePhase(fs, phaseIdx);
+    };
+
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     *        which depend on composition
+     *
+     * *Only* use this method if neither the pressure nor the
+     * temperature of the phase changed between two update*()
+     * calls. If more changed, call updatePhase()!
+     */
+    template <class FluidState>
+    void updatePhaseComposition(const FluidState &fs, int phaseIdx)
+    {
+        asImp_().updatePhase(fs, phaseIdx);
+    };
+
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     *        which depend on the mole fraction of a single component
+     *
+     * *Only* use this method if just a single component's
+     * concentration changed between two update*() calls. If more than
+     * one concentration changed, call updatePhaseComposition() of
+     * updatePhase()!
+     */
+    template <class FluidState>
+    void updatePhaseSingleMoleFraction(const FluidState &fs, 
+                                       int phaseIdx,
+                                       int compIdx)
+    {
+        asImp_().updatePhaseComposition(fs, phaseIdx);
+    };
+
+private:
+    Implementation &asImp_()
+    { return *static_cast<Implementation*>(this); } 
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/material/idealgas.hh b/dumux/material/idealgas.hh
index 8aaba763a05feb660d558b4697379da328172918..66a60313c50031436d7ddb9afda335ea64bf37ee 100644
--- a/dumux/material/idealgas.hh
+++ b/dumux/material/idealgas.hh
@@ -51,19 +51,24 @@ public:
 
     /*!
      * \brief The pressure of the gas in \f$\mathrm{[N/m^2]}\f$, depending on
-     *        concentration and temperature.
+     *        the molar density and temperature.
      */
     static Scalar pressure(Scalar temperature,
-                           Scalar concentration)
-    { return R*temperature*concentration; }
+                           Scalar rhoMolar)
+    { return R*temperature*rhoMolar; }
 
     /*!
-     * \brief The molar concentration of the gas in \f$\mathrm{[mol/m^3]}\f$, depending on
-     *        pressure and temperature.
+     * \brief The molar density of the gas \f$\mathrm{[mol/m^3]}\f$,
+     *        depending on pressure and temperature.
      */
-    static Scalar concentration(Scalar temperature,
+    static Scalar molarDensity(Scalar temperature,
                                 Scalar pressure)
     { return pressure/(R*temperature); }
+
+    DUMUX_DEPRECATED_MSG("use molarDensity()")
+    static Scalar concentration(Scalar temperature,
+                                Scalar pressure)
+    { return molarDensity(temperature, pressure); }
 };
 
 } // end namepace