diff --git a/dumux/porousmediumflow/mpnc/implicit/diffusion/fluxvariables.hh b/dumux/porousmediumflow/mpnc/implicit/diffusion/fluxvariables.hh
index 233f6f3e3bbb0774b9d3868f961d2be9e8fdc626..ce4b0ea710f2718e48db68fb0aacddab786bb07c 100644
--- a/dumux/porousmediumflow/mpnc/implicit/diffusion/fluxvariables.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/diffusion/fluxvariables.hh
@@ -105,12 +105,12 @@ public:
                 for (int compIdx = 0; compIdx < numComponents; ++compIdx){
 
                     // calculate mole fractions at the integration points of the face
-                    moleFraction_[phaseIdx][compIdx] += elemVolVars[volVarsIdx].fluidState().moleFraction(phaseIdx, compIdx)*
+                    moleFraction_[phaseIdx][compIdx] += elemVolVars[volVarsIdx].moleFraction(phaseIdx, compIdx)*
                     face.shapeValue[idx];
 
                     // calculate mole fraction gradients
                     tmp =  feGrad;
-                    tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(phaseIdx, compIdx);
+                    tmp *= elemVolVars[volVarsIdx].moleFraction(phaseIdx, compIdx);
                     moleFractionGrad_[phaseIdx][compIdx] += tmp;
                 }
             }
@@ -129,8 +129,8 @@ public:
         {
             // make sure to only calculate diffusion coefficents
             // for phases which exist in both finite volumes
-            if (elemVolVars[i].fluidState().saturation(phaseIdx) <= 1e-4 ||
-                elemVolVars[j].fluidState().saturation(phaseIdx) <= 1e-4)
+            if (elemVolVars[i].saturation(phaseIdx) <= 1e-4 ||
+                elemVolVars[j].saturation(phaseIdx) <= 1e-4)
             {
                 continue;
             }
@@ -141,11 +141,11 @@ public:
             // and j
             //
             Scalar red_i =
-                elemVolVars[i].fluidState().saturation(phaseIdx)/elemVolVars[i].porosity() *
-                pow(elemVolVars[i].porosity() * elemVolVars[i].fluidState().saturation(phaseIdx), 7.0/3);
+                elemVolVars[i].saturation(phaseIdx)/elemVolVars[i].porosity() *
+                pow(elemVolVars[i].porosity() * elemVolVars[i].saturation(phaseIdx), 7.0/3);
             Scalar red_j =
-                elemVolVars[j].fluidState().saturation(phaseIdx)/elemVolVars[j].porosity() *
-                pow(elemVolVars[j].porosity() * elemVolVars[j].fluidState().saturation(phaseIdx), 7.0/3);
+                elemVolVars[j].saturation(phaseIdx)/elemVolVars[j].porosity() *
+                pow(elemVolVars[j].porosity() * elemVolVars[j].saturation(phaseIdx), 7.0/3);
 
             if (phaseIdx == wPhaseIdx) {
                 // Liquid phase diffusion coefficients in the porous medium
diff --git a/dumux/porousmediumflow/mpnc/implicit/energy/fluxvariables.hh b/dumux/porousmediumflow/mpnc/implicit/energy/fluxvariables.hh
index dca4d60f006688a36de142b4e716ad4f023cd3f4..265bef607b35081d7e9590dfdc723afcde862665 100644
--- a/dumux/porousmediumflow/mpnc/implicit/energy/fluxvariables.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/energy/fluxvariables.hh
@@ -141,7 +141,7 @@ public:
             // index for the element volume variables
             int volVarsIdx = face.fapIndices[idx];
 
-            tmp *= elemVolVars[volVarsIdx].fluidState().temperature(/*phaseIdx=*/0);
+            tmp *= elemVolVars[volVarsIdx].temperature(/*phaseIdx=*/0);
             temperatureGradient += tmp;
         }
 
@@ -192,13 +192,13 @@ protected:
         if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
         {
             lambdaI =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].fluidState().saturation(wPhaseIdx),
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().solidThermalConductivity(element, fvGeometry, i),
                                                                    problem.spatialParams().porosity(element, fvGeometry, i));
             lambdaJ =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].fluidState().saturation(wPhaseIdx),
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().solidThermalConductivity(element, fvGeometry, j),
@@ -211,7 +211,7 @@ protected:
             fvGeometryI.subContVol[0].global = elementI.geometry().center();
 
             lambdaI =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].fluidState().saturation(wPhaseIdx),
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().solidThermalConductivity(elementI, fvGeometryI, 0),
@@ -222,7 +222,7 @@ protected:
             fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
 
             lambdaJ =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].fluidState().saturation(wPhaseIdx),
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().solidThermalConductivity(elementJ, fvGeometryJ, 0),
diff --git a/dumux/porousmediumflow/mpnc/implicit/energy/fluxvariableskinetic.hh b/dumux/porousmediumflow/mpnc/implicit/energy/fluxvariableskinetic.hh
index f1595f2afe2de229f77fc12969b325ab61bde1cb..4d45188a5e06c1ff738f4641f8896f5b3af4e44c 100644
--- a/dumux/porousmediumflow/mpnc/implicit/energy/fluxvariableskinetic.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/energy/fluxvariableskinetic.hh
@@ -246,13 +246,13 @@ protected:
         if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
         {
             lambdaI =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].fluidState().saturation(wPhaseIdx),
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().solidThermalConductivity(element, fvGeometry, i),
                                                                    problem.spatialParams().porosity(element, fvGeometry, i));
             lambdaJ =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].fluidState().saturation(wPhaseIdx),
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().solidThermalConductivity(element, fvGeometry, j),
@@ -265,7 +265,7 @@ protected:
             fvGeometryI.subContVol[0].global = elementI.geometry().center();
 
             lambdaI =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].fluidState().saturation(wPhaseIdx),
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[i].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().solidThermalConductivity(elementI, fvGeometryI, 0),
@@ -276,7 +276,7 @@ protected:
             fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
 
             lambdaJ =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].fluidState().saturation(wPhaseIdx),
+                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(wPhaseIdx),
                                                                    elemVolVars[j].thermalConductivity(nPhaseIdx),
                                                                    problem.spatialParams().solidThermalConductivity(elementJ, fvGeometryJ, 0),
diff --git a/dumux/porousmediumflow/mpnc/implicit/energy/localresidual.hh b/dumux/porousmediumflow/mpnc/implicit/energy/localresidual.hh
index 3d09f690b48edbda1df0b715b61ae599d7421962..dcc3cbee521b304c845de287b01ec9aa6fe3c6c7 100644
--- a/dumux/porousmediumflow/mpnc/implicit/energy/localresidual.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/energy/localresidual.hh
@@ -161,7 +161,7 @@ public:
 
         // heat stored in the rock matrix
         storage[energyEqIdx] +=
-            volVars.fluidState().temperature(/*phaseIdx=*/0)
+            volVars.temperature(/*phaseIdx=*/0)
             * volVars.solidDensity()
             * (1.0 - volVars.porosity())
             * volVars.solidHeatCapacity();
@@ -248,7 +248,7 @@ if (!std::isfinite(storage[energyEqIdx]))
         // use the phase enthalpy of the upstream vertex to calculate
         // the enthalpy transport
         const VolumeVariables &up = elemVolVars[upIdx];
-        flux[energyEqIdx] += up.fluidState().enthalpy(phaseIdx) * massFlux;
+        flux[energyEqIdx] += up.enthalpy(phaseIdx) * massFlux;
 #ifndef NDEBUG
 if (!std::isfinite(flux[energyEqIdx]) )
     DUNE_THROW(NumericalProblem, "Calculated non-finite energy flux");
diff --git a/dumux/porousmediumflow/mpnc/implicit/energy/localresidualkinetic.hh b/dumux/porousmediumflow/mpnc/implicit/energy/localresidualkinetic.hh
index 0f31f943244b9b40b33398889ae11ce20d4fe56d..2081a0b428b6ab49bef738662b67dbfb7f89d793 100644
--- a/dumux/porousmediumflow/mpnc/implicit/energy/localresidualkinetic.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/energy/localresidualkinetic.hh
@@ -209,8 +209,8 @@ public:
          * To be more precise this should be the components enthalpy.
          * In the same vein: Counter current diffusion is not accounted for here.
          */
-        const Scalar transportedThingUp =  up.fluidState().enthalpy(phaseIdx) ;
-        const Scalar transportedThingDn =  dn.fluidState().enthalpy(phaseIdx) ;
+        const Scalar transportedThingUp =  up.enthalpy(phaseIdx) ;
+        const Scalar transportedThingDn =  dn.enthalpy(phaseIdx) ;
 
 
         const Scalar massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
@@ -574,8 +574,8 @@ public:
         /* todo
          * CAUTION: this is not exactly correct: does diffusion carry the upstream phase enthalpy? To be more precise this should be the components enthalpy. In the same vein: Counter current diffusion is not accounted for here.
          */
-        const Scalar transportedThingUp =  up.fluidState().enthalpy(phaseIdx) ;
-        const Scalar transportedThingDn =  dn.fluidState().enthalpy(phaseIdx) ;
+        const Scalar transportedThingUp =  up.enthalpy(phaseIdx) ;
+        const Scalar transportedThingDn =  dn.enthalpy(phaseIdx) ;
 
         const Scalar massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
         flux[energyEq0Idx] += massFlux *
diff --git a/dumux/porousmediumflow/mpnc/implicit/energy/vtkwriter.hh b/dumux/porousmediumflow/mpnc/implicit/energy/vtkwriter.hh
index abab6bbae5816477b24b97ffc797a03fa63ba02b..51cabb317e1ff8b8eb14e964243fc13b5e364762 100644
--- a/dumux/porousmediumflow/mpnc/implicit/energy/vtkwriter.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/energy/vtkwriter.hh
@@ -100,7 +100,7 @@ public:
             const VolumeVariables &volVars = elemVolVars[scvIdx];
 
             if (temperatureOutput_)
-                temperature_[dofIdxGlobal] = volVars.fluidState().temperature(/*phaseIdx=*/0);
+                temperature_[dofIdxGlobal] = volVars.temperature(/*phaseIdx=*/0);
         }
     }
 
@@ -190,12 +190,12 @@ public:
             int gobalIdx = this->problem_.model().dofMapper().subIndex(element, scvIdx, dofCodim);
             const VolumeVariables &volVars = elemVolVars[scvIdx];
 
-            if (temperatureOutput_) temperature_[gobalIdx] = volVars.fluidState().temperature(/*phaseIdx=*/0);
+            if (temperatureOutput_) temperature_[gobalIdx] = volVars.temperature(/*phaseIdx=*/0);
             for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
                 if (enthalpyOutput_)
-                    enthalpy_[phaseIdx][gobalIdx] = volVars.fluidState().enthalpy(phaseIdx);
+                    enthalpy_[phaseIdx][gobalIdx] = volVars.enthalpy(phaseIdx);
                 if (internalEnergyOutput_)
-                    internalEnergy_[phaseIdx][gobalIdx] = volVars.fluidState().internalEnergy(phaseIdx);
+                    internalEnergy_[phaseIdx][gobalIdx] = volVars.internalEnergy(phaseIdx);
             }
         }
     }
diff --git a/dumux/porousmediumflow/mpnc/implicit/energy/vtkwriterkinetic.hh b/dumux/porousmediumflow/mpnc/implicit/energy/vtkwriterkinetic.hh
index a90707c50e7788a99678a9075a21b7ae2ee40ae1..b811f641b12cff750ffcc648df7a5bf762d407b7 100644
--- a/dumux/porousmediumflow/mpnc/implicit/energy/vtkwriterkinetic.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/energy/vtkwriterkinetic.hh
@@ -136,8 +136,8 @@ public:
             const VolumeVariables &volVars = elemVolVars[localVertexIdx];
 
             for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                enthalpy_[phaseIdx][vIdxGlobal]          = volVars.fluidState().enthalpy(phaseIdx);
-                internalEnergy_[phaseIdx][vIdxGlobal]    = volVars.fluidState().internalEnergy(phaseIdx);
+                enthalpy_[phaseIdx][vIdxGlobal]          = volVars.enthalpy(phaseIdx);
+                internalEnergy_[phaseIdx][vIdxGlobal]    = volVars.internalEnergy(phaseIdx);
                 reynoldsNumber_[phaseIdx][vIdxGlobal]    = volVars.reynoldsNumber(phaseIdx);
                 prandtlNumber_[phaseIdx][vIdxGlobal]     = volVars.prandtlNumber(phaseIdx);
                 nusseltNumber_[phaseIdx][vIdxGlobal]             = volVars.nusseltNumber(phaseIdx);
@@ -377,12 +377,12 @@ public:
             const unsigned int vIdxGlobal = this->problem_.vertexMapper().subIndex(element, localVertexIdx, dim);
             const VolumeVariables &volVars = elemVolVars[localVertexIdx];
 
-            qBoil_[vIdxGlobal] = LocalResidual::QBoilFunc(volVars, volVars.fluidState().saturation(wPhaseIdx));
+            qBoil_[vIdxGlobal] = LocalResidual::QBoilFunc(volVars, volVars.saturation(wPhaseIdx));
             qsf_[vIdxGlobal] = LocalResidual::qsf(volVars);
 
             for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                enthalpy_[phaseIdx][vIdxGlobal]          = volVars.fluidState().enthalpy(phaseIdx);
-                internalEnergy_[phaseIdx][vIdxGlobal]    = volVars.fluidState().internalEnergy(phaseIdx);
+                enthalpy_[phaseIdx][vIdxGlobal]          = volVars.enthalpy(phaseIdx);
+                internalEnergy_[phaseIdx][vIdxGlobal]    = volVars.internalEnergy(phaseIdx);
                 reynoldsNumber_[phaseIdx][vIdxGlobal]    = volVars.reynoldsNumber(phaseIdx);
                 prandtlNumber_[phaseIdx][vIdxGlobal]     = volVars.prandtlNumber(phaseIdx);
                 nusseltNumber_[phaseIdx][vIdxGlobal]             = volVars.nusseltNumber(phaseIdx);
diff --git a/dumux/porousmediumflow/mpnc/implicit/mass/localresidual.hh b/dumux/porousmediumflow/mpnc/implicit/mass/localresidual.hh
index cf2030600aa6989d2010b80433b6acc6e99cc2a5..7ad8af1e1b6c8473887a886a6fd90f96e4f320b1 100644
--- a/dumux/porousmediumflow/mpnc/implicit/mass/localresidual.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/mass/localresidual.hh
@@ -76,8 +76,8 @@ public:
         storage = 0;
         for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
             storage[compIdx] +=
-                volVars.fluidState().saturation(phaseIdx)*
-                volVars.fluidState().molarity(phaseIdx, compIdx);
+                volVars.saturation(phaseIdx)*
+                volVars.molarity(phaseIdx, compIdx);
 #ifndef NDEBUG
 if (!std::isfinite(storage[compIdx]))
     DUNE_THROW(NumericalProblem, "Calculated non-finite storage");
@@ -131,10 +131,10 @@ if (!std::isfinite(volumeFlux))
             if (enableSmoothUpwinding_) {
                 const Scalar kGradPNormal   = fluxVars.kGradPNormal(phaseIdx);
                 const Scalar mobUp          = up.mobility(phaseIdx);
-                const Scalar conUp          = up.fluidState().molarity(phaseIdx, compIdx);
+                const Scalar conUp          = up.molarity(phaseIdx, compIdx);
 
                 const Scalar mobDn  = dn.mobility(phaseIdx);
-                const Scalar conDn  = dn.fluidState().molarity(phaseIdx, compIdx);
+                const Scalar conDn  = dn.molarity(phaseIdx, compIdx);
 
                 const Scalar mobConUp   = mobUp*conUp;
                 const Scalar mobConDn   = mobDn*conDn;
@@ -144,8 +144,8 @@ if (!std::isfinite(volumeFlux))
                 const Scalar sign   = (kGradPNormal > 0)?-1:1;
 
                 // approximate the mean viscosity at the face
-                const Scalar meanVisc = (up.fluidState().viscosity(phaseIdx) +
-                                   dn.fluidState().viscosity(phaseIdx))/2;
+                const Scalar meanVisc = (up.viscosity(phaseIdx) +
+                                   dn.viscosity(phaseIdx))/2;
 
                 // put the mean viscosity and permeanbility in
                 // relation to the viscosity of water at
@@ -180,11 +180,11 @@ if (!std::isfinite(volumeFlux))
             {// not use smooth upwinding
                 flux[compIdx] =
                         volumeFlux *
-                        ((     massUpwindWeight)*up.fluidState().molarity(phaseIdx, compIdx)
+                        ((     massUpwindWeight)*up.molarity(phaseIdx, compIdx)
                                 +
-                        (  1. - massUpwindWeight)*dn.fluidState().molarity(phaseIdx, compIdx) );
+                        (  1. - massUpwindWeight)*dn.molarity(phaseIdx, compIdx) );
                         if (!std::isfinite(flux[compIdx]))
-                            DUNE_THROW(NumericalProblem, "Calculated non-finite normal flux in phase " <<  phaseIdx << " comp " << compIdx << "T: "<<  up.fluidState().temperature(phaseIdx) << "S "<<up.fluidState().saturation(phaseIdx)  ) ;
+                            DUNE_THROW(NumericalProblem, "Calculated non-finite normal flux in phase " <<  phaseIdx << " comp " << compIdx << "T: "<<  up.temperature(phaseIdx) << "S "<<up.saturation(phaseIdx)  ) ;
             }
         }
     }
@@ -211,8 +211,8 @@ if (!std::isfinite(volumeFlux))
 
         const VolumeVariables &volVarsI = fluxVars.volVars(fluxVars.face().i);
         const VolumeVariables &volVarsJ = fluxVars.volVars(fluxVars.face().j);
-        if (volVarsI.fluidState().saturation(phaseIdx) < 1e-4 ||
-            volVarsJ.fluidState().saturation(phaseIdx) < 1e-4)
+        if (volVarsI.saturation(phaseIdx) < 1e-4 ||
+            volVarsJ.saturation(phaseIdx) < 1e-4)
         {
             return; // phase is not present in one of the finite volumes
         }
@@ -221,8 +221,8 @@ if (!std::isfinite(volumeFlux))
         // integration point by the arithmetic mean of the
         // concentration of the sub-control volumes
         Scalar molarDensityAtIP;
-        molarDensityAtIP = volVarsI.fluidState().molarDensity(phaseIdx);
-        molarDensityAtIP += volVarsJ.fluidState().molarDensity(phaseIdx);
+        molarDensityAtIP = volVarsI.molarDensity(phaseIdx);
+        molarDensityAtIP += volVarsJ.molarDensity(phaseIdx);
         molarDensityAtIP /= 2;
 
         Diffusion::flux(flux, phaseIdx, fluxVars, molarDensityAtIP);
diff --git a/dumux/porousmediumflow/mpnc/implicit/mass/localresidualkinetic.hh b/dumux/porousmediumflow/mpnc/implicit/mass/localresidualkinetic.hh
index e2eb52905d4214743ad6a6429173ab28ea79cf2b..19cf008fae996420c74d3771ae1aab5bd0b4afe9 100644
--- a/dumux/porousmediumflow/mpnc/implicit/mass/localresidualkinetic.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/mass/localresidualkinetic.hh
@@ -195,8 +195,8 @@ public:
         Valgrind::CheckDefined(mu_nPhaseWComp);
 
         const Scalar characteristicLength   = volVars.characteristicLength()  ;
-        const Scalar temperature            = volVars.fluidState().temperature(wPhaseIdx);
-        const Scalar pn                     = volVars.fluidState().pressure(nPhaseIdx);
+        const Scalar temperature            = volVars.temperature(wPhaseIdx);
+        const Scalar pn                     = volVars.pressure(nPhaseIdx);
         const Scalar henry                  = FluidSystem::henry(temperature) ;
         const Scalar gradNinWApprox  = ( mu_wPhaseNComp - mu_nPhaseNCompEquil) / characteristicLength;    // very 2p2c // 1. / henry *
         const Scalar gradWinNApprox  = ( mu_nPhaseWComp - mu_wPhaseWCompEquil) / characteristicLength;    // very 2p2c // 1. / pn *
@@ -205,7 +205,7 @@ public:
         Scalar x[numPhases][numComponents]; // mass fractions in wetting phase
         for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
             for (int compIdx=0; compIdx< numComponents; ++ compIdx){
-                x[phaseIdx][compIdx] = volVars.fluidState().moleFraction(phaseIdx, compIdx);
+                x[phaseIdx][compIdx] = volVars.moleFraction(phaseIdx, compIdx);
             }
         }
         Valgrind::CheckDefined(x);
@@ -221,7 +221,7 @@ public:
 #endif
         Scalar phaseDensity[numPhases];
         for(int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
-            phaseDensity[phaseIdx] = volVars.fluidState().molarDensity(phaseIdx);
+            phaseDensity[phaseIdx] = volVars.molarDensity(phaseIdx);
         }
 
         // diffusion coefficients in wetting phase
diff --git a/dumux/porousmediumflow/mpnc/implicit/mass/vtkwriter.hh b/dumux/porousmediumflow/mpnc/implicit/mass/vtkwriter.hh
index 6fac1ef659c0f573621a8b9c3fff56799b144dfb..62b9f51d47f302902690b80ea25173a970552988 100644
--- a/dumux/porousmediumflow/mpnc/implicit/mass/vtkwriter.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/mass/vtkwriter.hh
@@ -99,7 +99,7 @@ public:
 
             if (fugacityOutput_) {
                 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-                    fugacity_[compIdx][dofIdxGlobal] = volVars.fluidState().fugacity(compIdx);
+                    fugacity_[compIdx][dofIdxGlobal] = volVars.fugacity(compIdx);
                 }
             }
         }
diff --git a/dumux/porousmediumflow/mpnc/implicit/mass/vtkwriterkinetic.hh b/dumux/porousmediumflow/mpnc/implicit/mass/vtkwriterkinetic.hh
index 2295e2d8c0df8225d509818204416eb60ea01997..73b0ba8c9d9c29338f17db1f91a115be691af5f0 100644
--- a/dumux/porousmediumflow/mpnc/implicit/mass/vtkwriterkinetic.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/mass/vtkwriterkinetic.hh
@@ -127,7 +127,7 @@ public:
 //            }
 
             if (deltaPOutput_) {
-                deltaP_[vIdxGlobal] = volVars.fluidState().pressure(nPhaseIdx) - 100000.;
+                deltaP_[vIdxGlobal] = volVars.pressure(nPhaseIdx) - 100000.;
             }
         }
     }
diff --git a/dumux/porousmediumflow/mpnc/implicit/volumevariables.hh b/dumux/porousmediumflow/mpnc/implicit/volumevariables.hh
index a02a636c0ab8ccc285ec909917e0dda1d18f549e..705813c7fb2797f32c6d99bede1f2b3228663933 100644
--- a/dumux/porousmediumflow/mpnc/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/volumevariables.hh
@@ -265,6 +265,114 @@ public:
     const FluidState &fluidState() const
     { return fluidState_; }
 
+    /*!
+     * \brief Returns the saturation of a given phase within
+     *        the control volume in \f$[-]\f$.
+     *
+     * \param phaseIdx The phase index
+     */
+    Scalar saturation(int phaseIdx) const
+    { return fluidState_.saturation(phaseIdx); }
+
+    /*!
+     * \brief Returns the mass fraction of a given component in a
+     *        given phase within the control volume in \f$[-]\f$.
+     *
+     * \param phaseIdx The phase index
+     * \param compIdx The component index
+     */
+    Scalar massFraction(const int phaseIdx, const int compIdx) const
+    { return fluidState_.massFraction(phaseIdx, compIdx); }
+
+    /*!
+     * \brief Returns the mole fraction of a given component in a
+     *        given phase within the control volume in \f$[-]\f$.
+     *
+     * \param phaseIdx The phase index
+     * \param compIdx The component index
+     */
+    Scalar moleFraction(const int phaseIdx, const int compIdx) const
+    { return fluidState_.moleFraction(phaseIdx, compIdx); }
+
+    /*!
+     * \brief Return concentration \f$\mathrm{[mol/m^3]}\f$  of a component in the phase.
+     *
+     * \param compIdx The index of the component
+     */
+    Scalar molarity(const int phaseIdx, int compIdx) const
+    { return fluidState_.molarity(phaseIdx, compIdx); }
+
+    /*!
+     * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ the of the fluid phase.
+     */
+    Scalar molarDensity(const int phaseIdx) const
+    { return fluidState_.molarDensity(phaseIdx);}
+
+    /*!
+     * \brief Return the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
+     *        the control volume.
+     */
+    Scalar pressure(const int phaseIdx) const
+    { return fluidState_.pressure(phaseIdx); }
+
+    /*!
+     * \brief Return density \f$\mathrm{[kg/m^3]}\f$ the of the fluid phase.
+     */
+    Scalar density(const int phaseIdx) const
+    { return fluidState_.density(phaseIdx); }
+
+    /*!
+     * \brief Returns the fluid/solid phase temperature
+     *        in the sub-control volume for the assumption of local thermal
+     *        non-equillibrium where there is more than one energy equation
+     *        and each phase and the matrix can have different temperatures
+     *
+     * \param phaseIdx The local index of the phases
+     */
+    template <class T = TypeTag>
+    typename std::enable_if<(GET_PROP_VALUE(T, NumEnergyEquations) > 1), Scalar>::type temperature(const unsigned int phaseIdx) const
+    {
+        return EnergyVolumeVariables::temperature(phaseIdx);
+    }
+
+     /*!
+     * \brief Returns the fluid/solid phase temperature
+     *        in the sub-control volume for the assumption of local thermal
+     *        equillibrium where there is only one or no energy equation
+     *        and all phases including the  matrix have the same temperature
+     *
+     * \param phaseIdx The local index of the phases
+     */
+    template <class T = TypeTag>
+    typename std::enable_if<(GET_PROP_VALUE(T, NumEnergyEquations) < 2), Scalar>::type temperature(const unsigned int phaseIdx) const
+    {
+        return fluidState_.temperature(phaseIdx);
+    }
+
+    /*!
+     * \brief Return enthalpy \f$\mathrm{[kg/m^3]}\f$ the of the fluid phase.
+     */
+    Scalar enthalpy(const int phaseIdx) const
+    { return fluidState_.enthalpy(phaseIdx); }
+
+    /*!
+     * \brief Return internal energy \f$\mathrm{[kg/m^3]}\f$ the of the fluid phase.
+     */
+    Scalar internalEnergy(const int phaseIdx) const
+    { return fluidState_.internalEnergy(phaseIdx); }
+
+    /*!
+     * \brief Return fugacity \f$\mathrm{[kg/m^3]}\f$ the of the component.
+     */
+    Scalar fugacity(const int compIdx) const
+    { return fluidState_.fugacity(compIdx); }
+
+    /*!
+     * \brief Return average molar mass \f$\mathrm{[kg/m^3]}\f$ the of the phase.
+     */
+    Scalar averageMolarMass(const int phaseIdx) const
+    { return fluidState_.averageMolarMass(phaseIdx); }
+
     /*!
      * \brief Returns the effective mobility of a given phase within
      *        the control volume.
diff --git a/dumux/porousmediumflow/mpnc/implicit/vtkwritercommon.hh b/dumux/porousmediumflow/mpnc/implicit/vtkwritercommon.hh
index b45128e4832cb6d8276b34e0640f05627e4a8efc..58b68332ebe5c0b8cd1dad1477d0f48d275be3b4 100644
--- a/dumux/porousmediumflow/mpnc/implicit/vtkwritercommon.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/vtkwritercommon.hh
@@ -151,15 +151,15 @@ public:
             if (boundaryTypesOutput_) boundaryTypes_[dofIdxGlobal] = tmp;
 
             for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-                if (saturationOutput_) saturation_[phaseIdx][dofIdxGlobal] = volVars.fluidState().saturation(phaseIdx);
-                if (pressureOutput_) pressure_[phaseIdx][dofIdxGlobal] = volVars.fluidState().pressure(phaseIdx);
-                if (densityOutput_) density_[phaseIdx][dofIdxGlobal] = volVars.fluidState().density(phaseIdx);
+                if (saturationOutput_) saturation_[phaseIdx][dofIdxGlobal] = volVars.saturation(phaseIdx);
+                if (pressureOutput_) pressure_[phaseIdx][dofIdxGlobal] = volVars.pressure(phaseIdx);
+                if (densityOutput_) density_[phaseIdx][dofIdxGlobal] = volVars.density(phaseIdx);
                 if (mobilityOutput_) mobility_[phaseIdx][dofIdxGlobal] = volVars.mobility(phaseIdx);
-                if (averageMolarMassOutput_) averageMolarMass_[phaseIdx][dofIdxGlobal] = volVars.fluidState().averageMolarMass(phaseIdx);
+                if (averageMolarMassOutput_) averageMolarMass_[phaseIdx][dofIdxGlobal] = volVars.averageMolarMass(phaseIdx);
                 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-                    if (moleFracOutput_) moleFrac_[phaseIdx][compIdx][dofIdxGlobal] = volVars.fluidState().moleFraction(phaseIdx, compIdx);
-                    if (massFracOutput_) massFrac_[phaseIdx][compIdx][dofIdxGlobal] = volVars.fluidState().massFraction(phaseIdx, compIdx);
-                    if (molarityOutput_) molarity_[phaseIdx][compIdx][dofIdxGlobal] = volVars.fluidState().molarity(phaseIdx, compIdx);
+                    if (moleFracOutput_) moleFrac_[phaseIdx][compIdx][dofIdxGlobal] = volVars.moleFraction(phaseIdx, compIdx);
+                    if (massFracOutput_) massFrac_[phaseIdx][compIdx][dofIdxGlobal] = volVars.massFraction(phaseIdx, compIdx);
+                    if (molarityOutput_) molarity_[phaseIdx][compIdx][dofIdxGlobal] = volVars.molarity(phaseIdx, compIdx);
                 }
             }
         }
diff --git a/test/porousmediumflow/mpnc/implicit/combustionproblem1c.hh b/test/porousmediumflow/mpnc/implicit/combustionproblem1c.hh
index 92ad6858d14a4f9008402c23ea4da88c68035798..c102e016524265cc15c5af41364e403adade0e7f 100644
--- a/test/porousmediumflow/mpnc/implicit/combustionproblem1c.hh
+++ b/test/porousmediumflow/mpnc/implicit/combustionproblem1c.hh
@@ -478,11 +478,11 @@ public:
 
         FluidState fluidState;
 
-        const Scalar pn = elemVolVars[scvIdx].fluidState().pressure(nPhaseIdx);
-        const Scalar pw = elemVolVars[scvIdx].fluidState().pressure(wPhaseIdx);
+        const Scalar pn = elemVolVars[scvIdx].pressure(nPhaseIdx);
+        const Scalar pw = elemVolVars[scvIdx].pressure(wPhaseIdx);
 
-        const Scalar Tn = elemVolVars[scvIdx].fluidState().temperature(nPhaseIdx);
-        const Scalar Tw = elemVolVars[scvIdx].fluidState().temperature(wPhaseIdx);
+        const Scalar Tn = elemVolVars[scvIdx].temperature(nPhaseIdx);
+        const Scalar Tw = elemVolVars[scvIdx].temperature(wPhaseIdx);
 
         fluidState.setPressure(nPhaseIdx, pn);
         fluidState.setPressure(wPhaseIdx, pw);