diff --git a/CHANGELOG b/CHANGELOG
index 39a81c7d526adfcb396d1d1430f4d55f2f0ed06b..c9b6b66730bd475f80ef494bc507f320bcd565f1 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -8,6 +8,10 @@ Differences Between DuMuX 2.1 and DuMuX 2.2
     has to be replaced by 
     SET...PROP(..., SpatialParams ...
 
+  - TwoPIndices, TwoPNIIndices, and RichardsIndices additionally get TypeTag 
+    as template parameter. If the Indices are not obtained via the property, 
+    this has to be adapted.
+
 * Deprecated CLASSES/FILES, to be removed after 2.2:
   - Model specific base box problems: The common functionality has been collected 
     in PorousMediaBoxProblem in dumux/boxmodels/common/porousmediaboxproblem.hh. 
diff --git a/dumux/boxmodels/2p/2pindices.hh b/dumux/boxmodels/2p/2pindices.hh
index e09647895939c5befe2276dd70e90c4bee0871b8..78d1157850300601aa33034fe6f7e43162ceed2a 100644
--- a/dumux/boxmodels/2p/2pindices.hh
+++ b/dumux/boxmodels/2p/2pindices.hh
@@ -38,26 +38,36 @@ namespace Dumux
  * \ingroup BoxIndices
  * \brief The common indices for the isothermal two-phase model.
  */
-struct TwoPCommonIndices
+
+struct TwoPFormulation
 {
-    // Formulations
     static const int pwSn = 0; //!< Pw and Sn as primary variables
     static const int pnSw = 1; //!< Pn and Sw as primary variables
+};
+
+template <class TypeTag>
+struct TwoPCommonIndices
+{
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
 
     // Phase indices
-    static const int wPhaseIdx = 0; //!< Index of the wetting phase in a phase vector
-    static const int nPhaseIdx = 1; //!< Index of the non-wetting phase in a phase vector
+    static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< Index of the wetting phase
+    static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< Index of the non-wetting phase
 };
 
 /*!
  * \brief The indices for the \f$p_w-S_n\f$ formulation of the
  *        isothermal two-phase model.
  *
+ * \tparam TypeTag The problem type tag
  * \tparam formulation The formulation, either pwSn or pnSw
  * \tparam PVOffset The first index in a primary variable vector.
  */
-template <int formulation = TwoPCommonIndices::pwSn, int PVOffset = 0>
-struct TwoPIndices : public TwoPCommonIndices
+template <class TypeTag, 
+          int formulation = TwoPFormulation::pwSn, 
+          int PVOffset = 0>
+struct TwoPIndices 
+: public TwoPCommonIndices<TypeTag>, TwoPFormulation
 {
     // Primary variable indices
     static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
@@ -78,9 +88,9 @@ struct TwoPIndices : public TwoPCommonIndices
  *
  * \tparam PVOffset The first index in a primary variable vector.
  */
-template <int PVOffset>
-struct TwoPIndices<TwoPCommonIndices::pnSw, PVOffset>
-    : public TwoPCommonIndices
+template <class TypeTag, int PVOffset>
+struct TwoPIndices<TypeTag, TwoPFormulation::pnSw, PVOffset>
+: public TwoPCommonIndices<TypeTag>, TwoPFormulation
 {
     // Primary variable indices
     static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
diff --git a/dumux/boxmodels/2p/2ppropertydefaults.hh b/dumux/boxmodels/2p/2ppropertydefaults.hh
index 6f6b2a3e5789ae2bd4aa0eae064e8e45bab40785..2db607f0469e6a9513131aad28c785181aa98555 100644
--- a/dumux/boxmodels/2p/2ppropertydefaults.hh
+++ b/dumux/boxmodels/2p/2ppropertydefaults.hh
@@ -59,7 +59,7 @@ SET_INT_PROP(BoxTwoP, NumPhases, 2); //!< The number of phases in the 2p model i
 //! Set the default formulation to pWsN
 SET_INT_PROP(BoxTwoP,
              Formulation,
-             TwoPCommonIndices::pwSn);
+             TwoPFormulation::pwSn);
 
 //! Use the 2p local jacobian operator for the 2p model
 SET_TYPE_PROP(BoxTwoP,
@@ -81,7 +81,7 @@ SET_SCALAR_PROP(BoxTwoP, MassUpwindWeight, 1.0);
 //! The indices required by the isothermal 2p model
 SET_TYPE_PROP(BoxTwoP,
               Indices,
-              TwoPIndices<GET_PROP_VALUE(TypeTag, Formulation), 0>);
+              TwoPIndices<TypeTag, GET_PROP_VALUE(TypeTag, Formulation), 0>);
 
 //! DEPRECATED TwoPIndices property
 SET_TYPE_PROP(BoxTwoP, TwoPIndices, typename GET_PROP_TYPE(TypeTag, Indices));
diff --git a/dumux/boxmodels/2p2c/2p2cindices.hh b/dumux/boxmodels/2p2c/2p2cindices.hh
index a9e1cb54a36e8a4ffb44025fe45639e9b2984f36..9f16c4b01091cc2545a94396be98cd089b3a2085 100644
--- a/dumux/boxmodels/2p2c/2p2cindices.hh
+++ b/dumux/boxmodels/2p2c/2p2cindices.hh
@@ -75,8 +75,8 @@ public:
     static const int DUMUX_DEPRECATED_MSG("use nPhaseIdx instead") gPhaseIdx = nPhaseIdx; //!< Index of the gas phase
 
     // Component indices
-    static const int wCompIdx = 0; //!< Index of the primary component of the wetting phase
-    static const int nCompIdx = 1; //!< Index of the primary component of the non-wetting phase
+    static const int wCompIdx = FluidSystem::wCompIdx; //!< Index of the primary component of the wetting phase
+    static const int nCompIdx = FluidSystem::nCompIdx; //!< Index of the primary component of the non-wetting phase
 
     static const int DUMUX_DEPRECATED_MSG("use wCompIdx instead") lCompIdx = wCompIdx; //!< Index of the liquid's primary component
     static const int DUMUX_DEPRECATED_MSG("use nCompIdx instead") gCompIdx = nCompIdx; //!< Index of the gas' primary component
@@ -127,8 +127,8 @@ public:
     static const int DUMUX_DEPRECATED_MSG("use nPhaseIdx instead") gPhaseIdx = nPhaseIdx; //!< Index of the gas phase
 
     // Component indices
-    static const int wCompIdx = 0; //!< Index of the primary component of the wetting phase
-    static const int nCompIdx = 1; //!< Index of the primary component of the non-wetting phase
+    static const int wCompIdx = FluidSystem::wCompIdx; //!< Index of the primary component of the wetting phase
+    static const int nCompIdx = FluidSystem::nCompIdx; //!< Index of the primary component of the non-wetting phase
 
     static const int DUMUX_DEPRECATED_MSG("use wCompIdx instead") lCompIdx = wCompIdx; //!< Index of the liquid's primary component
     static const int DUMUX_DEPRECATED_MSG("use nCompIdx instead") gCompIdx = nCompIdx; //!< Index of the gas' primary component
diff --git a/dumux/boxmodels/2pni/2pniindices.hh b/dumux/boxmodels/2pni/2pniindices.hh
index b5911e8df70f7d832ed323f3862a0020f4ed3d2f..76fe96066214027ec563e2051d99c7c9765db759 100644
--- a/dumux/boxmodels/2pni/2pniindices.hh
+++ b/dumux/boxmodels/2pni/2pniindices.hh
@@ -40,15 +40,12 @@ namespace Dumux
  * \ingroup BoxIndices
  * \brief Enumerations for the non-isothermal two-phase model
  */
-template <int PVOffset = 0>
-class TwoPNIIndices : public TwoPIndices<PVOffset>
+template <class TypeTag, int PVOffset = 0>
+class TwoPNIIndices : public TwoPIndices<TypeTag, PVOffset>
 {
 public:
     static const int temperatureIdx = PVOffset + 2; //! The primary variable index for temperature
     static const int energyEqIdx = PVOffset + 2; //! The equation index of the energy equation
-    // Phase indices
-       static const int wPhaseIdx = 0; //!< Index of the wetting phase in a phase vector
-       static const int nPhaseIdx = 1; //!< Index of the non-wetting phase in a phase vector
 };
 // \}
 }
diff --git a/dumux/boxmodels/2pni/2pnipropertydefaults.hh b/dumux/boxmodels/2pni/2pnipropertydefaults.hh
index a2cbfa68fa98726fb7194c3353143422a01b9a9c..9f5c36f0a8ba859eb354d4f210410dd90e49e77b 100644
--- a/dumux/boxmodels/2pni/2pnipropertydefaults.hh
+++ b/dumux/boxmodels/2pni/2pnipropertydefaults.hh
@@ -68,7 +68,7 @@ SET_TYPE_PROP(BoxTwoPNI, VolumeVariables, TwoPNIVolumeVariables<TypeTag>);
 SET_TYPE_PROP(BoxTwoPNI, FluxVariables, TwoPNIFluxVariables<TypeTag>);
 
 //! The indices required by the non-isothermal two-phase model
-SET_TYPE_PROP(BoxTwoPNI, Indices, TwoPNIIndices<0>);
+SET_TYPE_PROP(BoxTwoPNI, Indices, TwoPNIIndices<TypeTag, 0>);
 
 //! DEPRECATED TwoPIndices and TwoPNIIndices properties
 SET_TYPE_PROP(BoxTwoPNI, TwoPIndices, typename GET_PROP_TYPE(TypeTag, Indices));
diff --git a/dumux/boxmodels/3p3c/3p3cfluxvariables.hh b/dumux/boxmodels/3p3c/3p3cfluxvariables.hh
index 6535827ad730bd1565d4d5bbe3f4bde511f1546f..bb551a9b45b93ee9018c9ca33b0e1450dca516d8 100644
--- a/dumux/boxmodels/3p3c/3p3cfluxvariables.hh
+++ b/dumux/boxmodels/3p3c/3p3cfluxvariables.hh
@@ -83,9 +83,9 @@ class ThreePThreeCFluxVariables
         nPhaseIdx = Indices::nPhaseIdx,
         gPhaseIdx = Indices::gPhaseIdx,
 
-        comp0Idx = Indices::comp0Idx,
-        comp1Idx = Indices::comp1Idx,
-        comp2Idx = Indices::comp2Idx
+        wCompIdx = Indices::wCompIdx,
+        nCompIdx = Indices::nCompIdx,
+        gCompIdx = Indices::gCompIdx
     };
 
 public:
@@ -112,12 +112,12 @@ public:
             density_[phaseIdx] = Scalar(0);
             molarDensity_[phaseIdx] = Scalar(0);
             potentialGrad_[phaseIdx] = Scalar(0);
-            massFractionComp0Grad_[phaseIdx] = Scalar(0);
-            massFractionComp1Grad_[phaseIdx] = Scalar(0);
-            massFractionComp2Grad_[phaseIdx] = Scalar(0);
-            moleFractionComp0Grad_[phaseIdx] = Scalar(0);
-            moleFractionComp1Grad_[phaseIdx] = Scalar(0);
-            moleFractionComp2Grad_[phaseIdx] = Scalar(0);
+            massFractionCompWGrad_[phaseIdx] = Scalar(0);
+            massFractionCompNGrad_[phaseIdx] = Scalar(0);
+            massFractionCompGGrad_[phaseIdx] = Scalar(0);
+            moleFractionCompWGrad_[phaseIdx] = Scalar(0);
+            moleFractionCompNGrad_[phaseIdx] = Scalar(0);
+            moleFractionCompGGrad_[phaseIdx] = Scalar(0);
         }
 
         calculateGradients_(problem, element, elemVolVars);
@@ -154,78 +154,78 @@ private:
             // the concentration gradient of the components
             // component in the phases
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(wPhaseIdx, comp0Idx);
-            massFractionComp0Grad_[wPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(wPhaseIdx, wCompIdx);
+            massFractionCompWGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(nPhaseIdx, comp0Idx);
-            massFractionComp0Grad_[nPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(nPhaseIdx, wCompIdx);
+            massFractionCompWGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(gPhaseIdx, comp0Idx);
-            massFractionComp0Grad_[gPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(gPhaseIdx, wCompIdx);
+            massFractionCompWGrad_[gPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(wPhaseIdx, comp1Idx);
-            massFractionComp1Grad_[wPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(wPhaseIdx, nCompIdx);
+            massFractionCompNGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(nPhaseIdx, comp1Idx);
-            massFractionComp1Grad_[nPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(nPhaseIdx, nCompIdx);
+            massFractionCompNGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(gPhaseIdx, comp1Idx);
-            massFractionComp1Grad_[gPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(gPhaseIdx, nCompIdx);
+            massFractionCompNGrad_[gPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(wPhaseIdx, comp2Idx);
-            massFractionComp2Grad_[wPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(wPhaseIdx, gCompIdx);
+            massFractionCompGGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(nPhaseIdx, comp2Idx);
-            massFractionComp2Grad_[nPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(nPhaseIdx, gCompIdx);
+            massFractionCompGGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(gPhaseIdx, comp2Idx);
-            massFractionComp2Grad_[gPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(gPhaseIdx, gCompIdx);
+            massFractionCompGGrad_[gPhaseIdx] += tmp;
 
             // the molar concentration gradients of the components
             // in the phases
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(wPhaseIdx, comp0Idx);
-            moleFractionComp0Grad_[wPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(wPhaseIdx, wCompIdx);
+            moleFractionCompWGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, comp0Idx);
-            moleFractionComp0Grad_[nPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
+            moleFractionCompWGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(gPhaseIdx, comp0Idx);
-            moleFractionComp0Grad_[gPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(gPhaseIdx, wCompIdx);
+            moleFractionCompWGrad_[gPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(wPhaseIdx, comp1Idx);
-            moleFractionComp1Grad_[wPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(wPhaseIdx, nCompIdx);
+            moleFractionCompNGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, comp1Idx);
-            moleFractionComp1Grad_[nPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, nCompIdx);
+            moleFractionCompNGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(gPhaseIdx, comp1Idx);
-            moleFractionComp1Grad_[gPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(gPhaseIdx, nCompIdx);
+            moleFractionCompNGrad_[gPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(wPhaseIdx, comp2Idx);
-            moleFractionComp2Grad_[wPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(wPhaseIdx, gCompIdx);
+            moleFractionCompGGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, comp2Idx);
-            moleFractionComp2Grad_[nPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, gCompIdx);
+            moleFractionCompGGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(gPhaseIdx, comp2Idx);
-            moleFractionComp2Grad_[gPhaseIdx] += tmp;
+            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(gPhaseIdx, gCompIdx);
+            moleFractionCompGGrad_[gPhaseIdx] += tmp;
         }
 
         // correct the pressure gradients by the hydrostatic
@@ -329,9 +329,9 @@ private:
             if (volVarsI.saturation(phaseIdx) <= 0 ||
                 volVarsJ.saturation(phaseIdx) <= 0)
             {
-                porousDiffCoeff_[phaseIdx][comp0Idx] = 0.0;
-                porousDiffCoeff_[phaseIdx][comp1Idx] = 0.0;
-                porousDiffCoeff_[phaseIdx][comp2Idx] = 0.0;
+                porousDiffCoeff_[phaseIdx][wCompIdx] = 0.0;
+                porousDiffCoeff_[phaseIdx][nCompIdx] = 0.0;
+                porousDiffCoeff_[phaseIdx][gCompIdx] = 0.0;
                 continue;
             }
 
@@ -347,12 +347,12 @@ private:
             // Diffusion coefficient in the porous medium
 
             // -> harmonic mean
-            porousDiffCoeff_[phaseIdx][comp0Idx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][comp0Idx],
-                                                                volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][comp0Idx]);
-            porousDiffCoeff_[phaseIdx][comp1Idx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][comp1Idx],
-                                                                volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][comp1Idx]);
-            porousDiffCoeff_[phaseIdx][comp2Idx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][comp2Idx],
-                                                                volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][comp2Idx]);
+            porousDiffCoeff_[phaseIdx][wCompIdx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][wCompIdx],
+                                                                volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][wCompIdx]);
+            porousDiffCoeff_[phaseIdx][nCompIdx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][nCompIdx],
+                                                                volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][nCompIdx]);
+            porousDiffCoeff_[phaseIdx][gCompIdx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * diffusionCoefficientMatrix_i[phaseIdx][gCompIdx],
+                                                                volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * diffusionCoefficientMatrix_j[phaseIdx][gCompIdx]);
 
         }
     }
@@ -434,51 +434,51 @@ public:
     /*!
      * \brief The mass fraction gradients of the components in a phase.
      */
-	DUMUX_DEPRECATED_MSG("use massFractionComp0Grad instead")
+	DUMUX_DEPRECATED_MSG("use massFractionCompWGrad instead")
     const DimVector &wConcentrationGrad(int phaseIdx) const
-    { massFractionComp0Grad(phaseIdx); };
+    { massFractionCompWGrad(phaseIdx); };
 
-	const DimVector &massFractionComp0Grad(int phaseIdx) const
-	{return massFractionComp0Grad_[phaseIdx];}
+	const DimVector &massFractionCompWGrad(int phaseIdx) const
+	{return massFractionCompWGrad_[phaseIdx];}
 
-	DUMUX_DEPRECATED_MSG("use massFractionComp1Grad instead")
+	DUMUX_DEPRECATED_MSG("use massFractionCompNGrad instead")
     const DimVector &cConcentrationGrad(int phaseIdx) const
-    { massFractionComp1Grad(phaseIdx); };
+    { massFractionCompNGrad(phaseIdx); };
 
-	const DimVector &massFractionComp1Grad(int phaseIdx) const
-	{ return massFractionComp1Grad_[phaseIdx]; };
+	const DimVector &massFractionCompNGrad(int phaseIdx) const
+	{ return massFractionCompNGrad_[phaseIdx]; };
 
-	DUMUX_DEPRECATED_MSG("use massFractionComp2Grad instead")
+	DUMUX_DEPRECATED_MSG("use massFractionCompGGrad instead")
     const DimVector &aConcentrationGrad(int phaseIdx) const
-    {  massFractionComp2Grad(phaseIdx); };
+    {  massFractionCompGGrad(phaseIdx); };
 
-	const DimVector &massFractionComp2Grad(int phaseIdx) const
-	{ return massFractionComp2Grad_[phaseIdx]; };
+	const DimVector &massFractionCompGGrad(int phaseIdx) const
+	{ return massFractionCompGGrad_[phaseIdx]; };
 
 
     /*!
      * \brief The molar concentration gradients of the components in a phase.
      */
-	DUMUX_DEPRECATED_MSG("use moleFractionComp0Grad instead")
+	DUMUX_DEPRECATED_MSG("use moleFractionCompWGrad instead")
     const DimVector &molarWConcGrad(int phaseIdx) const
-    {  moleFractionComp0Grad(phaseIdx); };
+    {  moleFractionCompWGrad(phaseIdx); };
 
-	const DimVector &moleFractionComp0Grad(int phaseIdx) const
-	{ return moleFractionComp0Grad_[phaseIdx]; };
+	const DimVector &moleFractionCompWGrad(int phaseIdx) const
+	{ return moleFractionCompWGrad_[phaseIdx]; };
 
-	DUMUX_DEPRECATED_MSG("use moleFractionComp1Grad instead")
+	DUMUX_DEPRECATED_MSG("use moleFractionCompNGrad instead")
     const DimVector &molarCConcGrad(int phaseIdx) const
-    { moleFractionComp1Grad(phaseIdx); };
+    { moleFractionCompNGrad(phaseIdx); };
 
-	const DimVector &moleFractionComp1Grad(int phaseIdx) const
-	{ return moleFractionComp1Grad_[phaseIdx]; };
+	const DimVector &moleFractionCompNGrad(int phaseIdx) const
+	{ return moleFractionCompNGrad_[phaseIdx]; };
 
-	DUMUX_DEPRECATED_MSG("use moleFractionComp2Grad instead")
+	DUMUX_DEPRECATED_MSG("use moleFractionCompGGrad instead")
     const DimVector &molarAConcGrad(int phaseIdx) const
-    { moleFractionComp2Grad(phaseIdx); };
+    { moleFractionCompGGrad(phaseIdx); };
 
-	const DimVector &moleFractionComp2Grad(int phaseIdx) const
-    { return moleFractionComp2Grad_[phaseIdx]; };
+	const DimVector &moleFractionCompGGrad(int phaseIdx) const
+    { return moleFractionCompGGrad_[phaseIdx]; };
 
     const SCVFace &face() const
     {
@@ -495,12 +495,12 @@ protected:
 
     // gradients
     DimVector potentialGrad_[numPhases];
-    DimVector massFractionComp0Grad_[numPhases];
-    DimVector massFractionComp1Grad_[numPhases];
-    DimVector massFractionComp2Grad_[numPhases];
-    DimVector moleFractionComp0Grad_[numPhases];
-    DimVector moleFractionComp1Grad_[numPhases];
-    DimVector moleFractionComp2Grad_[numPhases];
+    DimVector massFractionCompWGrad_[numPhases];
+    DimVector massFractionCompNGrad_[numPhases];
+    DimVector massFractionCompGGrad_[numPhases];
+    DimVector moleFractionCompWGrad_[numPhases];
+    DimVector moleFractionCompNGrad_[numPhases];
+    DimVector moleFractionCompGGrad_[numPhases];
 
     // density of each face at the integration point
     Scalar density_[numPhases], molarDensity_[numPhases];
diff --git a/dumux/boxmodels/3p3c/3p3cindices.hh b/dumux/boxmodels/3p3c/3p3cindices.hh
index 4b88ac95206e27a2a3448fb22088a2849d0750f3..1d5869fdb67a68f450e9de84a77eb2beed135634 100644
--- a/dumux/boxmodels/3p3c/3p3cindices.hh
+++ b/dumux/boxmodels/3p3c/3p3cindices.hh
@@ -48,17 +48,19 @@ class ThreePThreeCIndices
 
 public:
     // Phase indices
-    static const int wPhaseIdx = FluidSystem::wPhaseIdx; // !< DEPRECATED index of the water phase
-    static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< DEPRECATED index of the NAPL phase
-    static const int gPhaseIdx = FluidSystem::gPhaseIdx; //< DEPRECATED index of the gas phase
+    static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< index of the wetting liquid phase
+    static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< index of the nonwetting liquid phase
+    static const int gPhaseIdx = FluidSystem::gPhaseIdx; //!< index of the gas phase
 
-    // Component indices
-    static const int wCompIdx = 0; //!< Index of the water component
-    static const int cCompIdx = 1; //!< Index of the NAPL component
-    static const int aCompIdx = 2; //!< Index of the air component
-    static const int comp0Idx = 0; //!< Index of the water component
-    static const int comp1Idx = 1; //!< Index of the NAPL component
-    static const int comp2Idx = 2; //!< Index of the gas component
+    // Component indices to indicate the main component 
+    // of the corresponding phase at atmospheric pressure 1 bar 
+    // and room temperature 20°C:
+    static const int wCompIdx = FluidSystem::wCompIdx; 
+    static const int nCompIdx = FluidSystem::nCompIdx; 
+    static const int gCompIdx = FluidSystem::gCompIdx; 
+
+    static const int cCompIdx = nCompIdx; //!< DEPRECATED index of the NAPL component
+    static const int aCompIdx = gCompIdx; //!< DEPRECATED index of the air component
 
     // present phases (-> 'pseudo' primary variable)
     static const int threePhases = 1; //!< All three phases are present
@@ -78,10 +80,9 @@ public:
     static const int SOrX2Idx = switch2Idx; //!< Index of the either the saturation of the gas phase or the mass fraction secondary component if a phase is not present
 
     // equation indices
-
-    static const int conti0EqIdx = PVOffset    + comp0Idx; //!< Index of the mass conservation equation for the water component
-    static const int conti1EqIdx = conti0EqIdx + comp1Idx; //!< Index of the mass conservation equation for the contaminant component
-    static const int conti2EqIdx = conti0EqIdx + comp2Idx; //!< Index of the mass conservation equation for the gas component
+    static const int conti0EqIdx = PVOffset    + wCompIdx; //!< Index of the mass conservation equation for the water component
+    static const int conti1EqIdx = conti0EqIdx + nCompIdx; //!< Index of the mass conservation equation for the contaminant component
+    static const int conti2EqIdx = conti0EqIdx + gCompIdx; //!< Index of the mass conservation equation for the gas component
 
     static const int contiWEqIdx = conti0EqIdx + wCompIdx; //!< DEPRECATED index of the mass conservation equation for the water component
     static const int contiCEqIdx = conti0EqIdx + cCompIdx; //!< DEPRECATED index of the mass conservation equation for the contaminant component
diff --git a/dumux/boxmodels/3p3c/3p3clocalresidual.hh b/dumux/boxmodels/3p3c/3p3clocalresidual.hh
index 72c8570584782b23fcc2e4a6b3baf44d0600a099..8c83931e0dbc3712211cfeb8f7c701ef24e33d04 100644
--- a/dumux/boxmodels/3p3c/3p3clocalresidual.hh
+++ b/dumux/boxmodels/3p3c/3p3clocalresidual.hh
@@ -76,9 +76,9 @@ protected:
         nPhaseIdx = Indices::nPhaseIdx,
         gPhaseIdx = Indices::gPhaseIdx,
 
-        comp0Idx = Indices::comp0Idx,
-        comp1Idx = Indices::comp1Idx,
-        comp2Idx = Indices::comp2Idx
+        wCompIdx = Indices::wCompIdx,
+        nCompIdx = Indices::nCompIdx,
+        gCompIdx = Indices::gCompIdx
     };
 
     typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
@@ -205,51 +205,39 @@ public:
     {
         // TODO: reference!?  Dune::FieldMatrix<Scalar, numPhases, numComponents> averagedPorousDiffCoeffMatrix = fluxVars.porousDiffCoeff();
         // add diffusive flux of gas component in liquid phase
-        Scalar tmp = - fluxVars.porousDiffCoeff()[wPhaseIdx][comp2Idx] * fluxVars.molarDensity(wPhaseIdx);
-        tmp *= (fluxVars.moleFractionComp2Grad(wPhaseIdx) * fluxVars.face().normal);
-        Scalar j2W = tmp;
-
-        tmp = - fluxVars.porousDiffCoeff()[wPhaseIdx][comp1Idx] * fluxVars.molarDensity(wPhaseIdx);
-        tmp *= (fluxVars.moleFractionComp1Grad(wPhaseIdx) * fluxVars.face().normal);
-        Scalar j1W = tmp;
-
-        Scalar j0W = -(j2W+j1W);
-
-        tmp = - fluxVars.porousDiffCoeff()[gPhaseIdx][comp0Idx] * fluxVars.molarDensity(gPhaseIdx);
-        tmp *= (fluxVars.moleFractionComp0Grad(gPhaseIdx) * fluxVars.face().normal);
-        Scalar j0G = tmp;
-
-        tmp = - fluxVars.porousDiffCoeff()[gPhaseIdx][comp1Idx] * fluxVars.molarDensity(gPhaseIdx);
-        tmp *= (fluxVars.moleFractionComp1Grad(gPhaseIdx) * fluxVars.face().normal);
-        Scalar j1G = tmp;
-
-        Scalar j2G = -(j0G+j1G);
-
-        tmp = - fluxVars.porousDiffCoeff()[nPhaseIdx][comp0Idx] * fluxVars.molarDensity(nPhaseIdx);
-        tmp *= (fluxVars.moleFractionComp0Grad(nPhaseIdx) * fluxVars.face().normal);
-        Scalar j0N = tmp;
-
-        tmp = - fluxVars.porousDiffCoeff()[nPhaseIdx][comp2Idx] * fluxVars.molarDensity(nPhaseIdx);
-        tmp *= (fluxVars.moleFractionComp2Grad(nPhaseIdx) * fluxVars.face().normal);
-        Scalar j2N = tmp;
-
-        Scalar j1N = -(j2N+j0N);
-
-        /*
-          j0W = 0;
-          j0G = 0;
-          j0N = 0;
-          j1W = 0;
-          j1G = 0;
-          j1N = 0;
-          j2W = 0;
-          j2G = 0;
-          j2N = 0;
-        */
-
-        flux[conti0EqIdx] += j0W+j0G+j0N;
-        flux[conti1EqIdx] += j1W+j1G+j1N;
-        flux[conti2EqIdx] += j2W+j2G+j2N;
+        Scalar tmp = - fluxVars.porousDiffCoeff()[wPhaseIdx][gCompIdx] * fluxVars.molarDensity(wPhaseIdx);
+        tmp *= (fluxVars.moleFractionCompGGrad(wPhaseIdx) * fluxVars.face().normal);
+        Scalar jGW = tmp;
+
+        tmp = - fluxVars.porousDiffCoeff()[wPhaseIdx][nCompIdx] * fluxVars.molarDensity(wPhaseIdx);
+        tmp *= (fluxVars.moleFractionCompNGrad(wPhaseIdx) * fluxVars.face().normal);
+        Scalar jNW = tmp;
+
+        Scalar jWW = -(jGW+jNW);
+
+        tmp = - fluxVars.porousDiffCoeff()[gPhaseIdx][wCompIdx] * fluxVars.molarDensity(gPhaseIdx);
+        tmp *= (fluxVars.moleFractionCompWGrad(gPhaseIdx) * fluxVars.face().normal);
+        Scalar jWG = tmp;
+
+        tmp = - fluxVars.porousDiffCoeff()[gPhaseIdx][nCompIdx] * fluxVars.molarDensity(gPhaseIdx);
+        tmp *= (fluxVars.moleFractionCompNGrad(gPhaseIdx) * fluxVars.face().normal);
+        Scalar jNG = tmp;
+
+        Scalar jGG = -(jWG+jNG);
+
+        tmp = - fluxVars.porousDiffCoeff()[nPhaseIdx][wCompIdx] * fluxVars.molarDensity(nPhaseIdx);
+        tmp *= (fluxVars.moleFractionCompWGrad(nPhaseIdx) * fluxVars.face().normal);
+        Scalar jWN = tmp;
+
+        tmp = - fluxVars.porousDiffCoeff()[nPhaseIdx][gCompIdx] * fluxVars.molarDensity(nPhaseIdx);
+        tmp *= (fluxVars.moleFractionCompGGrad(nPhaseIdx) * fluxVars.face().normal);
+        Scalar jGN = tmp;
+
+        Scalar jNN = -(jGN+jWN);
+
+        flux[conti0EqIdx] += jWW+jWG+jWN;
+        flux[conti1EqIdx] += jNW+jNG+jNN;
+        flux[conti2EqIdx] += jGW+jGG+jGN;
     }
 
     /*!
diff --git a/dumux/boxmodels/3p3c/3p3cmodel.hh b/dumux/boxmodels/3p3c/3p3cmodel.hh
index 91ac7f2f6b30bd4a32426a2deb10f0b483950bde..75567cb65eb9eaa20437b9402c1e2316b8ad29a6 100644
--- a/dumux/boxmodels/3p3c/3p3cmodel.hh
+++ b/dumux/boxmodels/3p3c/3p3cmodel.hh
@@ -126,9 +126,9 @@ class ThreePThreeCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
         nPhaseIdx = Indices::nPhaseIdx,
         gPhaseIdx = Indices::gPhaseIdx,
 
-        comp0Idx = Indices::comp0Idx,
-        comp1Idx = Indices::comp1Idx,
-        comp2Idx = Indices::comp2Idx,
+        wCompIdx = Indices::wCompIdx,
+        nCompIdx = Indices::nCompIdx,
+        gCompIdx = Indices::gCompIdx,
 
         threePhases = Indices::threePhases,
         wPhaseOnly  = Indices::wPhaseOnly,
@@ -559,7 +559,7 @@ protected:
                 newPhasePresence = wnPhaseOnly;
 
                 globalSol[globalIdx][switch1Idx]
-                    = volVars.fluidState().moleFraction(wPhaseIdx, comp2Idx);
+                    = volVars.fluidState().moleFraction(wPhaseIdx, gCompIdx);
             }
             else if (volVars.saturation(wPhaseIdx) <= Smin)
             {
@@ -571,7 +571,7 @@ protected:
                 newPhasePresence = gnPhaseOnly;
 
                 globalSol[globalIdx][switch1Idx]
-                    = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
+                    = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
             }
             else if (volVars.saturation(nPhaseIdx) <= Smin)
             {
@@ -583,91 +583,91 @@ protected:
                 newPhasePresence = wgPhaseOnly;
 
                 globalSol[globalIdx][switch2Idx]
-                    = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
+                    = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
             }
         }
         else if (phasePresence == wPhaseOnly)
         {
             bool gasFlag = 0;
-            bool NAPLFlag = 0;
+            bool nonwettingFlag = 0;
             // calculate fractions of the partial pressures in the
             // hypothetical gas phase
-            Scalar x0g = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
-            Scalar x2g = volVars.fluidState().moleFraction(gPhaseIdx, comp2Idx);
-            Scalar x1g = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
+            Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+            Scalar xgg = volVars.fluidState().moleFraction(gPhaseIdx, gCompIdx);
+            Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
             /* take care:
-               for x2g in case wPhaseOnly we compute x2g=henry_air*x2w
-               for x0g in case wPhaseOnly we compute x0g=pwsat
-               for x1g in case wPhaseOnly we compute x1g=henry_NAPL*x1w
+               for xgg in case wPhaseOnly we compute xgg=henry_air*x2w
+               for xwg in case wPhaseOnly we compute xwg=pwsat
+               for xng in case wPhaseOnly we compute xng=henry_NAPL*x1w
             */
 
             Scalar xgMax = 1.0;
-            if (x0g + x2g + x1g > xgMax)
+            if (xwg + xgg + xng > xgMax)
                 wouldSwitch = true;
             if (staticVertexDat_[globalIdx].wasSwitched)
                 xgMax *= 1.02;
 
             // if the sum of the mole fractions would be larger than
             // 100%, gas phase appears
-            if (x0g + x2g + x1g > xgMax)
+            if (xwg + xgg + xng > xgMax)
             {
                 // gas phase appears
                 std::cout << "gas phase appears at vertex " << globalIdx
-                          << ", coordinates: " << globalPos << ", x0g + x2g + x1g: "
-                          << x0g + x2g + x1g << std::endl;
+                          << ", coordinates: " << globalPos << ", xwg + xgg + xng: "
+                          << xwg + xgg + xng << std::endl;
                 gasFlag = 1;
             }
 
             // calculate fractions in the hypothetical NAPL phase
-            Scalar x1n = volVars.fluidState().moleFraction(nPhaseIdx, comp1Idx);
+            Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
             /* take care:
-               for x1n in case wPhaseOnly we compute x1n=henry_mesitylene*x1w,
+               for xnn in case wPhaseOnly we compute xnn=henry_mesitylene*x1w,
                where a hypothetical gas pressure is assumed for the Henry
                x0n is set to NULL  (all NAPL phase is dirty)
                x2n is set to NULL  (all NAPL phase is dirty)
             */
 
             Scalar xnMax = 1.0;
-            if (x1n > xnMax)
+            if (xnn > xnMax)
                 wouldSwitch = true;
             if (staticVertexDat_[globalIdx].wasSwitched)
                 xnMax *= 1.02;
 
             // if the sum of the hypothetical mole fractions would be larger than
             // 100%, NAPL phase appears
-            if (x1n > xnMax)
+            if (xnn > xnMax)
             {
                 // NAPL phase appears
                 std::cout << "NAPL phase appears at vertex " << globalIdx
-                          << ", coordinates: " << globalPos << ", x1n: "
-                          << x1n << std::endl;
-                NAPLFlag = 1;
+                          << ", coordinates: " << globalPos << ", xnn: "
+                          << xnn << std::endl;
+                nonwettingFlag = 1;
             }
 
-            if ((gasFlag == 1) && (NAPLFlag == 0))
+            if ((gasFlag == 1) && (nonwettingFlag == 0))
             {
                 newPhasePresence = wgPhaseOnly;
                 globalSol[globalIdx][switch1Idx] = 0.9999;
                 globalSol[globalIdx][switch2Idx] = 0.0001;
             }
-            else if ((gasFlag == 1) && (NAPLFlag == 1))
+            else if ((gasFlag == 1) && (nonwettingFlag == 1))
             {
                 newPhasePresence = threePhases;
                 globalSol[globalIdx][switch1Idx] = 0.9999;
                 globalSol[globalIdx][switch2Idx] = 0.0001;
             }
-            else if ((gasFlag == 0) && (NAPLFlag == 1))
+            else if ((gasFlag == 0) && (nonwettingFlag == 1))
             {
                 newPhasePresence = wnPhaseOnly;
                 globalSol[globalIdx][switch1Idx]
-                    = volVars.fluidState().moleFraction(wPhaseIdx, comp2Idx);
+                    = volVars.fluidState().moleFraction(wPhaseIdx, gCompIdx);
                 globalSol[globalIdx][switch2Idx] = 0.0001;
             }
         }
         else if (phasePresence == gnPhaseOnly)
         {
-            bool NAPLFlag = 0;
-            bool waterFlag = 0;
+            bool nonwettingFlag = 0;
+            bool wettingFlag = 0;
 
             Scalar Smin = 0.0;
             if (staticVertexDat_[globalIdx].wasSwitched)
@@ -680,58 +680,58 @@ protected:
                 std::cout << "NAPL phase disappears at vertex " << globalIdx
                           << ", coordinates: " << globalPos << ", Sn: "
                           << volVars.saturation(nPhaseIdx) << std::endl;
-                NAPLFlag = 1;
+                nonwettingFlag = 1;
             }
 
 
             // calculate fractions of the hypothetical water phase
-            Scalar x0w = volVars.fluidState().moleFraction(wPhaseIdx, comp0Idx);
+            Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
             /*
-              take care:, x0w, if no water is present, then take x0w=x0g*pg/pwsat .
+              take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
               If this is larger than 1, then water appears
             */
             Scalar xwMax = 1.0;
-            if (x0w > xwMax)
+            if (xww > xwMax)
                 wouldSwitch = true;
             if (staticVertexDat_[globalIdx].wasSwitched)
                 xwMax *= 1.02;
 
             // if the sum of the mole fractions would be larger than
             // 100%, gas phase appears
-            if (x0w > xwMax)
+            if (xww > xwMax)
             {
                 // water phase appears
                 std::cout << "water phase appears at vertex " << globalIdx
-                          << ", coordinates: " << globalPos << ", x0w=x0g*pg/pwsat : "
-                          << x0w << std::endl;
-                waterFlag = 1;
+                          << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
+                          << xww << std::endl;
+                wettingFlag = 1;
             }
 
-            if ((waterFlag == 1) && (NAPLFlag == 0))
+            if ((wettingFlag == 1) && (nonwettingFlag == 0))
             {
                 newPhasePresence = threePhases;
                 globalSol[globalIdx][switch1Idx] = 0.0001;
                 globalSol[globalIdx][switch2Idx] = volVars.saturation(nPhaseIdx);
             }
-            else if ((waterFlag == 1) && (NAPLFlag == 1))
+            else if ((wettingFlag == 1) && (nonwettingFlag == 1))
             {
                 newPhasePresence = wgPhaseOnly;
                 globalSol[globalIdx][switch1Idx] = 0.0001;
                 globalSol[globalIdx][switch2Idx]
-                    = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
+                    = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
             }
-            else if ((waterFlag == 0) && (NAPLFlag == 1))
+            else if ((wettingFlag == 0) && (nonwettingFlag == 1))
             {
                 newPhasePresence = gPhaseOnly;
                 globalSol[globalIdx][switch1Idx]
-                    = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
+                    = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
                 globalSol[globalIdx][switch2Idx]
-                    = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
+                    = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
             }
         }
         else if (phasePresence == wnPhaseOnly)
         {
-            bool NAPLFlag = 0;
+            bool nonwettingFlag = 0;
             bool gasFlag = 0;
 
             Scalar Smin = 0.0;
@@ -745,156 +745,156 @@ protected:
                 std::cout << "NAPL phase disappears at vertex " << globalIdx
                           << ", coordinates: " << globalPos << ", Sn: "
                           << volVars.saturation(nPhaseIdx) << std::endl;
-                NAPLFlag = 1;
+                nonwettingFlag = 1;
             }
 
             // calculate fractions of the partial pressures in the
             // hypothetical gas phase
-            Scalar x0g = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
-            Scalar x2g = volVars.fluidState().moleFraction(gPhaseIdx, comp2Idx);
-            Scalar x1g = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
+            Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+            Scalar xgg = volVars.fluidState().moleFraction(gPhaseIdx, gCompIdx);
+            Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
             /* take care:
-               for x2g in case wPhaseOnly we compute x2g=henry_air*x2w
-               for x0g in case wPhaseOnly we compute x0g=pwsat
-               for x1g in case wPhaseOnly we compute x1g=henry_NAPL*x1w
+               for xgg in case wPhaseOnly we compute xgg=henry_air*x2w
+               for xwg in case wPhaseOnly we compute xwg=pwsat
+               for xng in case wPhaseOnly we compute xng=henry_NAPL*x1w
             */
             Scalar xgMax = 1.0;
-            if (x0g + x2g + x1g > xgMax)
+            if (xwg + xgg + xng > xgMax)
                 wouldSwitch = true;
             if (staticVertexDat_[globalIdx].wasSwitched)
                 xgMax *= 1.02;
 
             // if the sum of the mole fractions would be larger than
             // 100%, gas phase appears
-            if (x0g + x2g + x1g > xgMax)
+            if (xwg + xgg + xng > xgMax)
             {
                 // gas phase appears
                 std::cout << "gas phase appears at vertex " << globalIdx
-                          << ", coordinates: " << globalPos << ", x0g + x2g + x1g: "
-                          << x0g + x2g + x1g << std::endl;
+                          << ", coordinates: " << globalPos << ", xwg + xgg + xng: "
+                          << xwg + xgg + xng << std::endl;
                 gasFlag = 1;
             }
 
-            if ((gasFlag == 1) && (NAPLFlag == 0))
+            if ((gasFlag == 1) && (nonwettingFlag == 0))
             {
                 newPhasePresence = threePhases;
                 globalSol[globalIdx][switch1Idx] = volVars.saturation(wPhaseIdx);
                 globalSol[globalIdx][switch2Idx] = volVars.saturation(nPhaseIdx);;
             }
-            else if ((gasFlag == 1) && (NAPLFlag == 1))
+            else if ((gasFlag == 1) && (nonwettingFlag == 1))
             {
                 newPhasePresence = wgPhaseOnly;
                 globalSol[globalIdx][switch1Idx] = volVars.saturation(wPhaseIdx);
                 globalSol[globalIdx][switch2Idx]
-                    = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
+                    = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
             }
-            else if ((gasFlag == 0) && (NAPLFlag == 1))
+            else if ((gasFlag == 0) && (nonwettingFlag == 1))
             {
                 newPhasePresence = wPhaseOnly;
                 globalSol[globalIdx][switch1Idx]
-                    = volVars.fluidState().moleFraction(wPhaseIdx, comp2Idx);
+                    = volVars.fluidState().moleFraction(wPhaseIdx, gCompIdx);
                 globalSol[globalIdx][switch2Idx]
-                    = volVars.fluidState().moleFraction(wPhaseIdx, comp1Idx);
+                    = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
             }
         }
         else if (phasePresence == gPhaseOnly)
         {
-            bool NAPLFlag = 0;
-            bool waterFlag = 0;
+            bool nonwettingFlag = 0;
+            bool wettingFlag = 0;
 
             // calculate fractions in the hypothetical NAPL phase
-            Scalar x1n = volVars.fluidState().moleFraction(nPhaseIdx, comp1Idx);
+            Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
             /*
-              take care:, x1n, if no NAPL phase is there, take x1n=x1g*pg/pcsat
+              take care:, xnn, if no NAPL phase is there, take xnn=xng*pg/pcsat
               if this is larger than 1, then NAPL appears
             */
 
             Scalar xnMax = 1.0;
-            if (x1n > xnMax)
+            if (xnn > xnMax)
                 wouldSwitch = true;
             if (staticVertexDat_[globalIdx].wasSwitched)
                 xnMax *= 1.02;
 
             // if the sum of the hypothetical mole fraction would be larger than
             // 100%, NAPL phase appears
-            if (x1n > xnMax)
+            if (xnn > xnMax)
             {
                 // NAPL phase appears
                 std::cout << "NAPL phase appears at vertex " << globalIdx
-                          << ", coordinates: " << globalPos << ", x1n: "
-                          << x1n << std::endl;
-                NAPLFlag = 1;
+                          << ", coordinates: " << globalPos << ", xnn: "
+                          << xnn << std::endl;
+                nonwettingFlag = 1;
             }
             // calculate fractions of the hypothetical water phase
-            Scalar x0w = volVars.fluidState().moleFraction(wPhaseIdx, comp0Idx);
+            Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
             /*
-              take care:, x0w, if no water is present, then take x0w=x0g*pg/pwsat .
+              take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
               If this is larger than 1, then water appears
             */
             Scalar xwMax = 1.0;
-            if (x0w > xwMax)
+            if (xww > xwMax)
                 wouldSwitch = true;
             if (staticVertexDat_[globalIdx].wasSwitched)
                 xwMax *= 1.02;
 
             // if the sum of the mole fractions would be larger than
             // 100%, gas phase appears
-            if (x0w > xwMax)
+            if (xww > xwMax)
             {
                 // water phase appears
                 std::cout << "water phase appears at vertex " << globalIdx
-                          << ", coordinates: " << globalPos << ", x0w=x0g*pg/pwsat : "
-                          << x0w << std::endl;
-                waterFlag = 1;
+                          << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
+                          << xww << std::endl;
+                wettingFlag = 1;
             }
-            if ((waterFlag == 1) && (NAPLFlag == 0))
+            if ((wettingFlag == 1) && (nonwettingFlag == 0))
             {
                 newPhasePresence = wgPhaseOnly;
                 globalSol[globalIdx][switch1Idx] = 0.0001;
                 globalSol[globalIdx][switch2Idx]
-                    = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
+                    = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
             }
-            else if ((waterFlag == 1) && (NAPLFlag == 1))
+            else if ((wettingFlag == 1) && (nonwettingFlag == 1))
             {
                 newPhasePresence = threePhases;
                 globalSol[globalIdx][switch1Idx] = 0.0001;
                 globalSol[globalIdx][switch2Idx] = 0.0001;
             }
-            else if ((waterFlag == 0) && (NAPLFlag == 1))
+            else if ((wettingFlag == 0) && (nonwettingFlag == 1))
             {
                 newPhasePresence = gnPhaseOnly;
                 globalSol[globalIdx][switch1Idx]
-                    = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
+                    = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
                 globalSol[globalIdx][switch2Idx] = 0.0001;
             }
         }
         else if (phasePresence == wgPhaseOnly)
         {
-            bool NAPLFlag = 0;
+            bool nonwettingFlag = 0;
             bool gasFlag = 0;
-            bool waterFlag = 0;
+            bool wettingFlag = 0;
 
             // get the fractions in the hypothetical NAPL phase
-            Scalar x1n = volVars.fluidState().moleFraction(nPhaseIdx, comp1Idx);
+            Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
 
             // take care: if the NAPL phase is not present, take
-            // x1n=x1g*pg/pcsat if this is larger than 1, then NAPL
+            // xnn=xng*pg/pcsat if this is larger than 1, then NAPL
             // appears
             Scalar xnMax = 1.0;
-            if (x1n > xnMax)
+            if (xnn > xnMax)
                 wouldSwitch = true;
             if (staticVertexDat_[globalIdx].wasSwitched)
                 xnMax *= 1.02;
 
             // if the sum of the hypothetical mole fraction would be larger than
             // 100%, NAPL phase appears
-            if (x1n > xnMax)
+            if (xnn > xnMax)
             {
                 // NAPL phase appears
                 std::cout << "NAPL phase appears at vertex " << globalIdx
-                          << ", coordinates: " << globalPos << ", x1n: "
-                          << x1n << std::endl;
-                NAPLFlag = 1;
+                          << ", coordinates: " << globalPos << ", xnn: "
+                          << xnn << std::endl;
+                nonwettingFlag = 1;
             }
 
             Scalar Smin = -1.e-6;
@@ -922,37 +922,37 @@ protected:
                 std::cout << "Water phase disappears at vertex " << globalIdx
                           << ", coordinates: " << globalPos << ", Sw: "
                           << volVars.saturation(wPhaseIdx) << std::endl;
-                waterFlag = 1;
+                wettingFlag = 1;
             }
 
-            if ((gasFlag == 0) && (NAPLFlag == 1) && (waterFlag == 1))
+            if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 1))
             {
                 newPhasePresence = gnPhaseOnly;
                 globalSol[globalIdx][switch1Idx]
-                    = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
+                    = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
                 globalSol[globalIdx][switch2Idx] = 0.0001;
             }
-            else if ((gasFlag == 0) && (NAPLFlag == 1) && (waterFlag == 0))
+            else if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 0))
             {
                 newPhasePresence = threePhases;
                 globalSol[globalIdx][switch1Idx] = volVars.saturation(wPhaseIdx);
                 globalSol[globalIdx][switch2Idx] = 0.0;
             }
-            else if ((gasFlag == 1) && (NAPLFlag == 0) && (waterFlag == 0))
+            else if ((gasFlag == 1) && (nonwettingFlag == 0) && (wettingFlag == 0))
             {
                 newPhasePresence = wPhaseOnly;
                 globalSol[globalIdx][switch1Idx]
-                    = volVars.fluidState().moleFraction(wPhaseIdx, comp2Idx);
+                    = volVars.fluidState().moleFraction(wPhaseIdx, gCompIdx);
                 globalSol[globalIdx][switch2Idx]
-                    = volVars.fluidState().moleFraction(wPhaseIdx, comp1Idx);
+                    = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
             }
-            else if ((gasFlag == 0) && (NAPLFlag == 0) && (waterFlag == 1))
+            else if ((gasFlag == 0) && (nonwettingFlag == 0) && (wettingFlag == 1))
             {
                 newPhasePresence = gPhaseOnly;
                 globalSol[globalIdx][switch1Idx]
-                    = volVars.fluidState().moleFraction(gPhaseIdx, comp0Idx);
+                    = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
                 globalSol[globalIdx][switch2Idx]
-                    = volVars.fluidState().moleFraction(gPhaseIdx, comp1Idx);
+                    = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
             }
         }
 
diff --git a/dumux/boxmodels/3p3c/3p3cvolumevariables.hh b/dumux/boxmodels/3p3c/3p3cvolumevariables.hh
index c43f398654bb0de74f000143872f6fd2811de95b..8498cb1b4f717fa68f39c2b877ae2f2ce296a2e4 100644
--- a/dumux/boxmodels/3p3c/3p3cvolumevariables.hh
+++ b/dumux/boxmodels/3p3c/3p3cvolumevariables.hh
@@ -79,9 +79,9 @@ class ThreePThreeCVolumeVariables : public BoxVolumeVariables<TypeTag>
         numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
         numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
 
-        comp0Idx = Indices::comp0Idx,
-        comp2Idx = Indices::comp2Idx,
-        comp1Idx = Indices::comp1Idx,
+        wCompIdx = Indices::wCompIdx,
+        gCompIdx = Indices::gCompIdx,
+        nCompIdx = Indices::nCompIdx,
 
         wPhaseIdx = Indices::wPhaseIdx,
         gPhaseIdx = Indices::gPhaseIdx,
@@ -230,14 +230,14 @@ public:
             // stored explicitly.
 
             // extract mole fractions in the water phase
-            Scalar xw2 = priVars[switch1Idx];
-            Scalar xw1 = priVars[switch2Idx];
-            Scalar xw0 = 1 - xw2 - xw1;
+            Scalar xwg = priVars[switch1Idx];
+            Scalar xwn = priVars[switch2Idx];
+            Scalar xww = 1 - xwg - xwn;
 
             // write water mole fractions in the fluid state
-            fluidState_.setMoleFraction(wPhaseIdx, comp0Idx, xw0);
-            fluidState_.setMoleFraction(wPhaseIdx, comp2Idx, xw2);
-            fluidState_.setMoleFraction(wPhaseIdx, comp1Idx, xw1);
+            fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
+            fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
+            fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
 
             // calculate the composition of the remaining phases (as
             // well as the densities of all phases). this is the job
@@ -257,16 +257,16 @@ public:
             // we have all (partly hypothetical) phase pressures
             // and temperature and the mole fraction of water in
             // the gas phase
-            Scalar partPressNAPL = fluidState_.fugacityCoefficient(nPhaseIdx, comp1Idx)*pn_;
+            Scalar partPressNAPL = fluidState_.fugacityCoefficient(nPhaseIdx, nCompIdx)*pn_;
 
-            Scalar xg0 = priVars[switch1Idx];
-            Scalar xg1 = partPressNAPL/pg_;
-            Scalar xg2 = 1.-xg0-xg1;
+            Scalar xgw = priVars[switch1Idx];
+            Scalar xgn = partPressNAPL/pg_;
+            Scalar xgg = 1.-xgw-xgn;
 
             // write mole fractions in the fluid state
-            fluidState_.setMoleFraction(gPhaseIdx, comp0Idx, xg0);
-            fluidState_.setMoleFraction(gPhaseIdx, comp2Idx, xg2);
-            fluidState_.setMoleFraction(gPhaseIdx, comp1Idx, xg1);
+            fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
+            fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
+            fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
 
             // calculate the composition of the remaining phases (as
             // well as the densities of all phases). this is the job
@@ -279,17 +279,17 @@ public:
         }
         else if (phasePresence == wnPhaseOnly) {
             // only water and NAPL phases are present
-            Scalar pPartialC = fluidState_.fugacityCoefficient(nPhaseIdx,comp1Idx)*pn_;
-            Scalar henryC = fluidState_.fugacityCoefficient(wPhaseIdx,comp1Idx)*pw_;
+            Scalar pPartialC = fluidState_.fugacityCoefficient(nPhaseIdx,nCompIdx)*pn_;
+            Scalar henryC = fluidState_.fugacityCoefficient(wPhaseIdx,nCompIdx)*pw_;
 
-            Scalar xw2 = priVars[switch1Idx];
-            Scalar xw1 = pPartialC/henryC;
-            Scalar xw0 = 1.-xw2-xw1;
+            Scalar xwg = priVars[switch1Idx];
+            Scalar xwn = pPartialC/henryC;
+            Scalar xww = 1.-xwg-xwn;
 
             // write mole fractions in the fluid state
-            fluidState_.setMoleFraction(wPhaseIdx, comp0Idx, xw0);
-            fluidState_.setMoleFraction(wPhaseIdx, comp2Idx, xw2);
-            fluidState_.setMoleFraction(wPhaseIdx, comp1Idx, xw1);
+            fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
+            fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
+            fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
 
             // calculate the composition of the remaining phases (as
             // well as the densities of all phases). this is the job
@@ -304,14 +304,14 @@ public:
             // only the gas phase is present, gas phase composition is
             // stored explicitly here below.
 
-            const Scalar xg0 = priVars[switch1Idx];
-            const Scalar xg1 = priVars[switch2Idx];
-            Scalar xg2 = 1 - xg0 - xg1;
+            const Scalar xgw = priVars[switch1Idx];
+            const Scalar xgn = priVars[switch2Idx];
+            Scalar xgg = 1 - xgw - xgn;
 
             // write mole fractions in the fluid state
-            fluidState_.setMoleFraction(gPhaseIdx, comp0Idx, xg0);
-            fluidState_.setMoleFraction(gPhaseIdx, comp2Idx, xg2);
-            fluidState_.setMoleFraction(gPhaseIdx, comp1Idx, xg1);
+            fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
+            fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
+            fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
 
             // calculate the composition of the remaining phases (as
             // well as the densities of all phases). this is the job
@@ -324,16 +324,16 @@ public:
         }
         else if (phasePresence == wgPhaseOnly) {
             // only water and gas phases are present
-            Scalar xg1 = priVars[switch2Idx];
-            Scalar partPressH2O = fluidState_.fugacityCoefficient(wPhaseIdx, comp0Idx)*pw_;
+            Scalar xgn = priVars[switch2Idx];
+            Scalar partPressH2O = fluidState_.fugacityCoefficient(wPhaseIdx, wCompIdx)*pw_;
 
-            Scalar xg0 = partPressH2O/pg_;
-            Scalar xg2 = 1.-xg1-xg0;
+            Scalar xgw = partPressH2O/pg_;
+            Scalar xgg = 1.-xgn-xgw;
 
             // write mole fractions in the fluid state
-            fluidState_.setMoleFraction(gPhaseIdx, comp0Idx, xg0);
-            fluidState_.setMoleFraction(gPhaseIdx, comp2Idx, xg2);
-            fluidState_.setMoleFraction(gPhaseIdx, comp1Idx, xg1);
+            fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
+            fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
+            fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
 
             // calculate the composition of the remaining phases (as
             // well as the densities of all phases). this is the job
@@ -373,34 +373,34 @@ public:
          */
 
         // diffusivity coefficents
-        diffusionCoefficient_[gPhaseIdx][comp0Idx] =
+        diffusionCoefficient_[gPhaseIdx][wCompIdx] =
             FluidSystem::diffusionCoefficient(fluidState_,
                                               paramCache,
                                               gPhaseIdx,
-                                              comp0Idx);
-        diffusionCoefficient_[gPhaseIdx][comp1Idx] =
+                                              wCompIdx);
+        diffusionCoefficient_[gPhaseIdx][nCompIdx] =
             FluidSystem::diffusionCoefficient(fluidState_,
                                               paramCache,
                                               gPhaseIdx,
-                                              comp1Idx);
-        diffusionCoefficient_[gPhaseIdx][comp2Idx] = 0.0; // dummy, should not be used !
+                                              nCompIdx);
+        diffusionCoefficient_[gPhaseIdx][gCompIdx] = 0.0; // dummy, should not be used !
 
-        diffusionCoefficient_[wPhaseIdx][comp2Idx] =
+        diffusionCoefficient_[wPhaseIdx][gCompIdx] =
             FluidSystem::diffusionCoefficient(fluidState_,
                                               paramCache,
                                               wPhaseIdx,
-                                              comp2Idx);
-        diffusionCoefficient_[wPhaseIdx][comp1Idx] =
+                                              gCompIdx);
+        diffusionCoefficient_[wPhaseIdx][nCompIdx] =
             FluidSystem::diffusionCoefficient(fluidState_,
                                               paramCache,
                                               wPhaseIdx,
-                                              comp1Idx);
-        diffusionCoefficient_[wPhaseIdx][comp0Idx] = 0.0; // dummy, should not be used !
+                                              nCompIdx);
+        diffusionCoefficient_[wPhaseIdx][wCompIdx] = 0.0; // dummy, should not be used !
 
         /* no diffusion in NAPL phase considered  at the moment */
-        diffusionCoefficient_[nPhaseIdx][comp1Idx] = 0.0;
-        diffusionCoefficient_[nPhaseIdx][comp0Idx] = 0.0;
-        diffusionCoefficient_[nPhaseIdx][comp2Idx] = 0.0;
+        diffusionCoefficient_[nPhaseIdx][nCompIdx] = 0.0;
+        diffusionCoefficient_[nPhaseIdx][wCompIdx] = 0.0;
+        diffusionCoefficient_[nPhaseIdx][gCompIdx] = 0.0;
 
         Valgrind::CheckDefined(diffusionCoefficient_);
 
diff --git a/dumux/boxmodels/richards/richardsindices.hh b/dumux/boxmodels/richards/richardsindices.hh
index 0980c46071fd90e91666a2b53232be1543e4406c..43a873075cd4f35eb53a28c4983d0e945a6b7b32 100644
--- a/dumux/boxmodels/richards/richardsindices.hh
+++ b/dumux/boxmodels/richards/richardsindices.hh
@@ -36,28 +36,30 @@ namespace Dumux
  * \ingroup BoxIndices
  * \brief Index names for the Richards model.
  */
+
+template <class TypeTag>
 struct RichardsIndices
 {
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    
     //////////
     // primary variable indices
     //////////
-
+    
     //! Primary variable index for the wetting phase pressure
     static const int pwIdx = 0;
-
+    
     //////////
     // equation indices
     //////////
     //! Equation index for the mass conservation of the wetting phase
     static const int contiEqIdx = 0;
-
+    
     //////////
     // phase indices
     //////////
-    //! Phase index for the wetting phase
-    static const int wPhaseIdx = 0;
-    //! Phase index for the wetting phase
-    static const int nPhaseIdx = 1;
+    static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< Index of the wetting phase;
+    static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< Index of the non-wetting phase;
 };
 // \}
 
diff --git a/dumux/boxmodels/richards/richardspropertydefaults.hh b/dumux/boxmodels/richards/richardspropertydefaults.hh
index cd757ed928fcae780d7fe30fbed5426d97da9082..599a14b5183b024c36ba242e88bdd235441ddb46 100644
--- a/dumux/boxmodels/richards/richardspropertydefaults.hh
+++ b/dumux/boxmodels/richards/richardspropertydefaults.hh
@@ -78,7 +78,7 @@ SET_SCALAR_PROP(BoxRichards,
                 1.0);
 
 //! The class with all index definitions for the model
-SET_TYPE_PROP(BoxRichards, Indices, Dumux::RichardsIndices);
+SET_TYPE_PROP(BoxRichards, Indices, RichardsIndices<TypeTag>);
 
 //! DEPRECATED RichardsIndices property
 SET_TYPE_PROP(BoxRichards, RichardsIndices, typename GET_PROP_TYPE(TypeTag, Indices));
diff --git a/dumux/linear/seqsolverbackend.hh b/dumux/linear/seqsolverbackend.hh
index eb6598c6b614d47a1c01e6081fd13e5deb4d021e..276a78f46b3093e787b14d28520415e2b924d2d7 100644
--- a/dumux/linear/seqsolverbackend.hh
+++ b/dumux/linear/seqsolverbackend.hh
@@ -420,7 +420,6 @@ public:
     Vector bTmp(b);
 
     const double relaxation = GET_PARAM(TypeTag, double, PreconditionerRelaxation);
-    const int precondIter = GET_PARAM(TypeTag, int, PreconditionerIterations);
 
     Preconditioner precond(A, relaxation);
 
diff --git a/dumux/material/fluidsystems/h2oairfluidsystem.hh b/dumux/material/fluidsystems/h2oairfluidsystem.hh
index f66949c0ead16f3c62da6fb62b59f87396cdf895..b926433b7057a579f92d70aaef54e9b4ebed6e3d 100644
--- a/dumux/material/fluidsystems/h2oairfluidsystem.hh
+++ b/dumux/material/fluidsystems/h2oairfluidsystem.hh
@@ -101,11 +101,11 @@ public:
 
     static constexpr int numPhases = 2;
 
-    static constexpr int lPhaseIdx = 0; // index of the liquid phase
-    static constexpr int gPhaseIdx = 1; // index of the gas phase
+    static constexpr int wPhaseIdx = 0; // index of the wetting phase
+    static constexpr int nPhaseIdx = 1; // index of the non-wetting phase
 
-    static constexpr int wPhaseIdx = lPhaseIdx; // index of the wetting phase
-    static constexpr int nPhaseIdx = gPhaseIdx; // index of the non-wetting phase
+    static constexpr int lPhaseIdx = wPhaseIdx; // DEPRECATED index of the liquid phase
+    static constexpr int gPhaseIdx = nPhaseIdx; // DEPRECATED index of the gas phase
 
     /*!
      * \brief Return the human readable name of a phase
diff --git a/dumux/material/fluidsystems/h2oairmesitylenefluidsystem.hh b/dumux/material/fluidsystems/h2oairmesitylenefluidsystem.hh
index 583d5de5be8061ad94358f80f4cc5cf6c0b207eb..9b29a41da7679380afb43366f2765ebab37971d2 100644
--- a/dumux/material/fluidsystems/h2oairmesitylenefluidsystem.hh
+++ b/dumux/material/fluidsystems/h2oairmesitylenefluidsystem.hh
@@ -80,6 +80,13 @@ public:
     static const int H2OIdx = 0;
     static const int NAPLIdx = 1;
     static const int airIdx = 2;
+    
+    // export component indices to indicate the main component 
+    // of the corresponding phase at atmospheric pressure 1 bar 
+    // and room temperature 20°C:
+    static const int wCompIdx = H2OIdx;
+    static const int nCompIdx = NAPLIdx;
+    static const int gCompIdx = airIdx;
 
     /*!
      * \brief Initialize the fluid system's static parameters generically
diff --git a/dumux/material/fluidsystems/h2oairxylenefluidsystem.hh b/dumux/material/fluidsystems/h2oairxylenefluidsystem.hh
index 22eac194bcb088cf2306ba0d82dfdc80e95207af..4274a183232edc6dac59ab7642e52515971938b2 100644
--- a/dumux/material/fluidsystems/h2oairxylenefluidsystem.hh
+++ b/dumux/material/fluidsystems/h2oairxylenefluidsystem.hh
@@ -73,6 +73,13 @@ public:
     static const int NAPLIdx = 1;
     static const int airIdx = 2;
 
+    // export component indices to indicate the main component 
+    // of the corresponding phase at atmospheric pressure 1 bar 
+    // and room temperature 20°C:
+    static const int wCompIdx = H2OIdx;
+    static const int nCompIdx = NAPLIdx;
+    static const int gCompIdx = airIdx;
+
     static void init()
     { }
 
diff --git a/dumux/material/fluidsystems/h2on2fluidsystem.hh b/dumux/material/fluidsystems/h2on2fluidsystem.hh
index f0edd2a61aef6320448bc7255b55ffc3b0fdeab0..725e1b4a53876ec84c97c01c291d5bd2796c7d54 100644
--- a/dumux/material/fluidsystems/h2on2fluidsystem.hh
+++ b/dumux/material/fluidsystems/h2on2fluidsystem.hh
@@ -84,13 +84,18 @@ public:
     //! Number of phases in the fluid system
     static constexpr int numPhases = 2;
 
-    //! Index of the liquid phase
-    static constexpr int lPhaseIdx = 0;
-    static constexpr int wPhaseIdx = lPhaseIdx;
-    //! Index of the gas phase
-    static constexpr int gPhaseIdx = 1;
-    static constexpr int nPhaseIdx = gPhaseIdx;
-
+    static constexpr int wPhaseIdx = 0; // index of the wetting phase
+    static constexpr int nPhaseIdx = 1; // index of the non-wetting phase
+    
+    static constexpr int lPhaseIdx = wPhaseIdx; // DEPRECATED index of the liquid phase
+    static constexpr int gPhaseIdx = nPhaseIdx; // DEPRECATED index of the gas phase
+
+    // export component indices to indicate the main component 
+    // of the corresponding phase at atmospheric pressure 1 bar 
+    // and room temperature 20°C:
+    static const int wCompIdx = wPhaseIdx;
+    static const int nCompIdx = nPhaseIdx;
+    
     /*!
      * \brief Return the human readable name of a fluid phase
      *
diff --git a/dumux/material/fluidsystems/h2on2liquidphasefluidsystem.hh b/dumux/material/fluidsystems/h2on2liquidphasefluidsystem.hh
index 7552d68c3fb12c7a1836e6e08dd4f89b8f86cb8a..155f7a6a06da0748dc23f12e2aa8b70df5de7b2c 100644
--- a/dumux/material/fluidsystems/h2on2liquidphasefluidsystem.hh
+++ b/dumux/material/fluidsystems/h2on2liquidphasefluidsystem.hh
@@ -85,9 +85,8 @@ public:
     //! Number of phases in the fluid system
     static constexpr int numPhases = 1;
 
-    //! Index of the liquid phase
-    static constexpr int lPhaseIdx = 0;
-    static constexpr int wPhaseIdx = lPhaseIdx;
+    static constexpr int wPhaseIdx = 0; // index of the wetting phase
+    static constexpr int lPhaseIdx = wPhaseIdx; // DEPRECATED index of the liquid phase
 
     /*!
      * \brief Return the human readable name of a fluid phase.