diff --git a/dumux/discretization/staggered/freeflow/fickslaw.hh b/dumux/discretization/staggered/freeflow/fickslaw.hh
index 1cd6fa434fdeb1c1e2b5d703c9f957acb29557e4..102ce6d88afbe84b9fe472bd4758bdf41b5c7828 100644
--- a/dumux/discretization/staggered/freeflow/fickslaw.hh
+++ b/dumux/discretization/staggered/freeflow/fickslaw.hh
@@ -70,8 +70,7 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::staggered >
     static_assert(ModelTraits::numPhases() == 1, "Only one phase allowed supported!");
 
     enum {
-        conti0EqIdx = Indices::conti0EqIdx,
-        mainCompIdx = Indices::mainCompIdx,
+        conti0EqIdx = Indices::conti0EqIdx
     };
 
 public:
@@ -97,7 +96,7 @@ public:
 
         for(int compIdx = 0; compIdx < numComponents; ++compIdx)
         {
-            if(compIdx == mainCompIdx)
+            if(compIdx == FluidSystem::getMainComponent(0))
                 continue;
 
             // get equation index
@@ -120,7 +119,7 @@ public:
         }
 
         const Scalar cumulativeFlux = std::accumulate(flux.begin(), flux.end(), 0.0);
-        flux[mainCompIdx] = - cumulativeFlux;
+        flux[FluidSystem::getMainComponent(0)] = - cumulativeFlux;
 
         // Fick's law (for binary systems) states that the net flux of moles within the bulk phase has to be zero:
         // If a given amount of molecules A travel into one direction, the same amount of molecules B have to
@@ -146,7 +145,7 @@ public:
         const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
 
         const Scalar insideDistance = (insideScv.dofPosition() - scvf.ipGlobal()).two_norm();
-        const Scalar insideD = insideVolVars.effectiveDiffusivity(compIdx);
+        const Scalar insideD = insideVolVars.effectiveDiffusivity(0, compIdx);
         const Scalar ti = calculateOmega_(insideDistance, insideD, insideVolVars.extrusionFactor());
 
         if(scvf.boundary())
@@ -156,7 +155,7 @@ public:
             const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
             const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
             const Scalar outsideDistance = (outsideScv.dofPosition() - scvf.ipGlobal()).two_norm();
-            const Scalar outsideD = outsideVolVars.effectiveDiffusivity(compIdx);
+            const Scalar outsideD = outsideVolVars.effectiveDiffusivity(0, compIdx);
             const Scalar tj = calculateOmega_(outsideDistance, outsideD, outsideVolVars.extrusionFactor());
 
             tij = scvf.area()*(ti * tj)/(ti + tj);
diff --git a/dumux/freeflow/compositional/volumevariables.hh b/dumux/freeflow/compositional/volumevariables.hh
index 286d6f456e42784f45515fc6f798e3e1128b42bd..90d576626a2babacaf83bdc9cf3af359c16f1a7c 100644
--- a/dumux/freeflow/compositional/volumevariables.hh
+++ b/dumux/freeflow/compositional/volumevariables.hh
@@ -43,7 +43,6 @@ class FreeflowNCVolumeVariables : public FreeFlowVolumeVariables< Traits, Freefl
 
     using Scalar = typename Traits::PrimaryVariables::value_type;
 
-    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
     static constexpr int numComponents = Traits::ModelTraits::numComponents();
 
     static constexpr bool useMoles = Traits::ModelTraits::useMoles();
@@ -77,6 +76,7 @@ public:
 
         typename FluidSystem::ParameterCache paramCache;
         paramCache.updateAll(fluidState_);
+
         for (unsigned int compIIdx = 0; compIIdx < numComponents; ++compIIdx)
         {
             for (unsigned int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
@@ -87,7 +87,7 @@ public:
                     diffCoefficient_[compIIdx][compJIdx]
                         = FluidSystem::binaryDiffusionCoefficient(fluidState_,
                                                                   paramCache,
-                                                                  fluidSystemPhaseIdx,
+                                                                  0,
                                                                   compIIdx,
                                                                   compJIdx);
                 }
@@ -106,65 +106,62 @@ public:
                                    FluidState& fluidState)
     {
         fluidState.setTemperature(ParentType::temperature(elemSol, problem, element, scv));
-        fluidState.setPressure(fluidSystemPhaseIdx, elemSol[0][Indices::pressureIdx]);
+        fluidState.setPressure(0, elemSol[0][Indices::pressureIdx]);
 
         // saturation in a single phase is always 1 and thus redundant
         // to set. But since we use the fluid state shared by the
         // immiscible multi-phase models, so we have to set it here...
-        fluidState.setSaturation(fluidSystemPhaseIdx, 1.0);
+        fluidState.setSaturation(0, 1.0);
 
         Scalar sumFracMinorComp = 0.0;
 
-        for(int compIdx = 0; compIdx < numComponents; ++compIdx)
+        for(int compIdx = 1; compIdx < numComponents; ++compIdx)
         {
-            if(compIdx == Indices::mainCompIdx)
-                continue;
-
             // temporary add 1.0 to remove spurious differences in mole fractions
             // which are below the numerical accuracy
             Scalar moleOrMassFraction = elemSol[0][Indices::conti0EqIdx+compIdx] + 1.0;
             moleOrMassFraction = moleOrMassFraction - 1.0;
             if(useMoles)
-                fluidState.setMoleFraction(fluidSystemPhaseIdx, compIdx, moleOrMassFraction);
+                fluidState.setMoleFraction(0, compIdx, moleOrMassFraction);
             else
-                fluidState.setMassFraction(fluidSystemPhaseIdx, compIdx, moleOrMassFraction);
+                fluidState.setMassFraction(0, compIdx, moleOrMassFraction);
             sumFracMinorComp += moleOrMassFraction;
         }
         if(useMoles)
-            fluidState.setMoleFraction(fluidSystemPhaseIdx, Indices::mainCompIdx, 1.0 - sumFracMinorComp);
+            fluidState.setMoleFraction(0, 0, 1.0 - sumFracMinorComp);
         else
-            fluidState.setMassFraction(fluidSystemPhaseIdx, Indices::mainCompIdx, 1.0 - sumFracMinorComp);
+            fluidState.setMassFraction(0, 0, 1.0 - sumFracMinorComp);
 
         typename FluidSystem::ParameterCache paramCache;
-        paramCache.updatePhase(fluidState, fluidSystemPhaseIdx);
+        paramCache.updateAll(fluidState);
 
-        Scalar value = FluidSystem::density(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setDensity(fluidSystemPhaseIdx, value);
+        Scalar value = FluidSystem::density(fluidState, paramCache, 0);
+        fluidState.setDensity(0, value);
 
-        value = FluidSystem::molarDensity(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setMolarDensity(fluidSystemPhaseIdx, value);
+        value = FluidSystem::molarDensity(fluidState, paramCache, 0);
+        fluidState.setMolarDensity(0, value);
 
-        value = FluidSystem::viscosity(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setViscosity(fluidSystemPhaseIdx, value);
+        value = FluidSystem::viscosity(fluidState, paramCache, 0);
+        fluidState.setViscosity(0, value);
 
         // compute and set the enthalpy
         const Scalar h = ParentType::enthalpy(fluidState, paramCache);
-        fluidState.setEnthalpy(fluidSystemPhaseIdx, h);
+        fluidState.setEnthalpy(0, h);
     }
 
     /*!
      * \brief Return the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
      *        the control volume.
      */
-    Scalar pressure(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.pressure(fluidSystemPhaseIdx); }
+    Scalar pressure(int phaseIdx = 0) const
+    { return fluidState_.pressure(0); }
 
     /*!
      * \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ of a given phase within the
      *        control volume.
      */
-    Scalar density(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.density(fluidSystemPhaseIdx); }
+    Scalar density(int phaseIdx = 0) const
+    { return fluidState_.density(0); }
 
     /*!
      * \brief Return temperature \f$\mathrm{[K]}\f$ inside the sub-control volume.
@@ -187,8 +184,8 @@ public:
      * \brief Return the dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
      *        control volume.
      */
-    Scalar viscosity(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.viscosity(fluidSystemPhaseIdx); }
+    Scalar viscosity(int phaseIdx = 0) const
+    { return fluidState_.viscosity(0); }
 
     /*!
      * \brief Returns the mass fraction of a component in the phase \f$\mathrm{[-]}\f$
@@ -198,7 +195,7 @@ public:
      */
     Scalar massFraction(int phaseIdx, int compIdx) const
     {
-        return fluidState_.massFraction(fluidSystemPhaseIdx, compIdx);
+        return fluidState_.massFraction(0, compIdx);
     }
 
     /*!
@@ -208,7 +205,7 @@ public:
      */
     Scalar massFraction(int compIdx) const
     {
-        return fluidState_.massFraction(fluidSystemPhaseIdx, compIdx);
+        return fluidState_.massFraction(0, compIdx);
     }
 
     /*!
@@ -219,7 +216,7 @@ public:
      */
     Scalar moleFraction(int phaseIdx, int compIdx) const
     {
-        return fluidState_.moleFraction(fluidSystemPhaseIdx, compIdx);
+        return fluidState_.moleFraction(0, compIdx);
     }
 
     /*!
@@ -229,15 +226,15 @@ public:
      */
     Scalar moleFraction(int compIdx) const
     {
-        return fluidState_.moleFraction(fluidSystemPhaseIdx, compIdx);
+        return fluidState_.moleFraction(0, compIdx);
     }
 
     /*!
      * \brief Returns the mass density of a given phase \f$\mathrm{[kg/m^3]}\f$
      */
-    Scalar molarDensity(int phaseIdx = fluidSystemPhaseIdx) const
+    Scalar molarDensity(int phaseIdx = 0) const
     {
-        return fluidState_.molarDensity(fluidSystemPhaseIdx);
+        return fluidState_.molarDensity(0);
     }
 
     /*!
@@ -250,29 +247,30 @@ public:
         return FluidSystem::molarMass(compIdx);
     }
 
-     /*!
-     * \brief Returns the diffusion coefficient \f$\mathrm{[m^2/s]}\f$
-     *
-     * \param compIIdx the index of the component which diffusive
-     * \param compJIdx the index of the component with respect to which compIIdx diffuses
-     */
-    Scalar diffusionCoefficient(int compIIdx, int compJIdx = fluidSystemPhaseIdx) const
-    {
-        if (compIIdx == compJIdx)
-            DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for fluidSystemPhaseIdx = compIdx");
-        return diffCoefficient_[compIIdx][compJIdx];
-    }
+    /*!
+    * \brief Returns the diffusion coefficient \f$\mathrm{[m^2/s]}\f$
+    *
+    * \param compIIdx the index of the component which diffusive
+    * \param compJIdx the index of the component with respect to which compIIdx diffuses
+    */
+   Scalar diffusionCoefficient(int compIIdx, int compJIdx = 0) const
+   {
+       if (compIIdx == compJIdx)
+           DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for compIIdx = compJIdx");
+       return diffCoefficient_[compIIdx][compJIdx];
+   }
+
+    /*!
+    * \brief Returns the effective diffusion coefficient \f$\mathrm{[m^2/s]}\f$
+    *
+    * \param compIIdx the index of the component which diffusive
+    * \param compJIdx the index of the component with respect to which compIIdx diffuses
+    */
+   Scalar effectiveDiffusivity(int compIIdx, int compJIdx = 0) const
+   {
+       return diffusionCoefficient(compIIdx, compJIdx);
+   }
 
-     /*!
-     * \brief Returns the effective diffusion coefficient \f$\mathrm{[m^2/s]}\f$
-     *
-     * \param compIIdx the index of the component which diffusive
-     * \param compJIdx the index of the component with respect to which compIIdx diffuses
-     */
-    Scalar effectiveDiffusivity(int compIIdx, int compJIdx = fluidSystemPhaseIdx) const
-    {
-        return diffusionCoefficient(compIIdx, compJIdx);
-    }
 
     /*!
      * \brief Return the fluid state of the control volume.
@@ -282,7 +280,7 @@ public:
 
 protected:
     FluidState fluidState_;
-    std::array<std::array<Scalar, numComponents>, numComponents> diffCoefficient_;
+    std::array<std::array<Scalar, numComponents>, numComponents> diffCoefficient_ = {};
 };
 
 } // end namespace Dumux
diff --git a/dumux/freeflow/compositional/vtkoutputfields.hh b/dumux/freeflow/compositional/vtkoutputfields.hh
index 78927b11bce0ab898c747b0ecb0c1225fe18f927..e0bbce332913632741dea8f53404eeee71e05a95 100644
--- a/dumux/freeflow/compositional/vtkoutputfields.hh
+++ b/dumux/freeflow/compositional/vtkoutputfields.hh
@@ -56,12 +56,13 @@ public:
             vtk.addVolumeVariable([j](const auto& v){ return v.moleFraction(j); }, "x^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(phaseIdx));
             if (j != phaseIdx)
             {
-                vtk.addVolumeVariable([j](const auto& v){ return v.diffusionCoefficient(j); }, "D^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(phaseIdx));
+                vtk.addVolumeVariable([j](const auto& v){ return v.diffusionCoefficient(0, j); }, "D^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(phaseIdx));
+
+                if (ModelTraits::usesTurbulenceModel())
+                // the eddy diffusivity is recalculated for an arbitrary component which is not the phase component
+                vtk.addVolumeVariable([j](const auto& v){ return v.effectiveDiffusivity(0, j) - v.diffusionCoefficient(0, j); }, "D_t");
             }
         }
-        if (ModelTraits::usesTurbulenceModel())
-            // the eddy diffusivity is recalculated for an arbitrary component which is not the phase component
-            vtk.addVolumeVariable([](const auto& v){ return v.effectiveDiffusivity(1-phaseIdx) - v.diffusionCoefficient(1-phaseIdx); }, "D_t");
     }
 };
 
diff --git a/dumux/freeflow/navierstokes/volumevariables.hh b/dumux/freeflow/navierstokes/volumevariables.hh
index 5cd11521ce91535bf0d72cab1062524d91b4a8f8..f4391b1733fc2758fb22ef72928ce34bd54784c5 100644
--- a/dumux/freeflow/navierstokes/volumevariables.hh
+++ b/dumux/freeflow/navierstokes/volumevariables.hh
@@ -41,8 +41,6 @@ class NavierStokesVolumeVariables : public FreeFlowVolumeVariables< Traits, Navi
 
     using Scalar = typename Traits::PrimaryVariables::value_type;
 
-    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
-
 public:
     //! export the underlying fluid system
     using FluidSystem = typename Traits::FluidSystem;
@@ -83,43 +81,43 @@ public:
         const Scalar t = ParentType::temperature(elemSol, problem, element, scv);
         fluidState.setTemperature(t);
 
-        fluidState.setPressure(fluidSystemPhaseIdx, elemSol[0][Indices::pressureIdx]);
+        fluidState.setPressure(0, elemSol[0][Indices::pressureIdx]);
 
         // saturation in a single phase is always 1 and thus redundant
         // to set. But since we use the fluid state shared by the
         // immiscible multi-phase models, so we have to set it here...
-        fluidState.setSaturation(fluidSystemPhaseIdx, 1.0);
+        fluidState.setSaturation(0, 1.0);
 
         typename FluidSystem::ParameterCache paramCache;
-        paramCache.updatePhase(fluidState, fluidSystemPhaseIdx);
+        paramCache.updateAll(fluidState);
 
-        Scalar value = FluidSystem::density(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setDensity(fluidSystemPhaseIdx, value);
+        Scalar value = FluidSystem::density(fluidState, paramCache, 0);
+        fluidState.setDensity(0, value);
 
-        value = FluidSystem::molarDensity(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setMolarDensity(fluidSystemPhaseIdx, value);
+        value = FluidSystem::molarDensity(fluidState, paramCache, 0);
+        fluidState.setMolarDensity(0, value);
 
-        value = FluidSystem::viscosity(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setViscosity(fluidSystemPhaseIdx, value);
+        value = FluidSystem::viscosity(fluidState, paramCache, 0);
+        fluidState.setViscosity(0, value);
 
         // compute and set the enthalpy
         value = ParentType::enthalpy(fluidState, paramCache);
-        fluidState.setEnthalpy(fluidSystemPhaseIdx, value);
+        fluidState.setEnthalpy(0, value);
     }
 
     /*!
      * \brief Return the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
      *        the control volume.
      */
-    Scalar pressure(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.pressure(fluidSystemPhaseIdx); }
+    Scalar pressure(int phaseIdx = 0) const
+    { return fluidState_.pressure(0); }
 
     /*!
      * \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ of a given phase within the
      *        control volume.
      */
-    Scalar density(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.density(fluidSystemPhaseIdx); }
+    Scalar density(int phaseIdx = 0) const
+    { return fluidState_.density(0); }
 
     /*!
      * \brief Return temperature \f$\mathrm{[K]}\f$ inside the sub-control volume.
@@ -135,26 +133,26 @@ public:
      * \brief Returns the mass density of a given phase within the
      *        control volume.
      */
-    Scalar molarDensity(int phaseIdx = fluidSystemPhaseIdx) const
+    Scalar molarDensity(int phaseIdx = 0) const
     {
-        return fluidState_.molarDensity(fluidSystemPhaseIdx);
+        return fluidState_.molarDensity(0);
     }
 
     /*!
      * \brief Returns the molar mass of a given phase within the
      *        control volume.
      */
-    Scalar molarMass(int phaseIdx = fluidSystemPhaseIdx) const
+    Scalar molarMass(int phaseIdx = 0) const
     {
-        return fluidState_.averageMolarMass(fluidSystemPhaseIdx);
+        return fluidState_.averageMolarMass(0);
     }
 
     /*!
      * \brief Return the dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
      *        control volume.
      */
-    Scalar viscosity(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.viscosity(fluidSystemPhaseIdx); }
+    Scalar viscosity(int phaseIdx = 0) const
+    { return fluidState_.viscosity(0); }
 
     /*!
      * \brief Return the effective dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
index 98bf14c686aff002488b1bb1d6e8aded562d0406..91fe3d55be0652977dee0776958d07f1450ea3ba 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
@@ -391,7 +391,7 @@ public:
                 continue;
 
             Scalar schmidtNumber = elemVolVars[scvf.insideScvIdx()].kinematicViscosity()
-                                   / elemVolVars[scvf.insideScvIdx()].diffusionCoefficient(compIdx);
+                                   / elemVolVars[scvf.insideScvIdx()].diffusionCoefficient(0, compIdx);
             Scalar moleToMassConversionFactor = useMoles
                                                 ? 1.0 : FluidSystem::molarMass(compIdx);
             wallFunctionFlux[compIdx] +=
diff --git a/dumux/freeflow/rans/volumevariables.hh b/dumux/freeflow/rans/volumevariables.hh
index 2920a851e8f39af2372dc396f91c245c0acd47a1..1461d93d73a4ece4399b6141b927a3edd61f8d7d 100644
--- a/dumux/freeflow/rans/volumevariables.hh
+++ b/dumux/freeflow/rans/volumevariables.hh
@@ -49,7 +49,6 @@ class RANSVolumeVariables
     using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
 
     static constexpr bool enableEnergyBalance = Traits::ModelTraits::enableEnergyBalance();
-    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
 
 public:
 
@@ -232,13 +231,13 @@ public:
     Scalar eddyThermalConductivity() const
     { return eddyThermalConductivity_; }
 
-     /*!
-     * \brief Returns the effective diffusion coefficient \f$\mathrm{[m^2/s]}\f$
-     *
-     * \param compIIdx the index of the component which diffusive
-     * \param compJIdx the index of the component with respect to which compIIdx diffuses
-     */
-    Scalar effectiveDiffusivity(int compIIdx, int compJIdx = fluidSystemPhaseIdx) const
+    /*!
+    * \brief Returns the effective diffusion coefficient \f$\mathrm{[m^2/s]}\f$
+    *
+    * \param compIIdx the index of the component which diffusive
+    * \param compJIdx the index of the component with respect to which compIIdx diffuses
+    */
+    Scalar effectiveDiffusivity(int compIIdx, int compJIdx) const
     {
         return NavierStokesParentType::diffusionCoefficient(compIIdx, compJIdx) + eddyDiffusivity();
     }
diff --git a/dumux/freeflow/volumevariables.hh b/dumux/freeflow/volumevariables.hh
index 91937670f2b6a17c0bffb83cdb3e4c405a772092..e60b4a697c7a115a7da152dec2b94c60113e4432 100644
--- a/dumux/freeflow/volumevariables.hh
+++ b/dumux/freeflow/volumevariables.hh
@@ -137,7 +137,6 @@ class FreeFlowVolumeVariablesImplementation<Traits, Impl, true>
 {
     using ParentType = FreeFlowVolumeVariablesImplementation<Traits, Impl, false>;
     using Scalar = typename Traits::PrimaryVariables::value_type;
-    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
 
 public:
     //! export the type used for the primary variables
@@ -171,28 +170,28 @@ public:
      * \brief Returns the total internal energy of a phase in the
      *        sub-control volume.
      */
-    Scalar internalEnergy(int phaseIdx = fluidSystemPhaseIdx) const
-    { return ParentType::asImp_().fluidState().internalEnergy(fluidSystemPhaseIdx); }
+    Scalar internalEnergy(int phaseIdx = 0) const
+    { return ParentType::asImp_().fluidState().internalEnergy(0); }
 
     /*!
      * \brief Returns the total enthalpy of a phase in the sub-control
      *        volume.
      */
-    Scalar enthalpy(int phaseIdx = fluidSystemPhaseIdx) const
-    { return ParentType::asImp_().fluidState().enthalpy(fluidSystemPhaseIdx); }
+    Scalar enthalpy(int phaseIdx = 0) const
+    { return ParentType::asImp_().fluidState().enthalpy(0); }
 
     /*!
      * \brief Returns the component enthalpy \f$\mathrm{[J/(kg*K)]}\f$ in the sub-control volume.
      */
     Scalar componentEnthalpy(unsigned int compIdx) const
-    { return FluidSystem::componentEnthalpy(ParentType::asImp_().fluidState(), fluidSystemPhaseIdx, compIdx); }
+    { return FluidSystem::componentEnthalpy(ParentType::asImp_().fluidState(), 0, compIdx); }
 
     /*!
      * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
      *        of the fluid phase in the sub-control volume.
      */
     Scalar thermalConductivity() const
-    { return FluidSystem::thermalConductivity(ParentType::asImp_().fluidState(), fluidSystemPhaseIdx); }
+    { return FluidSystem::thermalConductivity(ParentType::asImp_().fluidState(), 0); }
 
     /*!
      * \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
@@ -208,7 +207,7 @@ public:
      *        in the sub-control volume.
      */
     Scalar heatCapacity() const
-    { return FluidSystem::heatCapacity(ParentType::asImp_().fluidState(), fluidSystemPhaseIdx); }
+    { return FluidSystem::heatCapacity(ParentType::asImp_().fluidState(), 0); }
 
     //! The temperature is a primary variable for non-isothermal models
     template<class ElemSol, class Problem, class Element, class Scv>
@@ -226,7 +225,7 @@ public:
     static Scalar enthalpy(const FluidState& fluidState,
                            const ParameterCache& paramCache)
     {
-        return FluidSystem::enthalpy(fluidState, paramCache, fluidSystemPhaseIdx);
+        return FluidSystem::enthalpy(fluidState, paramCache, 0);
     }
 };
 }
diff --git a/test/freeflow/ransnc/flatplatetestproblem.hh b/test/freeflow/ransnc/flatplatetestproblem.hh
index 040bb514dbbd222e7fbcf9e4ba123a4dbf7c830a..723868298f52de7cab857c6f6a591eaf9286d0cb 100644
--- a/test/freeflow/ransnc/flatplatetestproblem.hh
+++ b/test/freeflow/ransnc/flatplatetestproblem.hh
@@ -27,6 +27,7 @@
 #include <dune/grid/yaspgrid.hh>
 
 #include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 #include <dumux/freeflow/turbulenceproperties.hh>
 
@@ -81,16 +82,16 @@ namespace Properties
   #endif
 #endif
 
-NEW_PROP_TAG(FluidSystem);
-
-// Select the fluid system
-SET_TYPE_PROP(FlatPlateNCTestTypeTag, FluidSystem,
-              FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-SET_INT_PROP(FlatPlateNCTestTypeTag, PhaseIdx,
-             GET_PROP_TYPE(TypeTag, FluidSystem)::gasPhaseIdx);
+// The fluid system
+SET_PROP(FlatPlateNCTestTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the air phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
-SET_INT_PROP(FlatPlateNCTestTypeTag, ReplaceCompEqIdx, GET_PROP_VALUE(TypeTag, PhaseIdx));
+// replace the main component balance eq with a total balance eq
+SET_INT_PROP(FlatPlateNCTestTypeTag, ReplaceCompEqIdx, 0);
 
 // Set the grid type
 SET_TYPE_PROP(FlatPlateNCTestTypeTag, Grid,
@@ -158,9 +159,8 @@ class FlatPlateNCTestProblem : public ZeroEqProblem<TypeTag>
     using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
 
     static constexpr auto dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
-    static const unsigned int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr auto transportEqIdx = Indices::conti0EqIdx;
-    static constexpr auto transportCompIdx = Indices::conti0EqIdx;
+    static constexpr auto transportEqIdx = Indices::conti0EqIdx + 1;
+    static constexpr auto transportCompIdx = Indices::conti0EqIdx + 1;
 
 public:
     FlatPlateNCTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
@@ -171,6 +171,7 @@ public:
         FluidSystem::init();
         Dumux::TurbulenceProperties<Scalar, dimWorld, true> turbulenceProperties;
         FluidState fluidState;
+        const auto phaseIdx = 0;
         fluidState.setPressure(phaseIdx, 1e5);
         fluidState.setTemperature(temperature());
         fluidState.setMassFraction(phaseIdx, phaseIdx, 1.0);
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/darcyproblem.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/darcyproblem.hh
index 2f16c830fa18853d5510ee20e21b4aa948ab8afb..f4d864dc0dcbda86f0408d097c90b6737d5869de 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/darcyproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/darcyproblem.hh
@@ -33,6 +33,7 @@
 
 #include "./../1pspatialparams.hh"
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 #include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
 
@@ -49,10 +50,12 @@ NEW_TYPE_TAG(DarcyOnePTwoCTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
 SET_TYPE_PROP(DarcyOnePTwoCTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
 
 // The fluid system
-SET_TYPE_PROP(DarcyOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-// Use water as phase
-SET_INT_PROP(DarcyOnePTwoCTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::liquidPhaseIdx);
+SET_PROP(DarcyOnePTwoCTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
 // Use moles
 SET_BOOL_PROP(DarcyOnePTwoCTypeTag, UseMoles, true);
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/stokesproblem.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/stokesproblem.hh
index ebc528deaeec1b3c580a27e59f4888bd055fad7e..a3627d489667487e5fc4963024fe83218da11051 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/stokesproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/stokesproblem.hh
@@ -26,6 +26,7 @@
 
 #include <dune/grid/yaspgrid.hh>
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 
 #include <dumux/freeflow/navierstokes/problem.hh>
@@ -41,8 +42,13 @@ namespace Properties
 {
 NEW_TYPE_TAG(StokesOnePTwoCTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
 
-// Set the fluid system
-SET_TYPE_PROP(StokesOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+// The fluid system
+SET_PROP(StokesOnePTwoCTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
 // Set the grid type
 SET_TYPE_PROP(StokesOnePTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/darcyproblem.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/darcyproblem.hh
index c07e21fac319d11cd0cad8dc4fcccd34ab444781..6fc7a42f49a39447a355d41104ba10f2cae61332 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/darcyproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/darcyproblem.hh
@@ -33,6 +33,7 @@
 
 #include "./../1pspatialparams.hh"
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 #include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
 
@@ -49,10 +50,12 @@ NEW_TYPE_TAG(DarcyOnePTwoCTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
 SET_TYPE_PROP(DarcyOnePTwoCTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
 
 // The fluid system
-SET_TYPE_PROP(DarcyOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-// Use water as phase
-SET_INT_PROP(DarcyOnePTwoCTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::liquidPhaseIdx);
+SET_PROP(DarcyOnePTwoCTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
 // Use moles
 SET_BOOL_PROP(DarcyOnePTwoCTypeTag, UseMoles, true);
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/stokesproblem.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/stokesproblem.hh
index 5d6fed0b9543207c88beb4ac1fa5763f23409cde..ade5e9b2e1f9d3c4168743c4297f0269d584eaf1 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/stokesproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/stokesproblem.hh
@@ -26,6 +26,7 @@
 
 #include <dune/grid/yaspgrid.hh>
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 
 #include <dumux/freeflow/navierstokes/problem.hh>
@@ -41,11 +42,13 @@ namespace Properties
 {
 NEW_TYPE_TAG(StokesOnePTwoCTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
 
-// Set the fluid system
-SET_TYPE_PROP(StokesOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-// Use water as phase
-SET_INT_PROP(StokesOnePTwoCTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::liquidPhaseIdx);
+// The fluid system
+SET_PROP(StokesOnePTwoCTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
 // Set the grid type
 SET_TYPE_PROP(StokesOnePTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/stokesproblem.hh b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/stokesproblem.hh
index 7e0f64317079d5ebc22fd7973df772686dde3fe8..dd402578eb847da3a6a81bca50d4eb1c9d68ff88 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/stokesproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/stokesproblem.hh
@@ -26,6 +26,7 @@
 
 #include <dune/grid/yaspgrid.hh>
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 
 #include <dumux/freeflow/navierstokes/problem.hh>
@@ -49,11 +50,13 @@ NEW_TYPE_TAG(StokesOnePTwoCTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, Navier
 // Set the grid type
 SET_TYPE_PROP(StokesOnePTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
 
-// set the fluid system
-SET_TYPE_PROP(StokesOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-// set phase index (air)
-SET_INT_PROP(StokesOnePTwoCTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::gasPhaseIdx);
+// The fluid system
+SET_PROP(StokesOnePTwoCTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the water phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
 SET_INT_PROP(StokesOnePTwoCTypeTag, ReplaceCompEqIdx, 3);
 
@@ -107,9 +110,6 @@ class StokesSubProblem : public NavierStokesProblem<TypeTag>
 
     static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
 
-    static constexpr auto phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr auto transportCompIdx = 1 - phaseIdx;
-
 public:
     StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
     : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
@@ -128,10 +128,6 @@ public:
      */
     // \{
 
-
-    bool shouldWriteRestartFile() const
-    { return false; }
-
    /*!
      * \brief Return the temperature within the domain in [K].
      */
@@ -178,18 +174,10 @@ public:
             values.setNeumann(Indices::conti0EqIdx + 1);
         }
 
-        if (onUpperBoundary_(globalPos))
-        {
-            values.setDirichlet(Indices::velocityXIdx);
-            values.setDirichlet(Indices::velocityYIdx);
-            values.setNeumann(Indices::conti0EqIdx);
-            values.setNeumann(Indices::conti0EqIdx + 1);
-        }
-
         if (onRightBoundary_(globalPos))
         {
             values.setDirichlet(Indices::pressureIdx);
-            values.setOutflow(Indices::conti0EqIdx);
+            values.setOutflow(Indices::conti0EqIdx + 1);
 
 #if NONISOTHERMAL
             values.setOutflow(Indices::energyBalanceIdx);
@@ -242,17 +230,17 @@ public:
         FluidState fluidState;
         updateFluidStateForBC_(fluidState, elemVolVars[scv].pressure());
 
-        const Scalar density = useMoles ? fluidState.molarDensity(phaseIdx) : fluidState.density(phaseIdx);
+        const Scalar density = useMoles ? fluidState.molarDensity(0) : fluidState.density(0);
         const Scalar xVelocity = xVelocity_(globalPos);
 
         if (onLeftBoundary_(globalPos))
         {
             // rho*v*X at inflow
-            values[Indices::conti0EqIdx] = -xVelocity * density * refMoleFrac();
-            values[Indices::conti0EqIdx + 1] = -xVelocity * density * (1.0 - refMoleFrac());
+            values[Indices::conti0EqIdx + 1] = -xVelocity * density * refMoleFrac();
+            values[Indices::conti0EqIdx] = -xVelocity * density * (1.0 - refMoleFrac());
 
 #if NONISOTHERMAL
-            values[Indices::energyBalanceIdx] = -xVelocity * fluidState.density(phaseIdx) * fluidState.enthalpy(phaseIdx);
+            values[Indices::energyBalanceIdx] = -xVelocity * fluidState.density(0) * fluidState.enthalpy(0);
 #endif
         }
 
@@ -301,11 +289,11 @@ public:
         FluidState fluidState;
         updateFluidStateForBC_(fluidState, refPressure());
 
-        const Scalar density = FluidSystem::density(fluidState, phaseIdx);
+        const Scalar density = FluidSystem::density(fluidState, 0);
 
         PrimaryVariables values(0.0);
         values[Indices::pressureIdx] = refPressure() + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
-        values[Indices::conti0EqIdx] = refMoleFrac();
+        values[Indices::conti0EqIdx + 1] = refMoleFrac();
         values[Indices::velocityXIdx] = xVelocity_(globalPos);
 
 #if NONISOTHERMAL
@@ -373,22 +361,22 @@ private:
     void updateFluidStateForBC_(FluidState& fluidState, const Scalar pressure) const
     {
         fluidState.setTemperature(refTemperature());
-        fluidState.setPressure(phaseIdx, pressure);
-        fluidState.setSaturation(phaseIdx, 1.0);
-        fluidState.setMoleFraction(phaseIdx, transportCompIdx, refMoleFrac());
-        fluidState.setMoleFraction(phaseIdx, phaseIdx, 1.0 - refMoleFrac());
+        fluidState.setPressure(0, pressure);
+        fluidState.setSaturation(0, 1.0);
+        fluidState.setMoleFraction(0, 1, refMoleFrac());
+        fluidState.setMoleFraction(0, 0, 1.0 - refMoleFrac());
 
         typename FluidSystem::ParameterCache paramCache;
-        paramCache.updatePhase(fluidState, phaseIdx);
+        paramCache.updatePhase(fluidState, 0);
 
-        const Scalar density = FluidSystem::density(fluidState, paramCache, phaseIdx);
-        fluidState.setDensity(phaseIdx, density);
+        const Scalar density = FluidSystem::density(fluidState, paramCache, 0);
+        fluidState.setDensity(0, density);
 
-        const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
-        fluidState.setMolarDensity(phaseIdx, molarDensity);
+        const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0);
+        fluidState.setMolarDensity(0, molarDensity);
 
-        const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
-        fluidState.setEnthalpy(phaseIdx, enthalpy);
+        const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
+        fluidState.setEnthalpy(0, enthalpy);
     }
 
     //! \brief set the profile of the inflow velocity (horizontal direction)
diff --git a/test/multidomain/boundary/stokesdarcy/1p_2p/stokesproblem.hh b/test/multidomain/boundary/stokesdarcy/1p_2p/stokesproblem.hh
index df3a7290e5e3c24fce2f646d6706fe4b515785e2..f4264647660c22e656d054b04be3066b15125d62 100644
--- a/test/multidomain/boundary/stokesdarcy/1p_2p/stokesproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p_2p/stokesproblem.hh
@@ -26,6 +26,7 @@
 
 #include <dune/grid/yaspgrid.hh>
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/2pimmiscible.hh>
 #include <dumux/material/fluidsystems/1pliquid.hh>
 #include <dumux/material/fluidsystems/1pgas.hh>
@@ -53,16 +54,16 @@ SET_TYPE_PROP(StokesOnePTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffset
 SET_PROP(StokesOnePTypeTag, FluidSystem)
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using H2OType = Dumux::Components::SimpleH2O<Scalar>;
-    using H2OPhase = Dumux::FluidSystems::OnePLiquid<Scalar, H2OType>;
-    using AirType = Dumux::Components::TabulatedComponent<Components::Air<Scalar>, false >;
-    using AirPhase = Dumux::FluidSystems::OnePGas<Scalar, AirType>;
-    using type = FluidSystems::TwoPImmiscible<Scalar, H2OPhase, AirPhase> ;
+    using H2OType = Components::SimpleH2O<Scalar>;
+    using H2OPhase = FluidSystems::OnePLiquid<Scalar, H2OType>;
+    using AirType = Components::TabulatedComponent<Components::Air<Scalar>, false >;
+    using AirPhase = FluidSystems::OnePGas<Scalar, AirType>;
+    using TwoPFluidSystem = FluidSystems::TwoPImmiscible<Scalar, H2OPhase, AirPhase> ;
+
+    static constexpr auto phaseIdx = 1; // simulate the air phase
+    using type = FluidSystems::OnePAdapter<TwoPFluidSystem, phaseIdx>;
 };
 
-// consider air as phase
-SET_INT_PROP(StokesOnePTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::phase1Idx);
-
 // Set the problem property
 SET_TYPE_PROP(StokesOnePTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );